Google地图的dismo::gmap()和ggplot2

Question

我使用函数 dismo::gmap() 获取到一张地图，并想用 ggplot2 绘制它，因为我想使用 geom_point 和其他 ggplot 函数添加不同的特征。我更喜欢使用 dismo::gmap，而不是 ggmap::get_map() 下载谷歌地图图层。这是因为 dismo::gmap() 返回一个包含完整 CRS 信息的栅格图层（来自 raster 包），因此应该可以修改图层的投影。

> 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相关联的值

- jl-blancopastor

---

尝试使用 ggfortify::fortify(mapImageData2) 生成一个数据框，然后可以用 ggplot 绘制。 - Richard Telford

Answer 1

这个问题有两个方面。一个是如何让ggplot2绘制Raster*对象，另一个是如何重投影栅格并保留其值。

原始问题中包含了以下代码：

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 <- 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)

现在我们看到了重新投影后的地图。进展！但是我们仍然无法使用ggplot2绘制mapImageData2_proj，因为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))

- Jacob Socolar

Answer 2

要在ggplot中直接绘制一个Raster*对象，您可以按照以下步骤使用库RasterVis：

r <- raster(system.file("external/test.grd", package="raster"))
s <- stack(r, r*2)
names(s) <- c('meuse', 'meuse x 2')

library(ggplot2)

theme_set(theme_bw())
gplot(s) + geom_tile(aes(fill = value)) +
          facet_wrap(~ variable) +
          scale_fill_gradient(low = 'white', high = 'blue') +
          coord_equal()

如您所见，我们使用的是 gplot 而不是 ggplot。这个函数适用于 Raster* 对象，因为它允许仅加载 Raster* 对象的一部分到 RAM 中。您甚至可以使用 gplot(s, maxpixels=50000) 来选择在地图上加载的最大像素数。
实际上，我建议不将您的 Raster 转换为点，因为如果您的栅格非常大，您将无法将其加载到 RAM 中...

- Sébastien Rochette