st_erase_robust: Erasing from geometry set 'x' all parts overlapped by...

In a-benini/sfhelpers: Helper Functions for Wrapping Simple Features-Related Routines

Erasing from geometry set x all parts overlapped by geometry set y

Description

Erasing from geometry set x all parts overlapped by geometry set y

Usage

st_erase_robust(x, y, ...)

Arguments

x	object of class sf, sfc or sfg
y	object of class sf, sfc or sfg
...	arguments passed on to s2_options

Details

The example section of the sf-package help page on geometric operations on pairs of simple feature geometry sets (geos_binary_ops) presents code for a helper function that erases all parts from geometry set x overlapped by geometry set y. This function sometimes works as expected, sometimes it doesn't. (s. examples below).

Even when both input layers x and y consist of valid geometries (which can be checked with st_is_valid, respectively fixed with st_make_valid), the a. m. helper function can still throw an error. This is often due to internally applying st_union(st_combine()) to y. st_combine is sometimes also the cause for returning incorrectly erased geometry. To avoid these issues, st_erase_robust() uses st_union only. Leaving out st_combine() may add to the complexity of the involved geometries, thus potentially increasing the run time. To counteract this, st_erase_robust() detects with the help of st_intersects which geometries from the input x and y overlap with those of the other layer and applies geometric operations only to such geometries.

When using the helper function from the a.m. sf help page with input having longlat degrees CRS, switching off spherical geometry (s2) by setting sf_use_s2 to FALSE can help to overcome issues. But this isn't a safe workaround for all issues caused by st_combine() when it comes to erasing.

Value

Returns all parts of geometry set x not overlapped by geometry set y

Examples

library(sf)

# find code of helper function st_erase():
# ?geos_binary_ops

# copy function code:
st_erase <- function(x, y) st_difference(x, st_union(st_combine(y)))

# get some demo data:
nc   <- st_read(system.file("gpkg/nc.gpkg", package = "sf"), quiet = TRUE)
ext  <- st_bbox(nc) + rep(c(-0.1, 0.1), each = 2)
grid <- st_make_grid(ext) %>% st_sf(id = seq_along(.), geom = ., agr = "constant")

st_is_longlat(nc) # demo data has a longlat degrees crs

sf_use_s2(TRUE)
# check if helper function works with demo data:
try(st_erase(grid, nc))

# internal processing of input y (nc) returns the same error as st_erase():
try(st_union(st_combine(nc)))

# st_erase_robust() can handle this:
st_erase_robust(grid, nc) %>% plot()

a-benini/sfhelpers documentation built on Aug. 28, 2024, 3:30 a.m.