3  Geoprocessing

3.1 Spatial selection and join

3.1.1 Spatial selection

st_filter() is used to perform spatial selections. The .predicate argument lets you pick the selection criterion by using one of the “geometric predicate” functions (e.g. st_intersects(), st_within(), st_crosses()…).
Here, we’ll select the roads that intersect the municipality of Gramat

mun <- st_read("data/lot.gpkg", layer = "communes", quiet = TRUE)
roads <- st_read("data/lot.gpkg", layer = "routes", quiet = TRUE)
gramat <-  mun[mun$NOM_COM == "Gramat", ]

road_gramat <-  st_filter(x = roads, 
                          y = gramat,
                          .predicate = st_intersects)
# Plot
library(mapsf)
mf_map(gramat, col = "lightblue")
mf_map(roads, add = TRUE)
mf_map(road_gramat, col = "tomato", lwd = 2, add = TRUE)  

3.1.2 Spatial join

st_join() is used to perform spatial joins. Use the join argument to select the geometric predicate.

road_gramat <-  st_join(x = roads,
                         y = mun[, "INSEE_COM"],
                         join = st_intersects,
                         left = FALSE)
road_gramat
#> Simple feature collection with 1247 features and 3 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 587147.6 ymin: 6394844 xmax: 608194.7 ymax: 6420006
#> Projected CRS: RGF93 v1 / Lambert-93
#> First 10 features:
#>     ID      CLASS_ADM INSEE_COM                           geom
#> 1    1 Départementale     46240 LINESTRING (590557.5 641181...
#> 2    2 Départementale     46240 LINESTRING (593733.2 641429...
#> 3    3 Départementale     46240 LINESTRING (590665 6412381,...
#> 4    4 Départementale     46128 LINESTRING (598940.9 640909...
#> 5    5 Départementale     46104 LINESTRING (603201.9 640181...
#> 6    6     Sans objet     46235 LINESTRING (598162.3 640108...
#> 7    7 Départementale     46090 LINESTRING (598887.3 639763...
#> 7.1  7 Départementale     46138 LINESTRING (598887.3 639763...
#> 7.2  7 Départementale     46233 LINESTRING (598887.3 639763...
#> 8    8 Départementale     46090 LINESTRING (601184.3 639697...
Exercice
  1. Import the layers of municipalities and restaurants.
  2. Perform a spatial join to find the name and identifier of the municipality where each restaurant is located.

3.2 Geometrical operations

3.2.1 Centroids

mun_c <- st_centroid(mun)
#> Warning: st_centroid assumes attributes are constant over geometries
mf_map(mun)
mf_map(mun_c, add = TRUE, cex = 1.2, col = "red", pch = 20)

3.2.2 Aggregate polygons

dep_46 <- st_union(mun)

mf_map(mun, col = "lightblue")
mf_map(dep_46, col = NA, border = "red", lwd = 2, add = TRUE)

3.2.3 Aggregate polygons using a grouping variable

# Grouping variable
i <- mun$STATUT 

mun_u <- st_sf(
  STATUT     = tapply(X = mun$STATUT     , INDEX = i, FUN = head, 1),
  POPULATION = tapply(X = mun$POPULATION , INDEX = i, FUN = sum), 
  geometry   = tapply(X = mun            , INDEX = i, FUN = st_union), 
  crs        = st_crs(mun)
) 

3.2.4 Buffers

st_buffer() is used to construct buffer zones. The distance is expressed in units of the projection (st_crs(x)$units).

# Select a municipality
gramat <- mun[mun$NOM_COM == "Gramat", ]

gramat_b <- st_buffer(x = gramat, dist = 5000)

mf_map(gramat_b, col = "lightblue", lwd=2, border = "red")
mf_map(gramat, add = TRUE, lwd = 2)

3.2.5 Intersection

Using st_intersection(), we can cut one layer by another.

# create a buffer zone around the centroid of Gramat
zone <- st_geometry(gramat) |> 
  st_centroid() |> 
  st_buffer(10000)

mf_map(mun)
mf_map(zone, border = "red", col = NA, lwd = 2, add = TRUE)

mun_z <- st_intersection(x = mun, y = zone)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
mf_map(mun)
mf_map(mun_z, col = "red", border = "green", add = TRUE)

mf_map(mun_z)

In this example, we’ve used pipes (|>). Pipes are used to concatenate a sequence of instructions.

3.2.6 Regular grid

st_make_grid() is used to create a regular grid. It produces an sfc object, and then you need to use st_sf() to transform this sfc object into an sf object. We also add a column of unique identifiers.

# create a grid
grid <- st_make_grid(x = mun, cellsize = 5000)

# transform to sf, add identifier
grid <- st_sf(ID = 1:length(grid), geom = grid)

mf_map(grid, col = "grey", border = "white")
mf_map(mun, col = NA, border = "grey50", add = TRUE)

3.2.7 Counting points in polygons

Select the grid tiles which intersect the department with st_filter().

grid <- st_filter(grid, dep_46, .predicate = st_intersects)

# Import restaurants
restaurant <- st_read("data/lot.gpkg", layer = "restaurants", quiet = TRUE)

mf_map(grid, col = "grey", border = "white")
mf_map(restaurant, pch = 20, col = "red", cex = .5, add = TRUE)

We then use st_intersects(..., sparse = TRUE), which will produce a list of the elements of the restaurant object inside each element of the grid object (via their indexes).

inter <- st_intersects(grid, restaurant, sparse = TRUE)

length(inter) == nrow(grid)
#> [1] TRUE

To count the number of restaurants, all we need to do is report the length of each of the elements in this list.

grid$nb_restaurant <- lengths(inter)

mf_map(grid)
mf_map(grid, var = "nb_restaurant", type = "prop")
#> 94 '0' values are not plotted on the map.

3.2.8 Geometries simplification

The rmapshaper package (Teucher and Russell 2023) builds on the Mapshaper JavaScript library (Bloch 2013) to offer several topology-friendly methods for simplifying geometries.
The keep argument is used to indicate the level of simplification. The keep_shapes argument is used to keep all polygons when the simplification level is high.

library("rmapshaper")
mun_simp1 <- ms_simplify(mun, keep = 0.01 , keep_shapes = TRUE)
mun_simp2 <- ms_simplify(mun, keep = 0.001, keep_shapes = TRUE)
mf_map(mun)
mf_map(mun_simp1)
mf_map(mun_simp2)

Exercice
  1. Compute the number of restaurants per municipality.
  2. Which municipalities have more than 10 restaurants and fewer than 1000 inhabitants?
  3. Create a map showing all the municipalities in grey and the municipalities selected above in red.