# Higher-order geometry differences and intersections

Suppose you have the following geometry, consisting of three overlapping square polygons:

```
library(sf)
## Loading required package: methods
## Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3
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):

```
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
you can get an idea how ugly this can get.

## st_intersection

Now, inspired by a meticulously prepared pull request by Jeffrey Hanson, it suffices to do

```
(i = st_intersection(b))
## Geometry set for 7 features
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 0 ymin: 0 xmax: 1.8 ymax: 1.8
## epsg (SRID): NA
## proj4string: NA
## First 5 geometries:
## POLYGON ((1 0.2, 1 0, 0 0, 0 1, 0.2 1, 0.2 0.8,...
## POLYGON ((1 0.8, 1 0.2, 0.8 0.2, 0.8 0.8, 1 0.8))
## POLYGON ((1.2 1.2, 1.8 1.2, 1.8 0.2, 1 0.2, 1 0...
## POLYGON ((0.2 1, 0.8 1, 0.8 0.8, 0.2 0.8, 0.2 1))
## POLYGON ((0.8 1, 1 1, 1 0.8, 0.8 0.8, 0.8 1))
```

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

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

```
st_overlaps(i)
## Sparse geometry binary predicate list of length 7, where the predicate was `overlaps'
## 1: (empty)
## 2: (empty)
## 3: (empty)
## 4: (empty)
## 5: (empty)
## 6: (empty)
## 7: (empty)
st_overlaps(d)
## Sparse geometry binary predicate list of length 3, where the predicate was `overlaps'
## 1: (empty)
## 2: (empty)
## 3: (empty)
```

## Further reading

Jeffrey’s pull request is worth reading; the sf pkgdown site contains some further examples with squares.