I am trying to plot a simple polygon on a map to denote the area I am interested in. To date, I have defined the polygon as and am able to plot it on it's own.
poly <- st_polygon(list(as.matrix(data.frame(lat = c(40, 40, 60, 60, 40),
lon = c(-60, -40, -40, -60, -60))))) #EDIT: these are lat/lon coordinates
ggplot() +
geom_sf(data = poly,
aes(geometry = geometry),
col = "red)
However, I need to add a basemap so that I can show the placement of this polygon in relation to other spatial elements.
# Attempt #1: stamenmaps basemap
grid_box <- c(left = -70,
right = -30,
bottom = 35,
top = 70)
base <- get_stamenmap(grid_box, zoom = 6, maptype = "terrain-background", color = "color")
ggmap(base) +
geom_sf(data = poly,
aes(geometry = geometry),
col = "red")
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Error in `geom_sf()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 4th layer.
Caused by error in `FUN()`:
! object 'lon' not found
Run `rlang::last_error()` to see where the error occurred.
The only approach I found for adding a crs to my polygon is below (st_transform & st_as_sf didn't work), but this dramatically changed the scale of the coordinates and still doesn't plot.
# Attempt #2: new CRS
poly <- st_polygon(list(as.matrix(data.frame(lon = c(-62, -43, -43, -62, -62),
lat = c(43, 43, 70, 70, 43))))) %>%
st_sfc(crs = 3857)
ggmap(base) +
geom_sf(data = poly,
aes(geometry = geometry),
col = "red")
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Error in `geom_sf()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 4th layer.
Caused by error in `FUN()`:
! object 'lon' not found
Run `rlang::last_error()` to see where the error occurred.
How can I get this polygon to plot over my basemaps?
I believe you are on the right track with sf::st_sfc()
call in your second try. To make it work you have to think a bit what your coordinates mean?
In EPSG:3857 the coordinates are defined in meters. This seems not to be the case with your original columns named lat & long, which usually stand for latitude and longitude in degrees.
Should your coordinates be in degrees (which seems likely, though I could not find it explicitly stated in your question) you will be better off with EPSG:4326.
A trick you have to work around is the ggmap being labeled in EPSG:4326 (unprojected degrees) but actually drawn in EPSG:3857 (projected meters in Mercator) a possible approach is described over here How to put a geom_sf produced map on top of a ggmap produced raster
Building upon the question linked I propose the following:
library(sf)
library(ggmap)
poly <- st_polygon(list(as.matrix(data.frame(lon = c(-62, -43, -43, -62, -62),
lat = c(43, 43, 70, 70, 43))))) %>%
st_sfc(crs = 4326)
grid_box <- c(left = -70,
right = -30,
bottom = 35,
top = 70)
base <- get_stamenmap(grid_box, zoom = 5,
maptype = "terrain-background",
color = "color")
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:
base <- ggmap_bbox(base)
ggmap(base) +
coord_sf(crs = st_crs(3857)) +
geom_sf(data = st_transform(poly, 3857),
col = "red", fill = NA,
inherit.aes = F) +
theme(axis.title = element_blank())