Cannot plot sf layer on top of ggmap

Question

Subject of some previous questions - e.g. this. I cannot get the sf layer to render over the ggmap layer.

require(ggmap)
require(raster)
require(sf)

target_bbox <- c(left = -0.65, bottom = 51.2, right = 0.45, top = 51.8)
map = get_stamenmap(target_bbox, zoom = 9, maptype = 'toner-2010')

# can test this with: ggmap(map)

e = extent(list(x = c(-0.65, 0.45), y = c(51.2, 51.8)))
r = raster(volcano)
crs(r) = CRS("+proj=longlat +datum=WGS84")
p = raster::rasterToPolygons(r)
sf = st_as_sf(p) %>% st_transform(3857)

# test with: ggplot() + geom_sf(data = sf, aes(fill = layer), col = NA)

ggmap(map) + 
   geom_sf(data = sf, aes(fill = layer), col = NA, alpha = 0.4, inherit.aes = FALSE)
Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Coordinate system clash but no errors, yet no layer is visible. Any suggestions much appreciated.

This github solution - extracting the bounding box from the ggmap object with custom function ggmap_bbox() - only succeeds in adding another warning to that above:

Warning message:
In st_is_longlat(x) :
  bounding box has potentially an invalid value range for longlat data

score 1 · Accepted Answer · answered Dec 24 '20 at 21:29

It looks like you're trying to plot the sf volcano object over London, but you've defined it to be at lon/lat 0,0. Other than that, your call to ggmap and then geom_sf works fine. The problem was just that the sf object wasn't in the bounding box of the ggmap.

I had to change it to a 'Spatial' object to move it using maptools::elide, and then back to an 'sf' object as I couldn't find an easy way to shift sf polygons. The shift is an eyeball approximation.

sf_moved <- sf %>% 
  st_transform(3857) %>% # <- may be unnecessary
  as('Spatial') %>% 
  maptools::elide(.,shift = c(-80000, 6650000)) %>% 
  as('sf') %>% 
  st_set_crs(3857) %>% 
  st_transform(4326) 

ggmap(map) + 
  geom_sf(data = sf_moved, inherit.aes = FALSE, 
          col = NA, aes(fill = layer), alpha = .4) + 
  scale_fill_viridis_c()

Oh dear. Rookie error - I forgot to assign the extent to the raster object `extent(r) = e` — geotheory, Dec 24 '20 at 21:33
Happens to everyone. Assigning the extent would've been much easier than my method! — mrhellmann, Dec 24 '20 at 21:41

Cannot plot sf layer on top of ggmap

1 Answers1