--- layout: post title: "Higher-order geometry differences and intersections" date: "`r format(Sys.time(), '%d %B, %Y')`" comments: true author: Edzer Pebesma categories: r --- TOC [DOWNLOADHERE] Suppose you have the following geometry, consisting of three overlapping square polygons: ```{r echo=TRUE, fig=TRUE, fig.path = "images/", label="geoms1"} library(sf) pol = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))) b = st_sfc(pol, pol + c(.8, .2), pol + c(.2, .8)) par(mar = rep(0, 4)) plot(b, col = NA) ``` and you are interested in the area where all squares overlap (green), or where exactly two squares overlap (orange): ```{r echo=TRUE, fig=TRUE, fig.path = "images/", label="geoms2"} i = st_intersection(st_sf(b)) par(mar = rep(0, 4)) cl = sf.colors(3, categorical = TRUE) plot(b) plot(i[i$n.overlaps == 3,2], col = cl[1], add = TRUE) plot(i[i$n.overlaps == 2,2], col = cl[2], add = TRUE) ``` So far, with package `sf` or `rgeos` you could only get pairwise intersections, meaning you would have to go through all pairwise intersections and see whether they are intersected by others geometries or intersections. In [this StackOverflow question](https://stackoverflow.com/questions/44631044/efficient-extraction-of-all-sub-polygons-generated-by-self-intersecting-features) you can get an idea how ugly this can get. ## st\_intersection Now, inspired by a meticulously prepared [pull request](https://github.com/r-spatial/sf/pull/598) by [Jeffrey Hanson](http://jeffrey-hanson.com/), it suffices to do ```{r} (i = st_intersection(b)) ``` to get all the unique pieces, for each unique piece the number of contributing geometries, and a list-column with indexes of the geometries that contribute (overlap) for a particular piece. ## st\_difference The pull request Jeffrey wrote was to remove (erase) all overlapping pieces, which you now get by ```{r echo=TRUE, fig=TRUE, fig.path = "images/", label="geoms3"} d = st_difference(b) plot(d, col = cl) ``` For this latter approach, obviously the input order matters: what is returned are non-empty geometries with $x_1$, $x_2 - x_1$, $x_3 - x_2 - x_1$ etc. To prove that these intersections or differences do not have any overlaps, we can compute overlaps by ```{r} st_overlaps(i) st_overlaps(d) ``` ## Further reading Jeffrey's [pull request](https://github.com/r-spatial/sf/pull/598) is worth reading; the sf [pkgdown site](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) contains some further examples with squares.