--- 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.