将 shapefile 添加到 ggmap

Question

与描述的问题有些相似here我在对齐 shapefile 和 ggmap 对象时遇到问题。

我的 shapefile 由澳大利亚维多利亚州的局部区域边界组成，我试图将它们覆盖在州（维多利亚）的 google 地图之上。

源 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]

我转换 shapefile 以匹配 google 的伪墨卡托投影（至少我认为我正在这样做）
```
sp <- spTransform(sp, CRS("+proj=longlat +init=epsg:3857"))
```
并将sp变成一个强化数据框df.sp。

然后我用

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

获取维多利亚的 google 地形图，并 (gg) 绘制它们

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

从结果图中可以清楚地看出，shapefile 坐标和 google地图坐标不重叠。我在坐标转换方面做错了吗？如何正确匹配 shapefile 和 googlemap 坐标？我将不胜感激 help/insight 参与此事。

Answer 1

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

使用 gene's wonderful answer in GIS StackExchange, I was directed to the Prj2EPSG website and the Spatial Reference 网站。我仔细检查了 ABS 的 proj4 字符串（我假设这是该数据的来源）并且正如您所说的那样（即 EPSG:4283）。
在空间参考中输入 EPSG:3857 时，我得到了这个字符串（使用 Proj4 link）+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=longlat 替换了 +proj=merc 以投影经纬度坐标。
使用替换了 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()

将 shapefile 添加到 ggmap

Add shapefile to ggmap

gis

r

ggplot2

ggmap