How to put a geom_sf produced map on top of a ggmap produced raster

Question

I tried the following code:

library(ggplot2)
library(ggmap)
library(sf)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))
str(nc)

 Classes ‘sf’ and 'data.frame':  100 obs. of  15 variables:
 $ AREA     : num  0.114 0.061 0.143 0.07 0.153 0.097 0.062 0.091 0.118 0.124 ...
 $ PERIMETER: num  1.44 1.23 1.63 2.97 2.21 ...
 $ CNTY_    : num  1825 1827 1828 1831 1832 ...
 $ CNTY_ID  : num  1825 1827 1828 1831 1832 ...
 $ NAME     : Factor w/ 100 levels "Alamance","Alexander",..: 5 3 86 27 66 46 15 37 93 85 ...
 $ FIPS     : Factor w/ 100 levels "37001","37003",..: 5 3 86 27 66 46 15 37 93 85 ...
 $ FIPSNO   : num  37009 37005 37171 37053 37131 ...
 $ CRESS_ID : int  5 3 86 27 66 46 15 37 93 85 ...
 $ BIR74    : num  1091 487 3188 508 1421 ...
 $ SID74    : num  1 0 5 1 9 7 0 0 4 1 ...
 $ NWBIR74  : num  10 10 208 123 1066 ...
 $ BIR79    : num  1364 542 3616 830 1606 ...
 $ SID79    : num  0 3 6 2 3 5 2 2 2 5 ...
 $ NWBIR79  : num  19 12 260 145 1197 ...
 $ geometry :sfc_MULTIPOLYGON of length 100; first list element: List of 1
 ..$ :List of 1
 .. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ...
   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
   - attr(*, "sf_column")= chr "geometry"
   - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
   ..- attr(*, "names")= chr  "AREA" "PERIMETER" "CNTY_" "CNTY_ID" ...

map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google")
a <- unlist(attr(map,"bb")[1, ])
bb <- st_bbox(nc)
ggplot() + 
  annotation_raster(map, xmin = a[2], xmax = a[4], ymin = a[1], ymax = a[3]) + 
  xlim(c(bb[1], bb[3])) + ylim(c(bb[2], bb[4])) + 
  geom_sf(data = nc, aes(fill = AREA))

The two layers do not overlap properly; I tried changing projection with coord_sf() but I did not have success.

any suggestion? thanks

Answer 1

I've struggled with this myself, and with the help of this post I've come up with a solution. The bounding box of the ggmap object is in WGS84 (EPSG:4326), but the actual raster is in EPSG:3857. You have to hack the bounding box of the ggmap object to be in the same CRS as the underlying data:

library(ggplot2)
library(ggmap)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.3.0, proj.4 5.1.0

nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

# Transform nc to EPSG 3857 (Pseudo-Mercator, what Google uses)
nc_3857 <- st_transform(nc, 3857)

map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google")

# Define a function to fix the bbox to be in EPSG:3857
ggmap_bbox <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))

  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))

  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}

# Use the function:
map <- ggmap_bbox(map)

ggmap(map) + 
  coord_sf(crs = st_crs(3857)) + # force the ggplot2 map to be in 3857
  geom_sf(data = nc_3857, aes(fill = AREA), inherit.aes = FALSE)

Created on 2018-06-13 by the reprex package (v0.2.0).

Answer 2

You can use the plotting method from the sf package to do this. You'll need to specify the coordinate reference system, which we will need to assume (and it looks like correctly so) is 3857.

library(ggplot2)
library(ggmap)
library(sf)

nc_shp <- st_read(system.file("shape/nc.shp", package = "sf"))
nc_map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google")

st_crs(nc_map)
# Coordinate Reference System: NA

# assume the coordinate refence system is 3857
plot(st_transform(nc_shp, crs = 3857)[1], bgMap = nc_map)