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

Lien extérieur : très intéressant - Opération et géotraitements.

Google Sites

Report abuse