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:

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:

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)