com_c <- st_centroid(com)#> Warning: st_centroid assumes attributes are constant over geometriesplot(st_geometry(com))
plot(st_geometry(com_c), add = TRUE, cex = 1.2, col = "red", pch = 20)
8 Opérations sur les géométries
8.1 Extraire des centroides
8.2 Agréger des polygones
dep_46 <- st_union(com)
plot(st_geometry(com), col = "lightblue")
plot(st_geometry(dep_46), add = TRUE, lwd = 2, border = "red")
8.2.1 Agréger des polygones en fonction d’une variable
com_u <- aggregate(x = com[,c("POPULATION")],
by = list(STATUT = com$STATUT),
FUN = "sum")
plot(com_u)
8.3 Construire une zone tampon
gramat <- com[com$NOM_COM == "Gramat", ]
gramat_b <- st_buffer(x = gramat, dist = 5000)
plot(st_geometry(gramat_b), col = "lightblue", lwd=2, border = "red")
plot(st_geometry(gramat), add = TRUE, lwd = 2)
8.4 Réaliser une intersection
En utilisant la fonction st_intersection() on va découper une couche par une autre.
# création d'une zone tampon autour du centroid de la commune de Gramat
# en utilisant le pipe
zone <- st_geometry(gramat) |>
st_centroid() |>
st_buffer(10000)
plot(st_geometry(com))
plot(zone, border = "red", lwd = 2, add = TRUE)
com_z <- st_intersection(x = com, y = zone)#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometriesplot(st_geometry(com))
plot(st_geometry(com_z), col="red", border="green", add=T)
plot(st_geometry(com_z))
8.5 Créer une grille régulière
La fonction st_make_grid() permet de créer une grille régulière. La fonction produit un objet sfc, il faut ensuite utiliser la fonction st_sf() pour transformer cet objet sfc en objet sf. Lors de cette transformation nous rajoutons ici une colonne d’identifiants uniques.
grid <- st_make_grid(x = com, cellsize = 2500)
grid <- st_sf(ID = 1:length(grid), geom = grid)
plot(st_geometry(grid), col = "grey", border = "white")
plot(st_geometry(com), border = "grey50", add = TRUE)
8.6 Compter des points dans un polygone
# selection des carreaux de la grille qui intersectent le département
inter <- st_intersects(grid, dep_46, sparse = FALSE)
grid <- grid[inter, ]
restaurant <- st_read("data/lot46.gpkg", layer = "restaurant", quiet = TRUE)
plot(st_geometry(grid), col = "grey", border = "white")
plot(st_geometry(restaurant), pch = 20, col = "red", add = TRUE, cex = .2)
inter <- st_intersects(grid, restaurant, sparse = TRUE)
length(inter)#> [1] 936Ici nous utilisons l’argument sparse = TRUE. L’objet inter est une liste de la longueur de l’objet grid et chaque élément de la liste contient l’index des éléments de l’objet restaurant qu’il intersecte.
Par exemple le carreau 795 intersecte les restaurants 2716 et 2718 :
inter[795]#> [[1]]
#> [1] 2716 2718plot(st_geometry(grid[795, ]))
plot(st_geometry(restaurant), add = T)
plot(st_geometry(restaurant[c(2716, 2718), ]),
col = "red", pch = 19, add = TRUE)
Pour compter le nombre de restaurants il suffit de reporter la longueur de chacun des éléments de cette liste.
grid$nb_restaurant <- lengths(inter)
plot(grid["nb_restaurant"])
8.7 Aggréger les valeurs de points dans des polygones
Dans cet exemple nous importons un fichier csv (pop.csv) qui contient les données d’une grille de population (x et y longitude et latitude, EPSG:3035; Ind, la population).
Une fois importé nous transformons le data.frame en objet sf.
# Import du ficier
pop_raw <- read.csv("data/pop.csv")
head(pop_raw)#> x y Ind
#> 1 3630500 2383500 23.0
#> 2 3629500 2384500 25.0
#> 3 3630500 2384500 18.0
#> 4 3631500 2384500 11.5
#> 5 3644500 2384500 5.0
#> 6 3645500 2384500 9.0# transformation en objet sf
pop <- st_as_sf(pop_raw, coords = c("x", "y"), crs = 3035)
# gestion de la projection
pop <- st_transform(pop, st_crs(com))
# Affichage
library(mapsf)
mf_map(com)
mf_map(pop, var = "Ind", type = "prop", add = TRUE,
inches = .2, col = "ivory", border = "ivory4")
L’objectif est d’agréger les valeur de ces points (les populations contenues dans le champ Ind) dans les communes du Lot.
En utilisant la fonction st_intersection() nous sélectionnons les points contenus dans une commune et rajoutons à chacun des points l’ensemble des attributs de la commune dans lequel il se trouve.
inter <- st_intersection(x = pop, y = com)#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometrieshead(inter)#> Simple feature collection with 6 features and 13 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 558290.8 ymin: 6372439 xmax: 561363.8 ymax: 6374337
#> Projected CRS: RGF93 / Lambert-93
#> Ind INSEE_COM NOM_COM STATUT POPULATION AGR_H AGR_F IND_H
#> 1069 3 46001 Albas Commune simple 522 4.978581 0 4.936153
#> 1070 7 46001 Albas Commune simple 522 4.978581 0 4.936153
#> 1129 11 46001 Albas Commune simple 522 4.978581 0 4.936153
#> 1130 9 46001 Albas Commune simple 522 4.978581 0 4.936153
#> 1131 13 46001 Albas Commune simple 522 4.978581 0 4.936153
#> 1189 18 46001 Albas Commune simple 522 4.978581 0 4.936153
#> IND_F BTP_H BTP_F TER_H TER_F geometry
#> 1069 0 9.957527 0 44.91714 34.6818 POINT (559462.5 6372439)
#> 1070 0 9.957527 0 44.91714 34.6818 POINT (560457.4 6372534)
#> 1129 0 9.957527 0 44.91714 34.6818 POINT (559374.1 6373435)
#> 1130 0 9.957527 0 44.91714 34.6818 POINT (560369 6373530)
#> 1131 0 9.957527 0 44.91714 34.6818 POINT (561363.8 6373625)
#> 1189 0 9.957527 0 44.91714 34.6818 POINT (558290.8 6374337)Nous pouvons ensuite utiliser la fonction aggregate() pour agréger la population par communes.
resultat <- aggregate(x = list(pop_from_grid = inter$Ind),
by = list(INSEE_COM = inter$INSEE_COM),
FUN = "sum")
head(resultat)#> INSEE_COM pop_from_grid
#> 1 46001 472.0
#> 2 46002 79.0
#> 3 46003 518.5
#> 4 46004 144.5
#> 5 46005 308.0
#> 6 46006 345.0On peut ensuite créer un nouvel objet avec ce résultat.
com_res <- merge(com, resultat, by = "INSEE_COM", all.x = TRUE)
head(com_res)#> Simple feature collection with 6 features and 13 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 556391.9 ymin: 6371852 xmax: 614866.5 ymax: 6418606
#> Projected CRS: RGF93 / Lambert-93
#> INSEE_COM NOM_COM STATUT POPULATION AGR_H AGR_F
#> 1 46001 Albas Commune simple 522 4.978581 0.000000
#> 2 46002 Albiac Commune simple 67 0.000000 9.589041
#> 3 46003 Alvignac Commune simple 706 10.419682 0.000000
#> 4 46004 Anglars Commune simple 219 0.000000 0.000000
#> 5 46005 Anglars-Juillac Commune simple 329 4.894895 4.894895
#> 6 46006 Anglars-Nozac Commune simple 377 4.840849 0.000000
#> IND_H IND_F BTP_H BTP_F TER_H TER_F pop_from_grid
#> 1 4.936153 0.000000 9.957527 0.000000 44.917145 34.681799 472.0
#> 2 0.000000 0.000000 4.794521 0.000000 4.794521 9.589041 79.0
#> 3 10.419682 5.209841 10.419682 0.000000 57.308249 78.147612 518.5
#> 4 20.000000 15.000000 10.000000 0.000000 20.000000 20.000000 144.5
#> 5 4.894895 0.000000 0.000000 0.000000 29.369369 29.369369 308.0
#> 6 0.000000 0.000000 9.681698 4.840849 43.567639 38.726790 345.0
#> geometry
#> 1 MULTIPOLYGON (((559262 6371...
#> 2 MULTIPOLYGON (((605540.7 64...
#> 3 MULTIPOLYGON (((593707.7 64...
#> 4 MULTIPOLYGON (((613211.3 64...
#> 5 MULTIPOLYGON (((556744.9 63...
#> 6 MULTIPOLYGON (((576667.2 64...plot(com_res["pop_from_grid"])
plot(com_res$POPULATION, com_res$pop_from_grid)
