st_or: 'sf' equivalent of QGIS Union

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

sf equivalent of QGIS Union

Description

sf equivalent of QGIS Union

Usage

st_or(x, y, dim = 2, suffix = c(".x", ".y"), suffix.all = FALSE, ...)

Arguments

x	object of class sf, sfc or sfg
y	object of class sf, sfc or sfg
dim	integer: A combination of 0, 1, and/or 2 (default) that constrains the dimension(s) of the returned geometries. 0 for points, 1 for lines, 2 for surfaces.
suffix	length 2 character vector (default c(".x", ".y")) attached as suffix to attribute headings inherited from argument x resp. from y
suffix.all	TRUE suffixes all attribute headings inherited from arguments x and y according to suffix; FALSE (default) suffixes only homonymous attribute headings.
...	arguments passed on to s2_options

Details

st_or() consists at its core of code presented by TimSalbim on gis.stackexchange . In addition to its precursor this version of st_or() includes:

• a more robust version of the internal function st_erase() equivalent to st_erase_robust

• the ability to handle homonymous attribute variables of both input geometry sets (s. below examples)

• the possibility to give customized suffixes to attribute variables corresponding to the geometry set they originated from (s. below examples)

• handling of input layers with differently named geometry columns and/or being totally overlapped by the other input layer.

• restricts geometric operations to geometries of the input layers x and y overlapping those of the other input layer (same as st_erase_robust).

Value

geometry set containing the intersection of x and y and the non-overlapping parts of x and y. The attribute table is filled with attribute values from the respective original geometry set for non-overlapping features, and attribute values from both geometry sets for overlapping features. The dimension(s) of the returned geometries is/are determined with the argument dim.

Examples

library(dplyr)
library(sf)

st_or(poly_1, poly_2) %>% plot()

st_or(poly_1, poly_2) %>% st_drop_geometry()

# to avoid repeated warning messages triggered by non-spatially constant ...
# ... attribute variables, set them all to "constant":
st_agr(poly_1) <- "constant"
st_agr(poly_2) <- "constant"

# Give customized suffixes to homonymous attributes of layers x and y:
st_or(poly_1, poly_2, suffix = c("_poly_1", "_poly_2")) %>% plot()

# If only homonymous attributes from one layer should get a suffix, set ...
# ... the suffix for the other layer to an empty string:
st_or(poly_1, poly_2, suffix = c("", ".y")) %>% names()

# If all attributes attributes of both layers x and y should get a
# ... layer-specific suffix, set suffix.all = TRUE:
st_or(poly_1, poly_2, suffix.all = TRUE) %>% names()

# If only all attributes from one layer should get a suffix, set the suffix ...
# ... for the other layer to an empty string and set suffix.all = TRUE:
st_or(poly_1, poly_2, suffix = c("", ".y"), suffix.all = TRUE) %>% names()

# create two layers with overlapping linestrings:
ls1 <- st_linestring(cbind(c(0, 1, 1, 0), c(0:3)))
ls2 <- st_linestring(cbind(c(2, 1, 1), c(0, 0, 3)))
ls3 <- st_linestring(cbind(c(0, 0.5, 0.5, 0), c(0, 0, 2.5, 2)))
A <- st_sf(id_A = 1, A = "A", geom = st_sfc(ls1), agr = "constant")
B <- st_sf(id_B = 1:2, B = "B", geom = st_sfc(ls2, ls3), agr = "constant")

plot(st_geometry(A), col = "gray", lwd = 4, extent = st_bbox_common(A, B))
plot(st_geometry(B), col = "red", lty = 2, add = TRUE)
legend("right", legend = c("A", "B"), col = c("gray", "red"), lwd = c(4, 1), lty = c(1, 2))

# when both input layers consist of linestings, and if the default specification ...
# ... dim = 2 (for surfaces / (multi)polygons) is used, a sf-object with zero ...
# ... rows will be returned:
st_or(A, B)

# to get lines returned set dim = 1:
comb_dim_1 <- st_or(A, B, dim = 1) %>%
  mutate(comb = ifelse(is.na(A), "B", ifelse(is.na(B), "A", "A+B")))
plot(comb_dim_1[, "comb"], lwd = 3, key.pos = 1, main = "dim = 1: only lines")

# for returning lines and points where lines cross or touch lines of the other ...
# ... input layer, set dim = c(0, 1):
comb_dim_0_1 <- st_or(A, B, dim = c(0, 1)) %>%
  mutate(comb = ifelse(is.na(A), "B", ifelse(is.na(B), "A", "A+B")))
plot(comb_dim_0_1[, "comb"], lwd = 3, cex = 2, key.pos = 1, main = "dim = c(0, 1): points & lines")

all.equal(
  st_or(A, B, dim = c(0, 1)), # returns points & lines
  st_or(A, B, dim = c(0, 1, 2)) # returns points, lines (& if available surfaces)
)

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