mapa google z dismo :: gmap() i ggplot2

Question

Otrzymałem mapę z funkcją dismo::gmap() i chcę ją narysować za pomocą ggplot2, ponieważ chcę dodać różne feautures za pomocą geom_point i innych funkcji ggplota. Wolę używać dismo::gmap zamiast ggmap::get_map(), aby pobrać warstwę mapy google. Dzieje się tak, ponieważ dismo::gmap(), w przeciwieństwie do ggmap::get_map(), zwraca warstwę rastrową z rastra pakietu zawierającego kompletne informacje CRS i dlatego powinna istnieć możliwość modyfikacji rzutu warstwy.

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

Po wypróbowaniu tego, pojawia się następujący błąd:

Error: ggplot2 doesn't know how to deal with data of class RasterLayer

Gdy próbuję to

plot(mapImageData2_proj) 

Error in .plotraster2(x, col = col, maxpixels = maxpixels, add = add, : no values associated with this RasterLayer

Answer 1

Istnieją dwa problemy w tej kwestii. Jednym z nich jest uzyskanie ggplot2 do wykreślenia obiektu Raster *. Drugi polega na tym, jak zresetować raster, zachowując jego wartości.

PO zawiera kod

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

Jeśli uruchomić to i wtedy zrobić plot(mapImageData2) możemy uzyskać ładne działki. PO następnie uruchamia

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

Teraz, jeśli prowadzimy plot(mapImageData2_proj) otrzymujemy błąd! To dlatego, że projectExtent zwraca RasterLayer bez wartości. Musimy zamiast tego użyć projectRaster(). Szczegóły: ?projectExtent. Tak więc uruchamiamy:

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

Teraz widzimy reprojektowaną mapę. Postęp! Ale nadal nie możemy wykreślić mapImageData2_proj używając ggplot2, ponieważ ggplot2 nie wie, jak obsługiwać obiekt Raster *. Musimy przekształcić nasz raster w ramkę danych. Można to zrobić na kilka sposobów, ale bez ładowania dodatkowych pakietów dobrym rozwiązaniem jest raster::rasterToPoints(). Na przykład:

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

więc umieścić je wszystkie razem w przykładzie PO za:

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

Answer 2

Aby bezpośrednio wykreślić Raster * Obiekt w ggplot, można korzystać z biblioteki RasterVis następująco:

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

Jak widać, używamy gplot, a nie ggplot. Ta funkcja jest przeznaczona dla obiektów Raster *, ponieważ umożliwia załadowanie tylko części obiektu Raster * w pamięci RAM. Możesz nawet wybrać maksymalną liczbę pikseli do załadowania na mapę za pomocą gplot(s, maxpixels=50000).
Rzeczywiście, nie zalecam przekształcania Rastra jako punktów, ponieważ, jeśli twój raster jest ogromny, nie będziesz mógł załadować go w RAM ...