4

我正在尝试将 ESRI shapefile (.shp) 添加到北卡罗来纳州的 ggmap 图，该图具有以下代码：

x <- get_map(location="North Carolina", zoom=6, maptype="terrain")

ggmap(x) + geom_polygon(data=citylim83.df, aes(x=long, y=lat), fill="red", alpha=0.2)

我已经加载并强化到citylim83.df 不显示. 这是用于将 shapefile 加载到 ggplot 中的代码：

    epsgs <- make_EPSG()
    citylim <- readOGR(dsn=".", layer="MunicipalBoundaries_polys")`

在 EPSG 中进行搜索后，MunicipalBoundaries 的投影单位是 ft-US 用于国家平面系统。即使这个 .shp 的地理坐标系为 NAD83，我也想将其投影到 NAD83 以摆脱国家平面系统（NAD83 (UTM-17N) 的 EPSG 代码为 26917）：

    citylim83 <- spTransform(citylim, CRS("+init=epsg:26917"))
    summary(citylim83)

    Object of class SpatialPolygonsDataFrame
    Is projected: TRUE 
    [+init=epsg:26917 +proj=utm +zone=17 +ellps=GRS80 +datum=NAD83 +units=m

    citylim83.df <- fortify(citylim83)

然后在上面显示的 ggmap 代码中使用了这个数据框。

即使它现在肯定会投影到 NAD83 中，它仍然不会显示在基础 ggmap 上。我导入的 get_map 对象的基本投影是什么？是否有一个命令可以找出这个，以便我可以将我的地图与我想要显示在它上面的 shapefile 相匹配？我是否必须“取消投影”我的 citylim 对象？仅供参考，如果不清楚，shapefile 是北卡罗来纳州每个城市的城市边界。任何帮助将不胜感激，因为我对 ggplot2/ggmap 社区非常陌生。

4

1 回答 1

7

在这里你有一个答案：

library(ggmap)
library(maptools)
shapefile <- readShapeSpatial('MunicipalBoundaries_polys.shp', proj4string = CRS("+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.2199999997 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"))
shp <- spTransform(shapefile, CRS("+proj=longlat +datum=WGS84"))
data <- fortify(shp)

nc <- get_map("North Carolina", zoom = 6, maptype = 'terrain')
ncmap <- ggmap(nc,  extent = "device")
ncmap +
  geom_polygon(aes(x = long, y = lat, group = group), data = data,
               colour = 'grey', fill = 'black', alpha = .4, size = .1)

坐标需要转换为 longlat 投影并且它正在工作。

于 2013-08-06T17:08:05.320 回答