CODEWITHSUNDEEP
rggplot2shapefiler-sf
SkiFlyer
SkiFlyer

Reputation: 33

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

Upvotes: 3

Views: 167

Answers (1)

margusl
margusl

Reputation: 17554

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

Upvotes: 6

Related Questions