Mention and demonstrate `st_relate()` in Chapter 4

[Robinlovelace]

Currently there is no mention of the important function st_relate() in Chapter 4 or anywhere in the book.

[defuneste]

Ok, so it was a bit too hard to work directly on chapter 4, so I decided to write a bit more before deciding what should be kept.

The intersection matrix is the foundation of most geometric relationships supported by the OGC SQL/MM standard [...]. (PostGIS in Action 3ed p269)

The intersection matrix is 3 x 3 matrix that defines all possible pair-wise combinations of exterior, boundary and interior when two geometries interact. (PostGIS in Action 3ed p269)

Geometries can be polygons, lines and points. Polygons, two dimensional objects, are delimited by their boundaries (a line) and they can have an interior (an area) and an exterior (another area). Other geometries will mostly have lower dimensions in boundary and interior.

TODO : complete the schema

Each part of the matrix can get a value. The value is related to the dimensions of the returned object.

• 0 they intersect only on points (st_dimension return 0)
• 1 they intersect with lines (st_dimension return 1)
• 2 they intersect with an area (st_dimension return 2)
• T they interset and number of dimensions $>=$ 0 (ie all of the above)
• F No intersection
• * it doesn't matter

Here this the matrix of all of the possible options (in the schema TODO:

	Interior	Boundary	Exterior
Interior	2	1	2
Boundary	1	0	1
Exterior	2	1	2

This flattened matrix will look like this: "212101212". This is called a pattern.

This matrix can be used in two ways: describe a relation between polygons or specify a type or relation you are interested in. As a result: you can use st_relate in two different ways:

• the first is to compute the DE-9IM relations between two objects
• the second is to check matching patterns

Compute the DE-9IM relations between two objects

To have a better time reading the pattern let's make a quick function:

pattern_de_9im =  "2FFF1FFF2"

matrix_de_9im = function(pattern) {
    string = unlist(strsplit(pattern , ""))
    matrix_de_9im = matrix(string, nrow = 3, byrow = TRUE)
    colnames(matrix_de_9im) = c("I", "B", "E")
    row.names(matrix_de_9im) = c("I", "B", "E")
    return(matrix_de_9im)
}

matrix_de_9im(pattern_de_9im)

	I	B	E
I	"2"	"F"	"F"
B	"F"	"1"	"F"
E	"F"	"F"	"2"

I will enhance an example from the sf documentation.

library(sf)
# 4 polygones to start 
pol1 = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)))) - 0.5
pol2 = pol1 + 1
pol3 = pol1 + 1.5
pol4 = st_polygon(list(rbind(c(-1,0.5),c(-1,1.5),c(0,1.5),c(0,0.5),c(-1,0.5))))

squares = sf::st_sf(data.frame(name = LETTERS[1:4],
           sf::st_sfc(pol1, pol2, pol3, pol4))
          )

plot(squares)

Here we can see how each square relates to the others:

my_first_st_relate = sf::st_relate(squares, squares)
colnames(my_first_st_relate) = squares$name
row.names(my_first_st_relate) = squares$name
my_first_st_relate

	A	B	C	D
A	"2FFF1FFF2"	"FF2F01212"	"FF2FF1212"	"FF2F11212"
B	"FF2F01212"	"2FFF1FFF2"	"212101212"	"FF2FF1212"
C	"FF2FF1212"	"212101212"	"2FFF1FFF2"	"FF2FF1212"
D	"FF2F11212"	"FF2FF1212"	"FF2FF1212"	"2FFF1FFF2"

We are producing a symmetrical matrix and you just need to focus on half of it. The first column describes how the A square is related to the others squares.

First line, first column [A,A] describes how the same square relates to itself:

matrix_de_9im(my_first_st_relate["A","A"])

	I	B	E
I	"2"	"F"	"F"
B	"F"	"1"	"F"
E	"F"	"F"	"2"

They share the same interior, same boundary and same exterior.

Let's see a more complicated case: square B with C:

matrix_de_9im(my_first_st_relate["B","C"])

	I	B	E
I	"2"	"1"	"2"
B	"1"	"0"	"1"
E	"2"	"1"	"2"

They share a part of their interior, their interior covers a part of their boundaries, a part of each interior is also an exterior of the other and their boundaries are related on two points.

Feel free to experiment with the other squares!

Check matching patterns

Now that we are familiar with the pattern, we can search the matrix for particular pattern. The fact that patterns are just strings allows us to also use regular expressions or other string shenanigans. Here, for example, we can return every pair-wise relation that shares a boundary.

# . represent every single character 
matrix(grepl(pattern = "....1....", my_first_st_relate),
       nrow = 4, byrow = TRUE) 

	[,1]	[,2]	[,3]	[,4]
[1,]	TRUE	FALSE	FALSE	TRUE
[2,]	FALSE	TRUE	FALSE	FALSE
[3,]	TRUE	FALSE	TRUE	FALSE
[4,]	TRUE	FALSE	FALSE	TRUE

We can also translate spatial predicates with DE-9IM pattern: see the wikipedia page.

We can achieve the same result thanks to the pattern argument of the st_relate function. The syntax of the pattern changes a bit, you need to replace . (that means every single character) to * to match DE-9IM rules.

sf::st_relate(squares, squares,  pattern = "****1****",  sparse = FALSE)

	[,1]	[,2]	[,3]	[,4]
[1,]	TRUE	FALSE	FALSE	TRUE
[2,]	FALSE	TRUE	FALSE	FALSE
[3,]	TRUE	FALSE	TRUE	FALSE
[4,]	TRUE	FALSE	FALSE	TRUE

With sparse = FALSE it will return the matrix. If you change it to TRUE you get a a sparse geometry binary predicate in the form of a list (cf. chapter 4)

sf::st_relate(squares, squares,  pattern = "****1****",  sparse = TRUE)

Sparse geometry binary predicate list of length 4, where the predicate was `relate_pattern'
1: 1, 4
2: 2
3: 3
4: 1, 4

Sources cited:

Sf's doc: https://r-spatial.github.io/sf/reference/st_relate.html

Wiki: https://en.wikipedia.org/wiki/DE-9IM

Postgis's doc: https://postgis.net/docs/using_postgis_query.html#eval_spatial_rel

PostGIS in action 3ed: p268-273

Geotools: https://docs.geotools.org/latest/userguide/library/jts/dim9.html

Egenhofer M, Mark (1995) Categorizing binary topoligical relationships between regions, lines, and points in geographic databases. Technical Report, Department of Surveying Engineering, University of Maine, Orono, ME.