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