rggplot2shapefiler-sf

R and ggplot with st_crop: map cropping does not work as wanted


I am trying to plot only a part of the world map, limited by a square with limits lon (-30, 90) and lat (30, 82). When I try to crop the map with sf_crop, it does not return the desired square when plotting with ggplot, even though i make sure the used world map has the "normal" crs with specifying crs = 4326. I know I could probably introduce map limits with coord_sf, but I need an alternative solution. I would really appreciate any help.

I tried the following code that produces the described problem:

lapply(c("rnaturalearthdata", "sf", "ggplot2","dplyr"), require, character.only = TRUE)

world <- ne_coastline(scale = "small", returnclass = "sf") %>%
  st_transform(world, crs = 4326)
world_cropped <- st_crop(world, xmin = -30, xmax = 90, ymin = 30, ymax = 82)

ggplot() +
  geom_sf(data = world, fill = NA, color = "black")

ggplot() +
  geom_sf(data = world_cropped, fill = NA, color = "black")

First map, as expected

Second cropped map, did not cut a "square" as I expected, strange result

Highlighted (by hand) the map square I want to be shown approximatelly


Solution

  • When working with geographic coordinates (WGS84 / EPSG:4326), sf by default uses spherical geometry operations through s2. So your bounding box defined by 4 corners is not really a rectangle; as sides follow shortest paths between 2 points on a sphere, 2 sides running east-west become arcs when projected on a plane. We can visualize it by adding more points to "straight" edges of bbox with st_segmentize():

    library(rnaturalearth)
    library(ggplot2)
    library(sf)
    
    world <- ne_coastline(scale = "small") 
    bb_wgs84 <-  
      st_bbox(c(xmin = -30, xmax = 90, ymin = 30, ymax = 82), crs = "WGS84") |> 
      st_as_sfc()
    
    world_cropped <- st_crop(world, bb_wgs84)
    ggplot() +
      geom_sf(data = st_segmentize(bb_wgs84, 1000), fill = "lightgreen", alpha = .2) +
      geom_sf(data = bb_wgs84, fill = NA, color = "darkgreen", linewidth = 1) +
      geom_sf(data = world_cropped, fill = NA)
    

    Probably easiest workaround is switching s2 off for the crop operation so coordinates are treated as planar coordinates, which allows that rectangular cropping you are probably after.

    sf_use_s2(use_s2 = FALSE)
    #> Spherical geometry (s2) switched off
    
    world_cropped <- st_crop(world, xmin = -30, xmax = 90, ymin = 30, ymax = 82)
    #> although coordinates are longitude/latitude, st_intersection assumes that they
    #> are planar
    #> Warning: attribute variables are assumed to be spatially constant throughout
    #> all geometries
    sf_use_s2(use_s2 = TRUE)
    #> Spherical geometry (s2) switched on
    
    ggplot() +
      geom_sf(data = st_segmentize(bb_wgs84, 1000), fill = "lightgreen", alpha = .2) +
      geom_sf(data = bb_wgs84, fill = NA, color = "darkgreen", linewidth = 1) +
      geom_sf(data = world_cropped, fill = NA)
    

    Created on 2024-07-09 with reprex v2.1.0