install.packages(c('readxl', 'sf', 'terra', 'leaflet', 'geodata', 'wbstats',
'mapsf', 'rnaturalearth', 'osmextract', 'osrm', 'plotly'))GEO UNIV’R Tunisie 2024
LUN3 - Acquisition de données géographiques et visualisations de base
Note
Ce support est inspiré en de nombreux points du manuel Géomatique avec R (Giraud et Pecout 2024). Et on les en remercie chaleureusement 🙏.
Ce document montre les solutions techniques pour importer, explorer, manipuler, visualiser (simplement) et exporter des données géographiques avec R :
- En partant de rien.
- Avec des tableaux de données comprenant une dimension géographique (identifiant géographique).
- Avec des couches d’information géographique (vecteur et raster).
- En utilisant des packages interfaçant des API (Application Programming Interface).
- En mobilisant des packages interfaçant OpenStreetMap (OSM).
Packages utilisés
readxl: importer des fichiers Excel.plotly: créer des graphiques interactifs.sf: importer, manipuler et exporter des données géographiques vectorielles.terra: importer, manipuler et exporter des fichiers raster.leaflet: réaliser une carte interactive.geodata: accéder à des jeux de données géographiques de référence dans le monde.rnaturalearth: accéder à des fonds de carte du Monde.wbstats: utiliser les jeux de données de la Banque Mondiale via son API.osmextract: télécharger des données OpenStreetMap.osrm: calculer des temps de parcours routiers via l’engin de routage OSRM.mapsf: réaliser des cartes thématiques (présenté plus en détail dans une autre session de formation).
Si vous n’avez pas ces packages, vous pouvez les installer en exécutant la ligne suivante dans la console.
Données mobilisées
Télécharger les données
Les exemples proposés dans ce support de formation sont tantôt adaptés au contexte tunisien ou au continent africain. Elles sont accessibles ici. Le dossier /data_intro du dépôt est organisé comme suit.
“From scratch”
En partant de rien et avec quelques lignes de code on peut créer un jeu de données spatial pour en faire, par exemple, une carte interactive.
Trois vecteurs sont créés, puis transformés en data.frame et enfin convertis en objet spatial avec la fonction st_as_sf() du package sf. Cela permet de créer rapidement une carte interactive avec le package leaflet.
# 3 vecteurs avec nom, longitude et latitude
name <- c("Sentido Bellevue Park", "Université de Sfax", "Université Paris Cité")
long <- c(10.579, 10.742, 2.382)
lat <- c(35.913, 34.736, 48.827)
# Transformer en data.frame
loc <- data.frame(name, long, lat)
# Transformer en objet spatial
library(sf)
loc <- st_as_sf(loc,
coords = c("long", "lat"),
crs = 4326)
# Visualiser simplement avec leaflet
library(leaflet)
leaflet(loc) |>
addTiles() |>
addCircleMarkers(popup = loc$name) # HôpitauxTableaux de données
Import .csv
Import des données des délégations au format csv. On lui donne le nom de del_df.
del_df <- read.csv("data/data_intro/tun/data/don_del.csv", sep = ";", dec = ",")Avec R, il est très important de maîtriser la nature des objets importés ou transformés, et les convertir le cas échéant. La fonction class() nous permet de constater qu’il s’agit d’un data.frame. La fonction str() permet de détailler son contenu.
class(del_df)[1] "data.frame"str(del_df)'data.frame': 264 obs. of 30 variables:
$ del_code : chr "TN.AN.AR" "TN.AN.ET" "TN.AN.KA" "TN.AN.LS" ...
$ del_nom_fr : chr "Ariana Ville" "Cité Ettathamen" "kalaât El Andalous" "Soukra" ...
$ del_nom_ar : chr "أريانة المدينة" "حي التضامن" "قلعة الاندلس" "سكرة" ...
$ gou_code : chr "AN" "AN" "AN" "AN" ...
$ gou_nom : chr "Ariana" "Ariana" "Ariana" "Ariana" ...
$ gou_cap : int 1 0 0 0 0 0 0 1 0 0 ...
$ gou_cap_dist: num 1 8.8 19.8 6.7 8 8.8 15 1 9.5 7.4 ...
$ reg_code : chr "NE" "NE" "NE" "NE" ...
$ reg_nom : chr "Nord-est" "Nord-est" "Nord-est" "Nord-est" ...
$ popto_2004 : int 97687 78311 23045 89151 53752 60896 19404 32329 27977 31792 ...
$ popco_2004 : int 97687 78311 15313 89151 40176 53911 8909 32329 27977 31792 ...
$ immig_2004 : int 16961 5651 728 19129 7258 14053 1298 3963 5831 5508 ...
$ emigr_2004 : int 15426 5245 528 2832 985 1443 723 9381 1135 3991 ...
$ mobil_2004 : int 32387 10896 1256 21961 8243 15496 2021 13344 6966 9499 ...
$ menag_2004 : int 27468 11950 4709 21590 17119 14276 4215 7941 6498 8462 ...
$ ordin_2004 : int 9751 430 188 2979 685 2327 181 834 1176 1963 ...
$ porta_2004 : int 22524 4505 1865 13321 7087 9022 1775 4908 4178 6397 ...
$ telef_2004 : int 18596 3824 829 9262 7019 6339 1176 4002 2755 4747 ...
$ popto_2014 : int 114486 84312 26796 129693 89884 106414 24503 31128 40101 34962 ...
$ popco_2014 : int 114486 84312 18211 129693 58641 94961 11351 31128 40101 34962 ...
$ immig_2014 : int 15637 5028 1104 20786 14400 20128 1997 2392 5855 3957 ...
$ emigr_2014 : int 20448 6752 752 5528 1828 2521 822 8495 1883 3803 ...
$ mobil2014 : int 36085 11780 1856 26314 16228 22649 2819 10887 7738 7760 ...
$ menag_2014 : int 32498 22087 6554 33981 22781 27574 5922 8506 10066 9996 ...
$ ordin_2014 : int 25474 6836 1701 18191 8326 14284 1723 4339 5494 6471 ...
$ porta_2014 : int 32308 21715 6305 33683 22531 27338 5811 8363 9970 9901 ...
$ telef_2014 : int 19942 3496 428 9831 4249 7505 637 2918 3110 4901 ...
$ popto_2010 : num 109500 82922 24367 107829 69247 ...
$ surfa_2010 : num 18.61 3.38 188.21 27.89 22.91 ...
$ idr_2011 : num 0.638 0.386 0.383 0.557 0.466 0.531 0.358 0.489 0.51 0.523 ...On visualise les premières lignes avec la fonction head().
head(del_df, 5) del_code del_nom_fr del_nom_ar gou_code gou_nom gou_cap
1 TN.AN.AR Ariana Ville أريانة المدينة AN Ariana 1
2 TN.AN.ET Cité Ettathamen حي التضامن AN Ariana 0
3 TN.AN.KA kalaât El Andalous قلعة الاندلس AN Ariana 0
4 TN.AN.LS Soukra سكرة AN Ariana 0
5 TN.AN.MN Mnihla المنيهلة AN Ariana 0
gou_cap_dist reg_code reg_nom popto_2004 popco_2004 immig_2004 emigr_2004
1 1.0 NE Nord-est 97687 97687 16961 15426
2 8.8 NE Nord-est 78311 78311 5651 5245
3 19.8 NE Nord-est 23045 15313 728 528
4 6.7 NE Nord-est 89151 89151 19129 2832
5 8.0 NE Nord-est 53752 40176 7258 985
mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014
1 32387 27468 9751 22524 18596 114486 114486
2 10896 11950 430 4505 3824 84312 84312
3 1256 4709 188 1865 829 26796 18211
4 21961 21590 2979 13321 9262 129693 129693
5 8243 17119 685 7087 7019 89884 58641
immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014
1 15637 20448 36085 32498 25474 32308 19942
2 5028 6752 11780 22087 6836 21715 3496
3 1104 752 1856 6554 1701 6305 428
4 20786 5528 26314 33981 18191 33683 9831
5 14400 1828 16228 22781 8326 22531 4249
popto_2010 surfa_2010 idr_2011
1 109500 18.612 0.638
2 82922 3.376 0.386
3 24367 188.206 0.383
4 107829 27.895 0.557
5 69247 22.907 0.466Import Excel
La package readxl permet l’import de fichiers Excel. La fonction read_excel a plusieurs arguments utiles pour spécifier en entrée le format des colonnes (texte, numérique), ne pas considérer les n premières lignes du fichier, etc.
Nous importons ici un fichier Excel dérivé du World Population Prospects des Nations Unies, qui met à disposition toute sorte d’indicateurs démographiques (population, natalité, mortalité). On dispose ici des dénombrements de population par tranche d’âge quinquennale de 1950 à 2021 (feuille estimates).
library(readxl)
world <- read_excel("data/data_intro/world/data/unpp_POP_5y.xlsx",
sheet = "Estimates",
skip = 16,
col_types = c(rep("text", 10), rep("numeric", 22))) Manipulation de data.frame
Dans les faits ces fichiers nécessitent généralement d’être reformatés et réorganisés pour pouvoir être interprétés dans un logiciel quel qu’il soit… Commençons par comprendre la structure des données : noms de colonnes, années disponibles etc.
# Conversion en data.frame
world <- data.frame(world)
# Nom des colonnes
colnames(world) [1] "Index"
[2] "Variant"
[3] "Region..subregion..country.or.area.."
[4] "Notes"
[5] "Location.code"
[6] "ISO3.Alpha.code"
[7] "ISO2.Alpha.code"
[8] "SDMX.code.."
[9] "Type"
[10] "Parent.code"
[11] "Year"
[12] "X0.4"
[13] "X5.9"
[14] "X10.14"
[15] "X15.19"
[16] "X20.24"
[17] "X25.29"
[18] "X30.34"
[19] "X35.39"
[20] "X40.44"
[21] "X45.49"
[22] "X50.54"
[23] "X55.59"
[24] "X60.64"
[25] "X65.69"
[26] "X70.74"
[27] "X75.79"
[28] "X80.84"
[29] "X85.89"
[30] "X90.94"
[31] "X95.99"
[32] "X100." # Années disponibles
unique(world$Year) [1] 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964
[16] 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979
[31] 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994
[46] 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009
[61] 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 NAOn peut alors filtrer les dimensions du tableau, calculer des indicateurs, etc. qui pourront être utilisés pour des statistiques de base, des graphiques voire de futures représentations cartographiques.
# Somme de plusieurs colonnes
world$POP_TOT <- rowSums(world[,c(12:32)]) # Population totale
# Renommer une colonne
names(world)[3] <- "Name"
# Créer des indicateurs
# Part des jeunes
world$YOUNG_RT <- (world$X0.4 + world$X5.9 + world$X10.14) / world$POP_TOT * 100
# Part des personnes âgées
world$OLD_RT <- (world$X65.69 + world$X70.74 + world$X75.79 + world$X80.84 +
world$X85.89 + world$X90.94 + world$X95.99) / world$POP_TOT * 100
# Sélection d'une ligne
world_2021 <- world[world$Year == 2021,] # Année de référence
# Sélectionner plusieurs lignes
afr <- c("910", "911", "912", "913", "914")
afr_2021 <- world_2021[world_2021$Parent.code %in% afr,] # Pays africains
# Ordonner selon les valeurs d'une colonne
afr_2021 <- afr_2021[order(afr_2021$YOUNG_RT, decreasing = TRUE),]Il existe de nombreuses fonctions pour l’analyse statistique avec R. La plus basique étant probablement summary().
# Résumé stat
summary(afr_2021) Index Variant Name Notes
Length:58 Length:58 Length:58 Length:58
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
Location.code ISO3.Alpha.code ISO2.Alpha.code SDMX.code..
Length:58 Length:58 Length:58 Length:58
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
Type Parent.code Year X0.4
Length:58 Length:58 Min. :2021 Min. : 0.22
Class :character Class :character 1st Qu.:2021 1st Qu.: 299.04
Mode :character Mode :character Median :2021 Median : 1971.23
Mean :2021 Mean : 3576.71
3rd Qu.:2021 3rd Qu.: 4280.10
Max. :2021 Max. :34830.80
X5.9 X10.14 X15.19 X20.24
Min. : 0.275 Min. : 0.26 Min. : 0.197 Min. : 0.166
1st Qu.: 287.524 1st Qu.: 260.32 1st Qu.: 231.017 1st Qu.: 218.739
Median : 1799.815 Median : 1596.25 Median : 1340.744 Median : 1073.518
Mean : 3237.939 Mean : 2880.72 Mean : 2481.160 Mean : 2123.874
3rd Qu.: 3812.367 3rd Qu.: 3456.02 3rd Qu.: 3054.270 3rd Qu.: 2727.918
Max. :30567.689 Max. :26974.49 Max. :22929.125 Max. :18806.200
X25.29 X30.34 X35.39 X40.44
Min. : 0.252 Min. : 0.244 Min. : 0.287 Min. : 0.267
1st Qu.: 201.104 1st Qu.: 188.732 1st Qu.: 152.122 1st Qu.: 122.764
Median : 960.347 Median : 880.112 Median : 738.572 Median : 573.226
Mean : 1858.879 Mean : 1643.911 Mean : 1421.126 Mean :1171.518
3rd Qu.: 2325.300 3rd Qu.: 1923.398 3rd Qu.: 1581.452 3rd Qu.:1374.119
Max. :15675.790 Max. :13171.290 Max. :11591.915 Max. :9893.223
X45.49 X50.54 X55.59 X60.64
Min. : 0.358 Min. : 0.486 Min. : 0.452 Min. : 0.496
1st Qu.: 99.267 1st Qu.: 84.229 1st Qu.: 67.378 1st Qu.: 54.439
Median : 469.567 Median : 374.635 Median : 291.823 Median : 233.514
Mean : 938.585 Mean : 770.335 Mean : 616.303 Mean : 473.997
3rd Qu.:1124.716 3rd Qu.: 851.100 3rd Qu.: 667.341 3rd Qu.: 534.762
Max. :7846.818 Max. :6096.186 Max. :4876.240 Max. :3778.750
X65.69 X70.74 X75.79 X80.84
Min. : 0.399 Min. : 0.4345 Min. : 0.3185 Min. : 0.167
1st Qu.: 40.704 1st Qu.: 26.6259 1st Qu.: 16.2689 1st Qu.: 9.041
Median : 156.874 Median : 94.0440 Median : 58.9380 Median : 29.567
Mean : 349.206 Mean : 235.3130 Mean : 140.6908 Mean : 70.866
3rd Qu.: 377.993 3rd Qu.: 237.7533 3rd Qu.: 132.5434 3rd Qu.: 68.634
Max. :2757.823 Max. :1819.8435 Max. :1102.5085 Max. :494.341
X85.89 X90.94 X95.99 X100.
Min. : 0.079 Min. : 0.038 Min. : 0.0075 Min. :0.000000
1st Qu.: 4.053 1st Qu.: 0.968 1st Qu.: 0.0685 1st Qu.:0.004375
Median : 9.991 Median : 2.052 Median : 0.3372 Median :0.031750
Mean : 28.099 Mean : 7.973 Mean : 1.4798 Mean :0.219655
3rd Qu.: 27.383 3rd Qu.: 7.602 3rd Qu.: 1.2139 3rd Qu.:0.152125
Max. :178.652 Max. :81.603 Max. :22.9160 Max. :3.797500
POP_TOT YOUNG_RT OLD_RT
Min. : 5.4 Min. :13.97 Min. : 1.677
1st Qu.: 2388.4 1st Qu.:34.33 1st Qu.: 2.696
Median : 12774.0 Median :40.41 Median : 3.152
Mean : 24028.9 Mean :38.12 Mean : 4.298
3rd Qu.: 28556.3 3rd Qu.:43.45 3rd Qu.: 4.161
Max. :213401.3 Max. :48.90 Max. :26.714
L’intérêt de la programmation pour manipuler des données
Les fournisseurs de données institutionnels distribuent parfois des tableaux de données peu adaptés à l’intégration dans un logiciel quel qu’il soit. La mise en forme des données dans un langage de programmation permet :
- d’éviter d’avoir à manipuler ces fichiers manuellement, source d’erreurs.
- d’avoir un façon générique et rapide d’importer un grand nombre de tableaux de données.
Voici le fichier d’origineproduit par l’Institut National de la Statistique tunisien (INS) qui présente le nombre de chômeurs par sexe et niveau d’instruction par gouvernorat tunisien (INS) :
Un court bloc de code permet de réorganiser le fichier comme désiré, et pourrait être utilement étendu aux autres fichiers distribués par l’INS.
# Import
input_gouv <- read_excel("data/data_intro/tun/data/chomage.xls", # Chemin du fichier
sheet = "Sheet1", # Nom de la feuille Excel
skip = 3, # Retirer les trois premières lignes
col_types = c(rep("text", 3), rep("numeric", 6))) # Format des colonnes
# Table de passage Recensement / codes internationaux
gouv <- read.csv("data/data_intro/tun/data/code_tun.csv")
# Reformatage
input_gouv <- data.frame(input_gouv)
input_gouv[,1:2] <- NULL # Retirer les deux premières colonnes
# Noms de colonnes
cols <- c("CHOM_INS_NON_", "CHOM_INS_PRI_", "CHOM_INS_SEC_", "CHOM_INS_SUP_",
"CHOM_INS_UNK_", "CHOM_INS_TOT_")
colnames(input_gouv)[1] <- "id_TUN"
# Chômeurs
# Hommes
tmp1 <- input_gouv[c(1:25),]
names(tmp1)[2:length(tmp1)] <- paste0(cols, "M")
gouv <- merge(gouv, tmp1, by = "id_TUN", all.x = TRUE)
# Femmes
tmp2 <- input_gouv[c(26:50),]
names(tmp2)[2:length(tmp2)] <- paste0(cols, "F")
gouv <- merge(gouv, tmp2, by = "id_TUN", all.x = TRUE)Une fois mis en forme, on peut réaliser un graphique sans avoir modifié le fichier initial provenant de l’INS.
# 2 graphiques par ligne
par(mfrow = c(1,2), mar = c(2,2,2,2))
# Boxplot hommes
boxplot(gouv$CHOM_INS_NON_M, # Chômeurs hommes par niveau d'instruction
gouv$CHOM_INS_PRI_M,
gouv$CHOM_INS_SEC_M,
gouv$CHOM_INS_SUP_M,
ylim = c(0, 12000), # Bornes min/max de l'axe des ordonnées
main = "Hommes", # Titre plot
names = c("Rien", "Primaire", "Secondaire", "Tertiaire"), # Labels (X)
ylab = "Nombre de chômeurs, 2014", # Label (Y)
col = "#adcaf7", # Couleur des box-plots
cex.axis = .6, # Taille des labels des axes (réduit de 70 %)
cex.title = .6) # Taille du label du titre (réduit de 70 %)
# Boxplot femmes
boxplot(gouv$CHOM_INS_NON_F,
gouv$CHOM_INS_PRI_F,
gouv$CHOM_INS_SEC_F,
gouv$CHOM_INS_SUP_F,
ylim = c(0, 12000),
main = "Femmes",
names = c("Rien", "Primaire", "Secondaire", "Tertiaire"),
col = "#ed9fb0",
cex.axis = .6,
cex.title = .6)
Export
write.csv() exporte un data.frame selon un chemin spécifié.
write.csv(x = gouv, # Objet à exporter
file = "data/tun/data/gouv_chom.csv", # Chemin de fichier
row.names = FALSE) # Pour retirer les numéros de ligneGraphiques de base
Les possibilités offertes en terme de représentations graphiques sont nombreuses avec R ! En une ligne de code on peut créer des représentations graphiques variées.
On s’intéresse ici à la part des jeunes (0-14 ans, YOUNG_RT) et à la part des personnes âgées (65 ans et plus, OLD_RT) des pays africains dans la population totale en 2021.
par(mar = c(2,2,2,2), mfrow = c(2, 3))
barplot(afr_2021$YOUNG_RT, main = "Diagramme en barre")
boxplot(afr_2021$YOUNG_RT, main = "Boîtes à moustache")
hist(afr_2021$YOUNG_RT, main = "Histogrammes")
hist(afr_2021$YOUNG_RT, freq = FALSE, main = "Histogrammes et densité")
lines(density(afr_2021$YOUNG_RT), col = "blue")
stripchart(afr_2021$YOUNG_RT, method = "jitter", pch = 16,
main = "Diagrammes de dispersion")
plot(data = afr_2021, YOUNG_RT ~ OLD_RT, main = "Nuages de points")
Ces graphiques sont paramétrables avec une série d’arguments graphiques.
par(mar = c(4,4,0,4))
# Line plot
tun <- world[world$Name == "Tunisia",]
alg <- world[world$Name == "Algeria",]
mar <- world[world$Name == "Morocco",]
plot(alg$Year, # Abscisses
alg$YOUNG_RT, # Ordonnées
type = "l", # Type lignes
ylim = c(20, 50), # Bornes min/max des ordonnées
cex = .6, # Taille des points
col = "blue", # Couleur de la ligne
cex.lab = 0.7, # Réduit les labels d'un facteur de 0.7
cex.axis = 0.6, # Réduit les labels des graduations d'un facteur de 0.6
xlab = "Années", # Label abscisses
ylab = "Part des jeunes (%) dans la population totale") # Label ordonnées
lines(mar$Year,
mar$YOUNG_RT, # Rajouter une ligne (Maroc)
col = "darkgreen",
type = "l",
cex = .6)
lines(tun$Year, # Et Tunisie
tun$YOUNG_RT,
type = "l",
cex = .6,
col = "red")
# Organisation de la légende
legend(x = "topright",
legend = c("Algérie", "Maroc", "Tunisie"),
col = c("blue", "darkgreen", "red"), lty = 1,
cex = .8)
Graphiques interactifs
Le package plotly (Sievert et al. 2024) permet d’associer aux graphiques une dimension interactive. Pour utiliser les fonctions de ce package, il faut bien avoir en tête le format de données attendu (format long).
library(plotly)
# Importer jeu de données d'exemple de gapminder
gap <- read.csv("data/data_intro/world/data/gapminder.csv")
head(gap) country continent year lifeExp pop gdpPercap
1 Afghanistan Asia 1952 28.801 8425333 779.4453
2 Afghanistan Asia 1957 30.332 9240934 820.8530
3 Afghanistan Asia 1962 31.997 10267083 853.1007
4 Afghanistan Asia 1967 34.020 11537966 836.1971
5 Afghanistan Asia 1972 36.088 13079460 739.9811
6 Afghanistan Asia 1977 38.438 14880372 786.1134# Paramétrage des cercles proportionnels
area_max <- 2000 # Diamètre maximal
area_min <- area_max/(max(gap$pop)/min(gap$pop)) # Diamètre minimal
# Créer un graphique de type "Multiple Trace Animations"
fig <- gap %>%
plot_ly(
x = ~gdpPercap, # Variable X
y = ~lifeExp, # Variable Y
size = ~pop, # Taille des cercles
color = ~continent, # Couleur des cercles
frame = ~year, # Animation temporelle
sizes = c(area_min, area_max), # Gestion taille cercles
marker = list(opacity = 0.5, sizemode = 'area'), # Taille et opacité
text = ~country, # Nom du champ qui s'affichera interactivement
hoverinfo = "text",
type = 'scatter', # Type de graphique
mode = 'markers' # Mode de représentation des figurés
)
# Paramétrer l'échelle, les labels, etc.
fig <- fig %>% layout(
xaxis = list(type = "log", title = "PIB par habitant"),
yaxis = list(title = "Espérance de vie")
)
fig
Pour aller plus loin
r-graph-gallery présente les possibilités graphiques les plus courantes, en langage de base R. L’exploration de la syntaxe ggplotet son package de référence ggplot2 peut être utile pour la recherche de représentations graphiques avancées. Cette syntaxe est un peu différente du langage de base R et repose sur les principes de la “grammaire des graphiques”. Plusieurs manuels très bien construits permettent de rentrer dans l’univers de la visualisation de données avec ggplot2, comme Modern Data Visualization with R (Kabacoff 2023)
ggplot2 et ggallivial (R. Ysebaert, 2024)Données vectorielles
Le package sf permet d’importer et manipuler des couches géographiques vectorielles (points, lignes, polygones).
Import
Le geopackage est un format de données ouvert qui permet de stocker plusieurs couches géographiques dans un même fichier. La fonction st_layers() permet d’avoir un aperçu des couches présentes dans un geopackage.
library(sf)
st_layers("data/data_intro/tun/geom/tun_admin.gpkg")Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 delegation Multi Polygon 266 11 ETRS89-extended / LAEA Europe
2 gouvernorat Multi Polygon 24 4 ETRS89-extended / LAEA Europe
3 region Multi Polygon 6 2 ETRS89-extended / LAEA EuropeImporter des couches géographiques avec st_read().
del <- st_read("data/data_intro/tun/geom/tun_admin.gpkg", layer = "delegation", quiet = TRUE)
gouv <- st_read("data/data_intro/tun/geom/tun_admin.gpkg", layer = "gouvernorat", quiet = TRUE)
reg <- st_read("data/data_intro/tun/geom/tun_admin.gpkg", layer = "region", quiet = TRUE)Il s’agit bien d’objets sf.
class(del)[1] "sf" "data.frame"La fonction st_read() peut aussi être employée pour des formats de fichiers .geojson, .shapefiles, etc.
reg <- st_read("data/data_intro/tun/geom/map_reg.geojson", quiet = TRUE) Les objets sf sont des data.frame dont l’une des colonnes contient des géométries. Cette colonne est de la classe sfc (simple feature column) et chaque individu de la colonne est un sfg (simple feature geometry).
head(del)Simple feature collection with 6 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 4184398 ymin: 1282883 xmax: 4380989 ymax: 1539573
Projected CRS: ETRS89-extended / LAEA Europe
del_code del_nom_fr del_code_riate del_code_ins del_nom_ar gou_code
1 TN.SF.AG Agareb TS2146 TN347 عقارب TN.SF
2 TN.JE.AD Ain Drahem TS1224 TN225 عين دراهم TN.JE
3 TN.SS.AK Akouda TS2115 TN316 اكودة TN.SS
4 TN.KR.AL Alaa TS2216 TN417 العلاء TN.KR
5 TN.BJ.AM Amdoun TS1212 TN213 عمدون TN.BJ
6 TN.AN.AR Ariana Ville TS1120 TN121 أريانة المدينة TN.AN
gou_nom gou_cap gou_cap_dist reg_code reg_nom
1 Sfax 0 30.1 CE Centre-est
2 Jendouba 0 38.5 NO Nord-ouest
3 Sousse 0 14.3 CE Centre-est
4 Kairouan 0 48.6 CO Centre-ouest
5 Beja 0 18.9 NO Nord-ouest
6 Ariana 1 1.0 NE Nord-est
geom
1 MULTIPOLYGON (((4375308 130...
2 MULTIPOLYGON (((4224769 152...
3 MULTIPOLYGON (((4373533 142...
4 MULTIPOLYGON (((4306097 141...
5 MULTIPOLYGON (((4243838 153...
6 MULTIPOLYGON (((4338370 153...Affichage
Aperçu des variables avec plot() :
plot(del)
Afficher juste les géométries:
# 3 cartes par ligne
par(mfrow = c(1,3), mar = c(2,2,2,2))
# Délégations
plot(del$geom, # Géométries uniquement
col = "peachpuff", # Couleur du fond
border = "white", # Couleur de bordure
main = "Délégations") # Titre
# Gouvernorats
plot(gouv$geom,
col = "peachpuff",
border = "white",
main = "Gouvernorats")
# "Régions"
plot(reg$geom,
col = "peachpuff",
border = "white",
main = "Régions")
Manipulations de base
Jointures attributaires
Nous pouvons joindre un data.frame à un objet sf en utilisant la fonction merge() et en s’appuyant sur des identifiants communs aux deux objets.
Attention à l’ordre des arguments, l’objet retourné sera du même type que x. Il n’est pas possible de faire une jointure attributaire en utilisant deux objets sf.
del <- merge(
x = del[,"del_code"], # L'objet sf (seulement le champ del_code)
y = del_df, # le data.frame
by.x = "del_code", # identifiant dans x
by.y = "del_code", # identifiant dans y
all.x = TRUE # conserver toutes les lignes
)
head(del)Simple feature collection with 6 features and 30 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 4321914 ymin: 1528813 xmax: 4346090 ymax: 1562156
Projected CRS: ETRS89-extended / LAEA Europe
del_code del_nom_fr del_nom_ar gou_code gou_nom gou_cap
1 TN.AN.AR Ariana Ville أريانة المدينة AN Ariana 1
2 TN.AN.ET Cité Ettathamen حي التضامن AN Ariana 0
3 TN.AN.KA kalaât El Andalous قلعة الاندلس AN Ariana 0
4 TN.AN.LS Soukra سكرة AN Ariana 0
5 TN.AN.MN Mnihla المنيهلة AN Ariana 0
6 TN.AN.RA Raoued رواد AN Ariana 0
gou_cap_dist reg_code reg_nom popto_2004 popco_2004 immig_2004 emigr_2004
1 1.0 NE Nord-est 97687 97687 16961 15426
2 8.8 NE Nord-est 78311 78311 5651 5245
3 19.8 NE Nord-est 23045 15313 728 528
4 6.7 NE Nord-est 89151 89151 19129 2832
5 8.0 NE Nord-est 53752 40176 7258 985
6 8.8 NE Nord-est 60896 53911 14053 1443
mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014
1 32387 27468 9751 22524 18596 114486 114486
2 10896 11950 430 4505 3824 84312 84312
3 1256 4709 188 1865 829 26796 18211
4 21961 21590 2979 13321 9262 129693 129693
5 8243 17119 685 7087 7019 89884 58641
6 15496 14276 2327 9022 6339 106414 94961
immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014
1 15637 20448 36085 32498 25474 32308 19942
2 5028 6752 11780 22087 6836 21715 3496
3 1104 752 1856 6554 1701 6305 428
4 20786 5528 26314 33981 18191 33683 9831
5 14400 1828 16228 22781 8326 22531 4249
6 20128 2521 22649 27574 14284 27338 7505
popto_2010 surfa_2010 idr_2011 geom
1 109500 18.612 0.638 MULTIPOLYGON (((4338370 153...
2 82922 3.376 0.386 MULTIPOLYGON (((4330561 153...
3 24367 188.206 0.383 MULTIPOLYGON (((4333853 156...
4 107829 27.895 0.557 MULTIPOLYGON (((4342368 153...
5 69247 22.907 0.466 MULTIPOLYGON (((4332840 153...
6 77419 57.080 0.531 MULTIPOLYGON (((4338691 154...plot(del[,"idr_2011"])
Sélectionner des lignes, des colonnes
Les objets sf sont des data.frame, on peut donc sélectionner leur lignes et leur colonnes de la même manière.
# Sélection de lignes
sou <- del[del$gou_nom == "Sousse",]
# Sélection de colonnes
sou <- sou[,"idr_2011"]
# Ne conserver que les lignes avec une valeur
sou <- sou[!is.na(sou$idr_2011),]
plot(sou[,"idr_2011"])
Export
st_write("data/data_intro/tun/geom/sousse_deleg.geojson")Données raster
Le package terra (Hijmans 2024) permet d’importer et d’exporter des fichiers raster. Il repose sur la bibliothèque GDAL (GDAL/OGR contributors 2022) qui permet de lire et de traiter un très grand nombre de format d’images géographiques.
Import
La fonction rast() permet de créer et/ou d’importer des données raster. Nous importons ici un jeu de données au format .tif créé et distribué par WorldPop qui porte sur une estimation de la population tunisienne dans une résolution d’1 km (aussi accessible à 100 mètres).
library(terra)
pop <- rast("data/data_intro/tun/geom/pop.tif")Ce sont des objets de type SpatRaster.
popclass : SpatRaster
dimensions : 879, 488, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : 7.532083, 11.59875, 30.24125, 37.56625 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : pop.tif
name : pop
min value : 0.00
max value : 23068.07 La fonction summary() est toujours utile pour un résumé statistique des cellules. Vu le nombre important de cellules, ce résumé est effectué sur un échantillon.
summary(pop) pop
Min. : 0.00
1st Qu.: 0.95
Median : 9.16
Mean : 53.95
3rd Qu.: 25.72
Max. :18878.24
NA's :49070 Affichage
Par défaut et comme montré précédemment, la fonction plot() renvoie une légende continue sur les représentations cartographiques. On peut choisir de discrétiser cette information pour représenter des classes de valeur associée à une palette de couleurs personnalisée.
par(mfrow = c(1,2))
# Graphique en échelle continue
plot(pop)
# Avec discrétisation et paramétrage des couleurs
cuts <- c(0, 5, 25, 50, 250, 500, 2000, 18879)
cols <- colorRampPalette(c("yellow", "darkgoldenrod1", "brown2"))
plot(pop, breaks = cuts, col = cols(7))
Manipulation de base
Le découpage d’un raster en fonction de l’étendue d’un autre objet, SpatVector ou SpatRaster, est réalisable avec la fonction crop(). Les deux couches de données doivent être dans la même projection.
# Transformer en WGS 84
sou <- st_transform(sou, 4326)
# Crop avec la délégation de Sousse
pop_sou <- crop(pop, sou)
# Plot
plot(pop_sou, breaks = cuts, col = cols(8))
plot(sou$geom, col = NA, add = TRUE)
Export
La fonction writeRaster() permet d’enregistrer un objet SpatRaster.
writeRaster(x = pop_sou, filename = "data/data_intro/tun/geom/pop_sou.tif")Packages de données spatiales
De nombreux packages distribuent des données géographiques. Ils interfacent souvent des API qui permettent d’interroger des données mises à disposition sur le Web, directement avec R.
L’intérêt des API est de se connecter directement aux bases de données des fournisseurs de données, et de disposer des dernières mises à jour des données.
Principaux packages
Au niveau international
rnaturalearth(Massicotte et South 2023) : permet de récupérer les données cartographiques Natural Earth.
cshapes(Weidmann, Schvitz, et Girardin 2021) : met à disposition les frontières nationales, de 1886 à aujourd’hui.
geonames(Rowlingson 2019) : permet d’interroger la BD geonames, qui fournit notamment des localisations.wbstats(Piburn 2020) etWDI(Arel-Bundock 2022) : donnent accès aux données et statistiques de la Banque mondiale.
sen2r(Ranghetti et al. 2020) : permet de télécharger et prétraiter automatiquement les données du satellite Sentinel-2.
MODIStsp(Busetto et Ranghetti 2016) : permet de trouver, télécharger et traiter des images MODIS.geodata(Hijmans et al. 2023) : fournit un accès à des données sur le climat, l’altitude, le sol, la présence d’espèces et les limites administratives.elevatr(Hollister et al. 2023) : donne accès à des données d’élévation mises à disposition par Amazon Web Services Terrain Tiles, l’Open Topography Global Datasets API et l’USGS Elevation Point Query Service.rgee(Aybar 2023) : permet d’utiliser l’API de Google Earth Engine, catalogue de données publiques et infrastructure de calcul pour les images satellites.
nasapower(Sparks 2018) : API client NASA (prévision des ressources énergétiques mondiales, météorologie, énergie solaire de surface et climatologie).geoknife(Read et al. 2015) : permet de traiter (en ligne) des données matricielles volumineuses issues du Geo Data Portal de l’U.S. Geological Survey.rdhs(Watson, FitzJohn, et Eaton 2019) : API client et gestions de données de l’enquête démographique et de santé (DHS).giscoR(Hernangómez 2024a) : permet de télécharger des données cartographiques mondiales et européennes de la BD GISCO d’Eurostat (système d’information géographique de la Commission).eurostat(Lahti et al. 2017) : permet de télécharger des données de la BD Eurostat.
Au niveau national
- Brésil
geobr(Pereira et Goncalves 2024) : fournit un accès facile aux séries de données spatiales officielles du Brésil pour différentes années et découpage administratifs.
- Chili -
chilemapas(Vargas 2022) : donne accès aux divisions politiques et administratives du Chili. - Espagne
mapSpain(Hernangómez 2024b) : propose les limites administratives de l’Espagne à plusieurs niveaux (Communautés autonomes, Provinces, Municipalités), ainsi que des tuiles.
- États-Unis
tidycensus(Walker et Herman 2024) : permet de charger des données et géométries du recensement américain en formatsfettidyversetigris(Walker 2024b) : donne accès aux éléments cartographiques fournis par le US Census Bureau TIGER, y compris les limites cartographiques, les routes et l’eau.FedData(Bocinsky 2024) : automatise le téléchargement de données géospatiales disponibles à partir de plusieurs sources de données fédérées.acs(Glenn 2019) : permet de télécharger et manipuler les données de l’American Community Survey et les données décennales du recensement des États-Unis.censusapi(Recht 2022) : wrapper pour les API du Census Bureau des États-Unis.idbr(Walker 2024a) : interface avec l’API de la base de données internationale du US Census Bureau.
ipumsr(Greg Freedman Ellis, Derek Burk, et Finn Roberts 2024) : Permet d’importer des données de recensement, d’enquête et géographiques fournies par l’IPUMS.totalcensus(Li 2021) : permet d’extraire les données du recensement décennal et de l’American Community Survey au niveaux des block, block group et tract.
- Finland
mapsFinland(Haukka 2023) : donne un accès à des cartes et données concernant la Finland.
- France
happign(Carteron 2023) : accès à certaines données de l’IGN.insee(Leclerc 2022) : pour télécharger facilement les données de la base BDM (Banque de Données Macroéconomiques) de l’INSEE.
- Pologne
rgugik(Dyba et Nowosad 2021) : permet l’acquisition automatique de données ouvertes à partir des ressources du Bureau central polonais de la géodésie et de la cartographie (Główny Urząd Geodezji i Kartografii ).
- Uruguay
geouy(Detomasi 2023) : permet le chargement d’informations géographiques sur l’Uruguay.
- …
geodata
Ce package facilite l’accès à des données géographique de référence sur le climat, la couverture du sol, les limites administratives et plusieurs autres jeux de données de référence au niveau mondial.
Les couches géographiques inclues dans le package
geodata
| Function | Description |
|---|---|
bio_oracle |
Marine data from Bio-Oracle |
cmip6_world, cmip6_tile |
Downscaled and calibrated CMIP6 projected future climate data |
country_codes |
Country codes |
crop_calendar_sacks |
Crop calendar data by Sacks et al |
crop_monfreda |
Crop area and yield data for 175 crops by Monfreda et al. |
crop_spam |
MapSPAM crop data (area, yield, value) |
cropland |
Cropland density for the world derived from different sources (ESA, GLAD, QED) |
elevation_3s, elevation_30s, elevation_global |
Elevation data |
gadm, world |
Administrative boundaries for any country, or the entire world from GADM |
landcover |
Landcover data derived from ESA WorldCover |
footprint |
Human footprint data from the Last of the Wild project |
osm |
OpenStreetMap data by country (places and roads) |
population |
Population density Gridded Population of the World |
soil_af |
Chemical and physical soil properties data for Africa for different soil depths |
soil_af_elements |
Connect to or download chemical soil element concentration (for the 0-30 cm topsoil) data for Africa |
soil_af_water |
Physical soil properties data for Africa for water balance computation |
soil_af_isda |
Soil data for Africa derived from the iDSA data set |
soil_world_vsi |
Virtually connect to the global SoilGrids data |
soil_world |
Global soils data from SoilGrids |
sp_occurrence |
Species occurrence records from the Global Biodiversity Information Facility |
travel_time |
Travel time to and from cities and ports by Nelson et al. |
worldclim_global, worldclim_country, worldclim_tile |
WorldClim glocal climate data |
Extraction des couches géographiques d’altitude et de température en Tunisie.
library(geodata)
elev <- elevation_30s(country = "TUN", path = tempdir())
temp <- worldclim_country(country = "Tunisia",
res = 10,
var = "tavg",
path = tempdir())Ce sont des objets de type SpatRaster.
class(temp)[1] "SpatRaster"
attr(,"package")
[1] "terra"class(elev)[1] "SpatRaster"
attr(,"package")
[1] "terra"# Les afficher
par(mfrow = c(1,2))
# Altitude
cols <- colorRampPalette(c("#31ad37", "#f5f752", "#fca330", "#9c5903"))
plot(elev, main = "Altitude", col = cols(50))
# Température
cols <- colorRampPalette(c("#4575b4", "#91bfdb", "#e0f3f8",
"#fee090", "#fc8d59", "#d73027"))
plot(temp$TUN_wc2.1_30s_tavg_5, main = "Températures, Mai (1970-2000)",
col = cols(50))
Natural Earth et Banque Mondiale.
Objectif
Nous cherchons ici à créer une carte de la part de la couverture forestière en Afrique en utilisant les données de la banque mondiale et le fond de carte Natural Earth. Rien besoin de télécharger, mais il faut une connexion internet.
Plusieurs fonds de carte des pays du Monde à différents niveaux de généralisation cartographique sont disponibles avec la package rnaturalearth. Nous sélectionnons uniquement les pays africains.
library(rnaturalearth)
# Import des pays
country <- ne_countries(type = "countries", # pays
scale = "small", # niveau de généralisation
returnclass = "sf") # objet retourné
# Conversion en projection Mercator
country <- st_transform(country, crs = "EPSG:3857")
# Si on s'intéresse à l'Afrique (modèle carto)
afr <- country[country$continent == "Africa",]Le package wbstatspermet d’interroger l’API de la base de données de la Banque Mondiale. On peut faire une recherche pour trouver le nom des tables qui répondent à une requête par mots-clés avec la fonction wb_search().
library(wbstats)
wb_search("Forest.area")# A tibble: 6 × 3
indicator_id indicator indicator_desc
<chr> <chr> <chr>
1 AG.LND.FRST.HA Forest area (hectares) Forest area i…
2 AG.LND.FRST.K2 Forest area (sq. km) Forest area i…
3 AG.LND.FRST.ZS Forest area (% of land area) Forest area i…
4 ER.FST.DFST.ZG Annual deforestation (% of chan… Average annua…
5 IN.ENV.COASTALZONE.FOREST.AREA Forest (Non- tidal)/ Plantation… <NA>
6 IN.ENV.FOREST.AREA Forest Cover ('000 hectares) <NA>
Avertissement
Chaque package reposant sur une API a son fonctionnement propre. Il est donc conseillé de consulter attentivement la documentation associée à ces packages (souvent précise) pour apprendre à les utiliser. Par ailleurs les validité des extractions réalisées est dépendante des fournisseurs de données. Si le contenu de la base de donnée évolue (organisation, nom des indicateurs, etc.), le code R sur lequel repose l’extraction des données peut ne plus s’exécuter correctement avec le temps. Pour éviter ce type de désagrément, il peut être utile de sauvegarder le résultat de l’export des données.
La fonction wb_data() extrait la table de la banque Mondiale qui nous intéresse.
# Sélection des indicateurs
wb <- wb_data("AG.LND.FRST.ZS", return_wide = FALSE)
head(wb)# A tibble: 6 × 11
indicator_id indicator iso2c iso3c country date value unit obs_status
<chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
1 AG.LND.FRST.ZS Forest area (… AF AFG Afghan… 2023 NA <NA> <NA>
2 AG.LND.FRST.ZS Forest area (… AF AFG Afghan… 2022 1.85 <NA> <NA>
3 AG.LND.FRST.ZS Forest area (… AF AFG Afghan… 2021 1.85 <NA> <NA>
4 AG.LND.FRST.ZS Forest area (… AF AFG Afghan… 2020 1.85 <NA> <NA>
5 AG.LND.FRST.ZS Forest area (… AF AFG Afghan… 2019 1.85 <NA> <NA>
6 AG.LND.FRST.ZS Forest area (… AF AFG Afghan… 2018 1.85 <NA> <NA>
# ℹ 2 more variables: footnote <chr>, last_updated <date>La fonction reshape() est utilisée pour transformer le data.frame dans un format plus conforme à nos habitudes (long vers large) : les pays dans la première colonne, puis la série temporelle disponible en colonnes.
# Ordonner par année
wb <- wb[order(wb$date),]
# Sélectionner les colonnes d'intérêt
wb <- data.frame(wb[,c("iso3c", "date", "value")])
# Reformater au format long
wb <- reshape(wb,
idvar = "iso3c", # identifiants
timevar = "date", # variable contenant les dates
direction = "wide") # format long
# Résultat pour quelques colonnes
head(wb[,c("iso3c", "value.2020", "value.2021", "value.2022", "value.2023")]) iso3c value.2020 value.2021 value.2022 value.2023
1 AFG 1.852782 1.8527820 1.852782 NA
2 ALB 28.791971 28.7919708 28.791971 NA
3 DZA 0.818309 0.8222276 0.826333 NA
4 ASM 85.650000 85.5000000 85.350000 NA
5 AND 34.042553 34.0425532 34.042553 NA
6 AGO 53.426951 52.9817224 52.536497 NANous cherchons l’année la plus récente pour laquelle nous disposons d’un maximum de valeurs. Cette ligne de code permet de dénombrer le nombre de NA pour chaque colonne du jeu de données. 2021 est l’année la plus récente avec le moins de valeurs manquantes.
# Regarder le nombre de données manquantes
sapply(wb, function(x) sum(length(which(is.na(x))))) iso3c value.1960 value.1961 value.1962 value.1963 value.1964 value.1965
0 217 217 217 217 217 217
value.1966 value.1967 value.1968 value.1969 value.1970 value.1971 value.1972
217 217 217 217 217 217 217
value.1973 value.1974 value.1975 value.1976 value.1977 value.1978 value.1979
217 217 217 217 217 217 217
value.1980 value.1981 value.1982 value.1983 value.1984 value.1985 value.1986
217 217 217 217 217 217 217
value.1987 value.1988 value.1989 value.1990 value.1991 value.1992 value.1993
217 217 217 39 35 16 12
value.1994 value.1995 value.1996 value.1997 value.1998 value.1999 value.2000
12 12 12 12 12 12 10
value.2001 value.2002 value.2003 value.2004 value.2005 value.2006 value.2007
10 10 10 10 10 8 8
value.2008 value.2009 value.2010 value.2011 value.2012 value.2013 value.2014
8 8 8 5 3 3 3
value.2015 value.2016 value.2017 value.2018 value.2019 value.2020 value.2021
3 3 3 3 3 3 3
value.2022 value.2023
7 217 Nous réalisons pour une finir une jointure attributaire avec le fond de carte précédemment importé avec rnaturalearth. La cartographie est réalisée avec la fonction mf_map du package mapsf, dont le fonctionnement sera présenté ultérieurement dans la formation.
# Jointure attributaire
afr <- merge(x = afr[,c("adm0_a3", "name_fr")], # fond de carte
wb, # jeu de données de la banque mondiale
by.x = "adm0_a3", # champ de jointure fond de carte
by.y = "iso3c", # champ de jointure jeu de données
all.x = TRUE) # conserver toutes les géométries
# Carte avec la librairie map sf
library(mapsf)
mf_map(afr, # objet sf
var = "value.2021", # variable
type = "choro", # type choroplèthe
breaks = "quantile", # méthode de discrétisation
nbreaks = 5, # nombre de classes
leg_title = "Part de la surface totale (%)", # titre de légende
leg_adj = c(0, 2)) # ajustement position légende
mf_layout(title = "Surface forestière 2021", # titre de la carte
credits = "Source: Banque Mondiale, 2024", # sources
scale = FALSE) # pas d'échelle
OpenStreetMap
OpenStreetMap (OSM) est un projet de cartographie participative qui a pour but de constituer une base de données géographique libre à l’échelle mondiale. OSM vous permet de voir, modifier et utiliser des données géographiques dans le monde entier.
Plusieurs packages permettre d’extraire, interroger et visualiser des données issues d’OSM :
osmextract(Gilardi et Lovelace 2023) : permet d’importer des données OSM.
osmdata(Mark Padgham et al. 2017) : pour télécharger et utiliser les données d’OSM, en utilisant l’API Overpass turbo.maptiles(Giraud 2024) : Ce package télécharge, compose et affiche des tuiles à partir d’un grand nombre de fournisseurs (OpenStreetMap, Stadia, Esri, CARTO ou Thunderforest…).
nominatimlite(Hernangómez et Lacko 2024) ettidygeocoder(Cambon et al. 2021) : Des géocodeurs qui utilisent OSM, notamment.
Objectif
Nous cherchons à extraire d’OSM pour la Tunisie les géométries suivantes :
- Découpages administratifs : Délégations et secteurs.
- Équipements : objets décrits avec une clé OSM amenity (équipements utiles et importants).
Ces données nous permettront ensuite de calculer les temps de parcours routiers des centres de secteur vers l’hôpital ou la clinique la plus accessible.
Extraction et mise en forme
Le package osmextract (Gilardi et Lovelace 2023) permet d’extraire de télécharger des données OSM pour une zone géographique donnée. Cela permet aussi d’éviter de surcharger un serveur Overpass turbo (utilisé par le package osmdata) et ainsi de travailler sur des volumes de données plus importants.
Découpages administratifs
La fonction oe_get() permet de télécharger un extrait de la base de données OSM pour une zone particulière et un type d’objet géographique. L’argument place correspond au nom du fichier *.pbf accessible sur le site Geofabrik. L’argument extra_tag permet de sélectionner les objets de la base OSM correspondant à une clé particulière (se référer à la documentation d’OSM pour choisir les clés).
On commence par extraire tous les polygones et les points OSM inclus en Tunisie. L’import d’extra_tags permet d’obtenir des informations supplémentaires qui serviront à filtrer la base de données.
- N.B. : L’extraction peut durer quelques minutes, même avec une bonne connexion
library(osmextract)
osm_poly <- oe_get(place = "Tunisia",
layer = "multipolygons",
extra_tags = c("amenity", "ref:tn:hasc_2", "ref:tn:codegeo",
"name:fr", "name:ar"))
osm_pt <- oe_get(place = "Tunisia",
layer = "points",
extra_tags = c("amenity", "name:fr", "name:ar"))On peut ensuite filtrer l’extraction en fonction des objets / champs d’intérêt.
# Délégations et secteurs
admin <- osm_poly[!is.na(osm_poly$admin_level),] # Retirer pas d'attribut de niveau hiérarchique
deleg <- admin[admin$admin_level == 5,] # Délégations
sect <- admin[admin$admin_level == 6,] # Secteurs
# Ne garder que les champs utiles
deleg <- deleg[,c("osm_id", "ref_tn_codegeo", "ref_tn_hasc_2", "name_ar",
"name_fr", "admin_level")]
sect <- sect[,c("osm_id", "ref_tn_codegeo", "name_ar",
"name_fr", "admin_level")]Ce bloc de code, dont on ne détaillera pas le contenu ici, permet d’harmoniser le nom des champs entre les couches géographiques et de disposer des délégations et gouvernorats d’appartenance des couches géographiques des secteurs et délégations.
# Délégation d'appartenance du secteur
sect_pt <- st_centroid(sect) # Centroide du secteur
sect_pt <- st_intersection(sect_pt, deleg) # Intersection avec couche délégation
sect_pt <- st_set_geometry(sect_pt, NULL) # Retirer géométries
sect <- merge(sect, # Enrichir les secteurs du code d'appartenance de la délégation
sect_pt[,c("osm_id", "ref_tn_hasc_2")],
by = "osm_id",
all.x = TRUE)
sect <- sect[!duplicated(sect$osm_id),] # Retirer valeurs dupliquées
# Renommer colonnes
names(sect)[2] <- "id_tn"
names(sect)[6] <- "id_hasc_deleg"
sect$id_hasc_gouv <- substr(sect$id_hasc_deleg, 1, 5) # Gouvernorat d'appartenance
names(sect)[2,6] <- (c("tn_codegeo", "int_codegeo"))
deleg$id_hasc_gouv <- substr(deleg$id_hasc_deleg, 1, 5)Équipements
Les objets décrits par la clé amenity sont de nature hétérogènes (points, polygones). Ce bloc de code extrait les centroides des polygones et les associent à la couche de points initiale.
Ces couches géographiques sont finalement exportées dans un geopackage tun_osm.gpkg.
# Consolidation des géométries
sel_poly <- osm_poly[!is.na(osm_poly$amenity),] # Retirer les objets qui n'ont pas de clé amenity
sel_pt <- osm_pt[!is.na(osm_pt$amenity),]
sel_poly <- st_make_valid(sel_poly) # Consolider les géométries des polygones
sel_pt2 <- st_centroid(sel_poly) # extraire le centroide
cols <- intersect(names(sel_pt), names(sel_pt2)) # Garder les colonnes identiques
sel_pt <- rbind(sel_pt[,cols, sel_pt2[,cols]]) # Combiner points et polygones
sel_pt <- sel_pt[,c("osm_id", "amenity", "name_ar", "name_fr")]
# Exporter les couche ainsi créées
st_write(sel_pt, "data/tun/geom/tun_osm.gpkg", layer = "poi")
st_write(deleg, "data/tun/geom/tun_osm.gpkg", layer = "deleg")
st_write(sect, "data/tun/geom/tun_osm.gpkg", layer = "sect")
Pour aller plus loin
Ce document (Giraud 2017) montre comment créer un fond de carte avec des données OSM, extraire des objets d’intérêt (bars et restaurants) avec le package osmdata. Il propose également plusieurs pistes cartographiques pour manipuler et visualiser ces données. Notons simplement que la cartographie utilise les fonctions du package cartography, plus maintenu, et qu’il est conseillé d’utiliser dorénavant le package mapsf.
Résultats de l’extraction
Dénombrement des 15 objets OSM les plus fréquents décrits par la clé amenity en Tunisie.
# Import des aménités préparées en amont
osm_pt <- st_read("data/data_intro/tun/geom/tun_osm.gpkg", layer = "poi", quiet = TRUE)
# Nombre de points avec le tag "amenity".
tbl <- table(osm_pt$amenity)
tbl <- tbl[order(tbl, decreasing = TRUE)]
tbl[1:15]
cafe restaurant place_of_worship bank
2680 1272 1094 996
school pharmacy fast_food fuel
916 700 597 525
post_office police bar parking
327 153 135 135
doctors atm dentist
113 109 109 # On ne prend ici que les cliniques et hôpitaux
osm_pt <- osm_pt[osm_pt$amenity %in% c("clinic", "hospital"), ]Nombre de délégations, secteurs et cliniques/hôpitaux inclus dans OSM.
# Import des unités territoriales préparées en amont
deleg <- st_read("data/data_intro/tun/geom/tun_osm.gpkg", layer = "deleg", quiet = TRUE)
sect <- st_read("data/data_intro/tun/geom/tun_osm.gpkg", layer = "sect", quiet = TRUE)
# Nombre d'objets (lignes)
nrow(deleg)[1] 266nrow(sect)[1] 2104nrow(osm_pt)[1] 119Le résultat sous forme d’une carte interactive avec le package leaflet, qui sera présenté ultérieurement dans la formation.
Code
# Choisir la palette
library(leaflet)
pal <- colorFactor(palette = c("red", "gold"),
domain = c("clinic", "hospital"))
# Cartographie interactive
leaflet(osm_pt) |> # Emprise = hôpitaux
addProviderTiles("OpenStreetMap.HOT") |> # Type de tuiles chargées
addPolygons(data = sect, # Secteurs
col = "white",
fillColor = "lightgrey",
fillOpacity = 0.7,
weight = 1,
popup = paste0("<b>", sect$id_tn, "<br></b>",
sect$name_ar, "<br>", sect$name_fr),
group = "Secteurs")|>
addPolygons(data = deleg, # Délégations
col = "darkgrey",
fill = "lightgrey",
fillOpacity = 0,
weight = 1.2,
popup = paste0("<b>", deleg$id_tn, " / ",
deleg$id_hasc_deleg,"<br></b>",
deleg$name_ar, "<br>", deleg$name_fr),
group = "Délégations") |>
addCircleMarkers(radius = 4, # Hôpitaux
stroke = FALSE,
color = ~ pal(amenity),
fillOpacity = 1,
popup = paste0(osm_pt$name_ar, "<br>", osm_pt$name_fr),
group = "Cliniques et hôpitaux") |>
addLegend(pal = pal, # Légende pour différencier cliniques et hôpitaux
values = c("clinic", "hospital"),
opacity = 0.7,
title = "OSM amenity",
position = "bottomright") |>
addLayersControl(overlayGroups = c("Secteurs", "Délégations", "Cliniques et hôpitaux"),
options = layersControlOptions(collapsed = FALSE))