来自dismo :: gmap()和ggplot2的谷歌地图

Question

我获得了一个带有该功能的地图,dismo::gmap()并希望用ggplot2绘制它,因为我想使用geom_point和其他ggplot函数添加不同的功能.我更喜欢使用dismo::gmap而不是ggmap::get_map()下载谷歌地图图层.这是因为dismo::gmap(),与ggmap::get_map()包括完整CRS信息的包栅格返回栅格图层不同,因此应该可以修改图层的投影.

> head(data_info$latitude, 20)
#[1] 49.11306 49.39333 48.78083 51.85000 53.57361 50.67806 52.69083 52.21389 53.46361 50.99917 53.99750 53.54528 53.61417 48.00556 48.01306 53.45000
#[17] 51.93667 54.53083 51.95500 54.29639
> head(data_info$longitude, 20)
#[1] 13.134722 12.323056 13.803889 12.177778 14.143611 13.175833 12.649444 13.454167 11.629722 10.906111 11.415556  8.426944  7.160000 11.123889 10.786111
#[16] 12.766667 11.987222 13.091389 10.967500 13.684167


   e = extent(-14 , 58 , 28 , 64)
mapImageData2 <- gmap(e, type = c("terrain"), lonlat = TRUE, 
                      path = "&style=feature:all|element:labels|visibility:off&style=feature:administrative.country|element:geometry.stroke|visibility:off")

mapImageData2_proj <- projectExtent(mapImageData2, crs = "+proj=utm +zone=31 +datum=WGS84")

# plot the points on the map 
ggplot(mapImageData2_proj, extent = "device") + 
  geom_point(inherit.aes = FALSE, aes(x = data_info$longitude, y = data_info$latitude),
             data = gps, colour = "red", size = 1, pch = 20)

尝试此操作后,我收到以下错误:

错误:ggplot2不知道如何处理类RasterLayer的数据

如果我试试这个

plot(mapImageData2_proj)

.plotraster2中的错误(x,col = col,maxpixels = maxpixels,add = add,:没有与此RasterLayer关联的值

Answer 1

这个问题有两个问题.一个是如何让ggplot2绘制一个Raster*对象.另一种是如何在保留其值的同时重新投影栅格.

OP包含代码

library(dismo)
e = extent(-14 , 58 , 28 , 64)
mapImageData2 <- gmap(e, type = c("terrain"), lonlat = TRUE, 
                      path = "&style=feature:all|element:labels|visibility:off&style=feature:administrative.country|element:geometry.stroke|visibility:off")

如果我们运行这个,然后plot(mapImageData2)我们得到一个很好的情节.然后OP运行

mapImageData2_proj <- projectExtent(mapImageData2, crs = "+proj=utm +zone=31 +datum=WGS84")

现在,如果我们运行,plot(mapImageData2_proj)我们会收到错误!那是因为projectExtent返回没有值的RasterLayer.我们需要改为使用projectRaster().详情?projectExtent请见.所以我们运行:

mapImageData2_proj <- projectRaster(mapImageData2, crs = "+proj=utm +zone=31 +datum=WGS84")
plot(mapImageData2_proj)

现在我们看到重新投影的地图.进展!但是我们仍然无法mapImageData2_proj使用绘图ggplot2,因为ggplot2不知道如何处理Raster*对象.我们需要将栅格转换为数据帧.有几种方法可以做到这一点,但没有加载任何额外的包,一个很好的选择raster::rasterToPoints().例如:

myPoints <- raster::rasterToPoints(myRaster)
myDataFrame <- data.frame(myPoints)
colnames(myDataFrame) <- c("Longitude", "Latitude", "Values")

ggplot(data=myDataFrame, aes_string(y = "Latitude", x = "Longitude")) + 
   geom_raster(aes(fill = Values))

所以在OP的例子中将它们放在一起:

library(dismo)
e = extent(-14 , 58 , 28 , 64)
mapImageData2 <- gmap(e, type = c("terrain"), lonlat = TRUE, 
                      path = "&style=feature:all|element:labels|visibility:off&style=feature:administrative.country|element:geometry.stroke|visibility:off")

plot(mapImageData2)

mapImageData2_proj <- projectRaster(mapImageData2, crs = "+proj=utm +zone=31 +datum=WGS84")

plot(mapImageData2_proj)

myRaster <- mapImageData2_proj
myPoints <- raster::rasterToPoints(myRaster)
myDataFrame <- data.frame(myPoints)
colnames(myDataFrame) <- c("X", "Y", "Values")

ggplot(data=myDataFrame, aes_string(y = "Y", x = "X")) + 
  geom_raster(aes(fill = Values))