<- st_read("data/lot.gpkg", layer = "communes", quiet = TRUE)
mun <- st_read("data/lot.gpkg", layer = "routes", quiet = TRUE)
roads <- mun[mun$NOM_COM == "Gramat", ]
gramat
<- st_filter(x = roads,
road_gramat 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 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
3.1.2 Spatial join
st_join()
is used to perform spatial joins. Use the join
argument to select the geometric predicate.
<- st_join(x = roads,
road_gramat 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...
- Import the layers of municipalities and restaurants.
- 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
<- st_centroid(mun) mun_c
#> 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
<- st_union(mun)
dep_46
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
<- mun$STATUT
i
<- st_sf(
mun_u 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
<- mun[mun$NOM_COM == "Gramat", ]
gramat
<- st_buffer(x = gramat, dist = 5000)
gramat_b
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
<- st_geometry(gramat) |>
zone st_centroid() |>
st_buffer(10000)
mf_map(mun)
mf_map(zone, border = "red", col = NA, lwd = 2, add = TRUE)
<- st_intersection(x = mun, y = zone) mun_z
#> 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)
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
<- st_make_grid(x = mun, cellsize = 5000)
grid
# transform to sf, add identifier
<- st_sf(ID = 1:length(grid), geom = grid)
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()
.
<- st_filter(grid, dep_46, .predicate = st_intersects)
grid
# Import restaurants
<- st_read("data/lot.gpkg", layer = "restaurants", quiet = TRUE)
restaurant
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).
<- st_intersects(grid, restaurant, sparse = TRUE)
inter
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.
$nb_restaurant <- lengths(inter)
grid
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")
<- ms_simplify(mun, keep = 0.01 , keep_shapes = TRUE)
mun_simp1 <- ms_simplify(mun, keep = 0.001, keep_shapes = TRUE)
mun_simp2 mf_map(mun)
mf_map(mun_simp1)
mf_map(mun_simp2)
- Compute the number of restaurants per municipality.
- Which municipalities have more than 10 restaurants and fewer than 1000 inhabitants?
- Create a map showing all the municipalities in grey and the municipalities selected above in red.