Search code examples
rgisr-sfggmapstamen-maps

How to properly combine a polygon map on top of a raster map in r?

I would like to combine a shapefile map of some US states on top of a raster map obtained from using library ggmap.

This is the code I have

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

us.state.longlat <- us_states()
us.state.longlat.cropped.large <- sf::st_crop(us.state.longlat, 
                                              c(xmin=-100, xmax=-80, ymin=30, ymax=50))

bbox <- setNames(st_bbox(us.state.longlat.cropped.large), 
                 c("left", "bottom", "right", "top"))

ggmap(get_map(bbox, source = "stamen", zoom = 5)) +
  geom_sf(data = us.state.longlat.cropped.large,
          inherit.aes = FALSE,
          color = "red", 
          fill = "transparent")

This is what I got:

enter image description here

Also got the following warning that I am not sure what it means:

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

As you see there is a mismatch between the US state shapefile and the raster map. The coordinate reference system for the shapefile map us.state.longlat.cropped.large is this:

st_crs(us.state.longlat.cropped.large)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

The issue seems to be less worrisome when cropping a small area map. For example:

us.state.longlat.cropped.small <- sf::st_crop(us.state.longlat, c(xmin=-92, xmax=-85, ymin=32, ymax=36))
bbox <- setNames(st_bbox(us.state.longlat.cropped.small), 
                 c("left", "bottom", "right", "top"))

ggmap(get_map(bbox, source = "stamen", zoom = 7)) +
  geom_sf(data = us.state.longlat.cropped.small,
          inherit.aes = FALSE,
          color = "red", 
          fill = "transparent")

This is what I got for a reduced area of the map: enter image description here

In any case, I don't know what I am missing here but in both cases state borders should match between the shapefile and the raster map.

Solution

  • Probably this is due to different CRS on the tile (most likely EPSG:3857) and the shapefile (probably EPSG:4326).

    I would recommend an approach with two newer packages:

    • maptiles allows to download map tiles as true SpatRaster objects, so they can be modified accordingly.
    • tidyterra (shameless plug since I am the developer) allows to plot SpatRaster objects on the ggplot2 system. Also tidyterra is compatible with ggplot2::coord_sf(), meaning that the SpatRaster and the sf object would be internally projected to the same CRS thus avoding the inconsistencies your are facing with ggmap

    An additional word of caution: almost all the WMTS services (tiles, as OSM, Stamen, GoogleMaps, etc) provides the tiles on EPSG:3857 (WebMercator). What maptiles does is returning the tile on the CRS of the object you provide. This means that if you provide an object on a CRS that is not EPSG:3857 (e.g. lonlat EPSG:4326 on your example) maptiles would reproject the raster. As reprojecting rasters is quite different than reprojection vectors (sf) the reprojected raster would look deformed.

    See a reprex and also this related question: setting CRS for plotting with stamenmaps

    library(USAboundaries)
library(ggplot2)
library(sf)


us.state.longlat <- us_states()
us.state.longlat.cropped.large <- sf::st_crop(
  us.state.longlat,
  c(xmin = -100, xmax = -80, ymin = 30, ymax = 50)
)

# Since tiles are usually in EPSG:3857 this would avoid deformations
us.state.mercator.cropped.large <- st_transform(us.state.longlat.cropped.large, 3857)

library(maptiles)
library(tidyterra)

tile <- get_tiles(us.state.mercator.cropped.large, "Stamen.Terrain",
  zoom = 5,
  crop = TRUE
)

# Crop exactly to extent
tile_crop <- terra::crop(tile, us.state.mercator.cropped.large)

ggplot() +
  geom_spatraster_rgb(data = tile_crop) +
  geom_sf(
    data = us.state.mercator.cropped.large,
    color = "red",
    fill = "transparent"
  ) +
  coord_sf(expand = FALSE)

    enter image description here