Issue projecting polygon that crosses dateline R

Question

I have a polygon in a North America equal area projection and I would like to transform to unprojected long/lat. I would like to use it to crop a global raster in unprojected longlat.

The problem is that it crosses the dateline, and returns a problematic polygon. Is there a way to instead return two polygons, one on either side of the dateline?

I'm also open to other suggestions.

Note: The raster I am aiming to crop is global at 1km resolution. So I am trying to avoid needing to project that raster, as that would be computationally heavy.

library(terra)
library(sf)
library(rnaturalearth)

ee <- ext(c(-4342031, 3457969, -3389091, 4710909))
ee <- as.polygons(ee, crs = "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

ee <- st_transform(st_as_sf(ee), 4326)

world <- ne_coastline()
plot(world, lwd =0.5)
plot(ee, add=T)

Answer 1

There is a function called sf::st_break_antimeridian() that is great for this type of thing, however, it only works for non-projected coordinates (GCS).

As your ee object is projected and becomes 'problematic' when converted to WGS84/EPSG:4326, a simple solution is to:

• create a modified version of the WGS84/EPSG:4326 CRS, offsetting the central meridian to ensure the extent of your ee object does not cross the anti-meridian (in this example. 180 works).
• use sf::st_break_antimeridian() to split/break ee at the anti-meridian
• project back to the standard WGS84/EPSG:4326 CRS

As your SpatRaster is high resolution, e.g. 1km cell size, you may need to crop() and mask() in two steps. When I tried in a single step, the values in the result were not correct.

library(terra)
library(sf)
library(rnaturalearth)
library(ggplot2)
library(tidyterra)

# CRS projection strings
prj_aea <- "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
prj_wgs84_180 <- "+proj=longlat +datum=WGS84 +lon_0=180"

# Create sf object, project to modified WGS84/EPSG:4326 CRS, split at anti-meridian,
# project back to standard WGS84/EPSG:4326
ee <- as.polygons(ext(c(-4342031, 3457969, -3389091, 4710909)), crs = prj_aea) |>
  st_as_sf() |>
  st_transform(prj_wgs84_180) |>
  st_break_antimeridian(lon_0 = 180) |>
  st_transform(4326)

# Example SpatRaster with a global extent at approximately 1km resolution
r <- rast(ext(-180, 180, -90, 90), res = 1000/111320, crs = "EPSG:4326")
set.seed(1)
r[] <- sample(1:5, ncell(r), replace = TRUE)

# Crop, then mask r using ee
r1 <- crop(r, ee)
r1 <- mask(r1, ee)

# Plot
ggplot() +
  geom_spatraster(data = r1) +
  scale_fill_gradient2(low = "#0072B2",
                       mid = "#ffff6d", midpoint = 3,
                       high = "#C74A8F",
                       na.value = "transparent") +
  geom_sf(data = ne_coastline(), colour = "grey45", linewidth = 0.1) +
  coord_sf(expand = FALSE) +
  theme(panel.background = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

Below is the original answer, practical only if your raster data are relatively low resolution

library(terra)
library(sf)
library(rnaturalearth)
library(ggplot2)

# CRS projection strings
prj_aea <- "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
prj_wgs84_180 <- "+proj=longlat +datum=WGS84 +lon_0=180"
 
# Create sf object, project to modified WGS84 / EPSG:4326 CRS
ee <- as.polygons(ext(c(-4342031, 3457969, -3389091, 4710909)), crs = prj_aea) |>
  st_as_sf() |>
  st_transform(prj_wgs84_180)

# Create world sf, split at anti-meridian, project to modified WGS84 / EPSG:4326
# Will throw a warning, which you can safely ignore
world <- ne_coastline() |>
  st_break_antimeridian(lon_0 = 180) |>
  st_transform(prj_wgs84_180)

ggplot() +
  geom_sf(data = world, colour = "grey") +
  geom_sf(data = ee, fill = NA, linewidth = 1) +
  theme_void()

Here is the workflow for cropping your raster. This assumes your raster is a SpatRaster with WGS84/EPSG:4326 as its CRS. Ensure the nrow and ncol values for the r1 object match the dimensions of your actual data.

# Create example SpatRaster with global extent to represent your data
r <- rast(ext(-180, 180, -90, 90), nrow = 180, ncol = 360, crs = "EPSG:4326")
set.seed(1)
r[] <- sample(1:5, ncell(r), replace = TRUE)

# Project r to same CRS as ee
r1 <- rast(ext(0, 360, -90, 90), nrow = 180, ncol = 360, crs = prj_wgs84_180)
r1 <- project(r, r1)
r1 <- rotate(r1, 180)

# Crop r1 using ee, project back to standard WGS84/EPSG:4326
r2 <- crop(r1, ee, mask = TRUE)
r3 <- project(r2, r)

Plot the results to compare:

library(tidyterra)
library(patchwork)

p1 <- ggplot() +
  geom_spatraster(data = r2) +
  scale_fill_gradient2(low = "#0072B2", 
                       mid = "#ffff6d", midpoint = 3, 
                       high = "#C74A8F",
                      na.value = "transparent") +
  geom_sf(data = world, colour = "grey45", linewidth = 0.1) +
  labs(title = "Projected to custom CRS and cropped") +
  coord_sf(expand = FALSE) +
  theme(panel.background = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

p2 <- ggplot() +
  geom_spatraster(data = r3) +
  scale_fill_gradient2(low = "#0072B2", 
                       mid = "#ffff6d", midpoint = 3, 
                       high = "#C74A8F",
                       na.value = "transparent") +
  geom_sf(data = ne_coastline(), colour = "grey45", linewidth = 0.1) +
  labs(title = "Cropped data projected back to standard WGS84") +
  coord_sf(expand = FALSE) +
  theme(legend.position = "none",
        panel.background = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

p1 + p2

Note that the above image is relatively low resolution, so if you zoom in, there are discrepancies in some of the like-for-like cell colours (particularly around the edges). This issue does not persist at 600dpi.