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

Answer 2

Although already answered, I would like to share another approach for this issue. You can extract easily tiles with maptiles and plot that seamlessly with tidyterra. The only thing to bear in mind is that tiles are provided usually in EPSG:3857, so you may want to get and plot the tiles on that projection in order to be displayed with full detail. Otherwise the projection of a raster (that involves resampling) would create some deformations on the final images.

See:

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(maptiles)
library(ggplot2)
library(tidyterra)
#> 
#> Attaching package: 'tidyterra'
#> The following object is masked from 'package:stats':
#> 
#>     filter


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 to polygon 

grid_box <- c(-70, 35, -30, 70)
class(grid_box) <- "bbox"
grid_poly <- st_as_sfc(grid_box) %>% st_set_crs(4326) %>%
  st_transform(3857)
 

basemap <- get_tiles(grid_poly, provider = "Stamen.Terrain", crop = TRUE)


ggplot() +
  geom_spatraster_rgb(data = basemap) +
  geom_sf(data = poly, fill = NA, color = "red") +
  coord_sf(crs=3857, expand = FALSE)

^{Created on 2022-12-10 with reprex v2.0.2}