In this chapter, we will briefly describe a few geometric operations made available by the sf package, on vector data.
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 select the roads that intersect the municipality of Gramat. We do not modify their geometry.
library(sf)
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) 
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...mun_c <- st_centroid(mun)#> Warning: st_centroid assumes attributes are constant over geometriesmf_map(mun)
mf_map(mun_c, add = TRUE, cex = 1.2, col = "red", pch = 20)
dep_46 <- st_union(mun)
mf_map(mun, col = "lightblue")
mf_map(dep_46, col = NA, border = "red", lwd = 2, add = TRUE)
In this example, we aggregate polygons according to their hierarchical status (simple municipality, subprefecture, prefecture). We can aggregate the geometries (union) as well as their attributes (sum of population).
# 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)
)
mun_u#> Simple feature collection with 3 features and 2 fields
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93
#> STATUT POPULATION geometry
#> 1 Commune simple 140259 POLYGON ((554732.7 6356281,...
#> 2 Préfecture 19907 MULTIPOLYGON (((580994.5 63...
#> 3 Sous-préfecture 13763 MULTIPOLYGON (((626296.4 63...mf_map(mun_u, col = c("wheat", "aquamarine4", "yellow3"))
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)
Using st_intersection(), we can cut one layer by another. Here, geometries are modified.
# 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 geometriesmf_map(mun)
mf_map(mun_z, col = "red", border = "green", add = TRUE)
mf_map(mun_z)
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)
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] TRUETo 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.
The rmapshaper package (Teucher and Russell 2025) 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)
