r - 将 shapefile 添加到 ggmap

Question

有点类似于描述的问题 here我在对齐 shapefile 和 ggmap 对象时遇到问题。

1. 我的 shapefile 由澳大利亚维多利亚州的局部区域边界组成，我试图将它们覆盖在该州(维多利亚州)的谷歌地图之上。

   源 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

这是我的第一篇文章，如果样式不尽如人意，请原谅我。我发现无论出于何种原因，rgdal 包中的默认 ESPG:3857 并未在 spTransform 函数中正确转换投影。

• 使用 gene's wonderful answer在 GIS StackExchange 中，我被定向到 Prj2EPSG网站和Spatial Reference网站。我仔细检查了 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()