I would like to combine a shapefile map of some US states on top of a raster map obtained from using library
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, ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER, 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.
ggmapalso. No luck. I looked at the documentation in a couple of places to find if there is a default projection for the
ggmapoptions and could find none. Google tends to use 3857 with a lot of their products. Transformed the states vector to that and it didn't help. This says to use 4326 with vector data. I tried use a different U.S. States polygon set I have on hand and got the same problem.