rr-sfggmap

setting CRS for plotting with stamenmaps


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?

Solution

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

    map of area of interest - red box on stamen map