将 shapefile 添加到 ggmap

Question

与此处描述的问题有些类似，我在对齐 shapefile 和 ggmap 对象时遇到问题。

1. 我的形状文件由澳大利亚维多利亚州的局部区域边界组成，我试图将它们覆盖在该州（维多利亚州）的谷歌地图上。

   源 shapefile 具有以下 PROJ4 字符串（从 prj 文件中提取）

   [+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs]
   

   对应于EPSG:4283。

   这是我的 shapefile 对象的摘要sp：

   > summary(sp)
   Object of class SpatialPolygonDataFrame
   Coordinates:
           min        max
   x  96.81677 159.109219
   y -43.74051  -9.142176
   Is projected: FALSE
   proj4string:
   [+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs]
   

2. 我转换 shapefile 以匹配谷歌的伪墨卡托投影（至少我认为我正在做的）

   sp <- spTransform(sp, CRS("+proj=longlat +init=epsg:3857"))
   

   并变成sp一个强化的数据框df.sp。

3. 然后我用

   map <- get_map("Victoria", zoom = 7, maptype = "terrain", source = "google")
   

   获取维多利亚的谷歌地形图，并（gg）绘制它们

   ggmap(map) + geom_polygon(data = df.sp, aes(x = long, y = lat, group = group)) + coord_equal() + theme_map()
   

从结果图中可以清楚地看出，shapefile 坐标和 googlemap 坐标不重叠。我的坐标变换有问题吗？如何正确匹配 shapefile 和 googlemap 坐标？我将不胜感激对此事的任何帮助/见解。

Answer 1

这是我的第一篇文章，所以如果样式不太好，请原谅我。我发现无论出于何种原因，包ESPG:3857中的默认设置rgdal都没有在spTransform函数中正确转换投影。

• 使用gene在GIS StackExchange中的精彩答案，我被引导到该Prj2EPSG网站和空间参考网站。我仔细检查了 ABS 的proj4字符串（我假设这就是该数据的来源），正如您所说（即EPSG:4283）。

• 当在空间参考中输入 EPSG:3857 时，我得到了这个字符串（使用 Proj4 链接）+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs。

• 我替换+proj=merc为+proj=longlat投影经度和纬度坐标。

• 使用该字符串merc替换似乎比仅使用该字符串给我带来更好的结果+init=epsg:3857，这向我表明 rgdal 包处理它的方式存在一些差异。

我很想知道进一步改进的方法，但我希望这可以帮助您同时获得更好的结果。

编辑：我还发现，在绘制时使用ggplot2添加coord_map()似乎可以修复形状文件在近距离时的对齐方式，但在进一步缩小时会使它们稍微偏离。在变焦比例 7 及以上时，我似乎获得了完美的对齐，但在 6 及以下时却没有。似乎无论出于何种原因，不同的缩放级别都会稍微改变投影。

所以尝试一下

ggmap(map) + 
  geom_polygon(data = df.sp, aes(x = long, y = lat, group = group)) +
  coord_equal() + 
  theme_map() + 
  coord_map()