setting CRS for plotting with stamenmaps

Question

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?

Answer 1

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())