Package cartography avec R
L'essentiel de cette page
Voici l'exposition de quelques fonctions bien utile de la librairie cartography.
Comment agréger des polygones, identifier les régions voisins, les fusionner...
0- Partons d'un exemple. Celui de la France
Le package raster permet charger très rapidement la France au format sp. Il est possible de convertir ensuite ce fond de carte au format sf afin de pouvoir le manipuler avec la librairie cartography.
library(raster)
fr <- getData('GADM', country='FRA', level=1)
plot(fr)
is(fr)
frsf <- st_as_sf(fr, "sf") # Conversion au format sf
st_bbox(frsf) # Identifier les limites
xmin ymin xmax ymax
-5.143751 41.333752 9.560416 51.089397
Affichage de la carte au format sf :
names(frsf)
plot(frsf) # Toutes les cartes possibles
plot(st_geometry(frsf)) # Seulemet les contours géométriques/géographiques
plot(frsf[4]) # Afficher les régions
plot(frsf["NAME_1"],key.pos=NULL) # Autre façon d'afficher les régions en inactivant la légende automatique
Convertir un objet sf en sp
Sp <- as(Sf, Class = "Spatial")
Convertir un objet sp en sf
Sf <- st_as_sf(Sp, "sf")
Cliquer sur ce lien pour aller dans la page sur les changements de projections.
st_crs() # Pour conntraître le système de projection
st_transform(avignon, crs = st_crs(avignon)) -> avignon
st_transform(avignon, "+proj=longlat +ellps=WGS84 +datum=WGS84") -> avignon
Convertir un jeu de données sf pour le passer en sp et le manipuler avec d'autres librairies comme raster
Ouverture d'un format GPKG
library(sf)
avignon <- st_read(dsn ="https://trouver.datasud.fr/dataset/99a53f28-5251-456e-a937-253a275a5808/resource/c3b6cf68-742e-4179-8092-02f69f617150/download/84007-monumentshistoriques.gpkg", quiet = TRUE)
# Geopackage GPKG Avignon - https://trouver.datasud.fr/organization/7033295c-7295-461b-9943-26bfb774a23b?res_format=GPKG
# Conversion en objet sp
sp <- as(avignon, Class = "Spatial")
is(sp)
plot(sp,xlim=c(844152,846250),ylim=c(6317526,6319009))
names(sp)
1- Travailler sur les polygones, les fusionner, les catégoriser
Réunir toutes les régions pour tracer le contours de la France (ici, en bleu)
mon_union <- st_union(frsf) # Fusionner des polygones
plot(st_geometry(frsf),border="red",lwd=3,legend=F)
plot(st_geometry(mon_union),border="blue",lwd=3,add=T) # Contours de la France
cf. aussi la fonction getBorders() de cartography.
Identifier toutes les régions limitrophes à une autre (intersection,s régions voisines)
aquitaine <- frsf[frsf$NAME_1 == "Nouvelle-Aquitaine", ]
intersections <- st_intersects(x = frsf, y = aquitaine, sparse = F)
plot(st_geometry(frsf[intersections,]), col = "#F7DC6F", add = TRUE)
Agréger des régions en fonction d'une variable
occitan <- c(rep(0,9),rep(1,2),rep(0,2))
frsf_u <- aggregate(x = frsf["geometry"],
by = list(occitan),
FUN = "sum")
plot(frsf_u,key.pos=NULL,main="Régions occitanes")
Découper la carte sur une fenêtre donné : intersection entre 2 formes un carré et le fond de carte
Contenu de contours : 5 points qui définissent un polygone rectangulaire.
x y
1 -1.788472 42.77767
2 2.609396 42.77767
3 2.609396 47.17576
4 -1.788472 47.17576
5 -1.788472 42.77767
cadre_aq <- st_bbox(aquitaine)
contours <- data.frame(
x=c(cadre_aq[1],cadre_aq[3],cadre_aq[3],
cadre_aq[1],cadre_aq[1]),
y=c(cadre_aq[2],cadre_aq[2],cadre_aq[4],
cadre_aq[4],cadre_aq[2])
)
cadre <- st_sf(st_sfc(st_polygon(list(as.matrix(contours)))), crs = st_crs(frsf))
intersection <- st_intersection(x = frsf, y = cadre) # ne pas confondre avec st_intersects
plot(st_geometry(intersection))
2- Centroïdes et distances inter-points
Extraire et afficher les centroïdes
frsf_c <- st_centroid(frsf)
plot(st_geometry(frsf))
plot(st_geometry(frsf_c), add=TRUE, cex=2, col="red", pch=20)
Établir une matrice de distances entre les centroïdes
library("lwgeom")
matrice_dist <- st_distance(x = frsf_c, y = frsf_c)
colnames(matrice_dist) <- frsf$NAME_1
rownames(matrice_dist) <- frsf$NAME_1
round(matrice_dist/1000,0) # en km
Cet extrait nous permet de voir les distances qui séparent quelques régions...
Units: [m]
Auvergne-Rhône-Alpes Bourgogne-Franche-Comté
Auvergne-Rhône-Alpes 0 192
Bourgogne-Franche-Comté 192 0
Bretagne 636 584
Centre-Val de Loire 310 238
Corse 524 659
Comptage des points par région (par polygone)
Etape 1
Commençons par générer des points aléatoires qui couvrent la carte de France
pts <- st_sample(x = frsf, size = 100)
Cette commande est très facile d'usage, mais dans les faits, j'aurais plus souvent un tableau de coordonnées à traiter. Voici donc ce qu'il faut faire si j'ai un tableau de coordonnées pour mes points :
# Je simule d'abord un tableau
st_bbox(frsf)->limites
x <- rnorm(100,(limites[1]+limites[3])/2,2)
y <- rnorm(100,(limites[2]+limites[4])/2,3)
# Je vais devoir ensuite le convertir en objet spatial
# Ne pas utiliser la fonction st_multipoint (elle donnera des erreurs)
# Utiliser plutôt st_as_sf en utilisant le même CRS que le fond de carte (st_crs())
pts<- st_as_sf(data.frame(x,y), coords = c("x", "y"),crs=st_crs(frsf))
Visualiser :
plot(st_geometry(frsf))
plot(pts,col="orange",cex=2,pch=16,add=T)
Etape 2
inter <- st_intersects(frsf, pts)
# Résultat du comptage
frsf$inter <- sapply(X = inter, FUN = length)
plot(st_geometry(frsf))
plot(st_geometry(frsf[frsf$inter>12,]), col = "red", add=TRUE)plot(st_geometry(frsf[frsf$inter<=12&frsf$inter>6,]), col = "orange", add=TRUE)plot(st_geometry(frsf[frsf$inter<=6,]), col = "green", add=TRUE)plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)3- Des fonctions cartography pour simplifier l'affichage (library(cartography))
# Fond de carte de la France et conversion
# en objet sf
library(raster)
fr <- getData('GADM', country='FRA', level=1)
frsf <- st_as_sf(fr, "sf") # Conversion au format sf
variable_de_couleur <- round(rnorm(length(frsf$NAME_1),6,3),0)
frsf2 <- cbind(frsf,variable_de_couleur)
# Palette de couleurs
colfunc<-colorRampPalette(c("white","red"))
colors <- (colfunc(5))
# Affichage
choroLayer(frsf2 , var="variable_de_couleur", col=colors)
# Fond de carte de la France et conversion
# en objet sf
library(raster)
fr <- getData('GADM', country='FRA', level=1)
frsf <- st_as_sf(fr, "sf") # Conversion au format sf
# Simuler les coordonnées de points
library("magrittr")
pts <- st_sample(frsf,50) %>% st_sf(.)
variable_de_taille <- round(rnorm(50,5,2),0)
new <- cbind(pts, variable_de_taille) # Regrouper le tout dans un objet sf
# Tracer
plot(st_geometry(frsf))
propSymbolsLayer(new, var="variable_de_taille",col="#E84923AA")
Aller plus loin dans le package cartography