Environnement de travail

Pour exécuter les programmes proposés par la formation, installez / mettez à jour R (version 3.4.4 minimum), R Studio. Lancez ensuite R studio et installez les packages au moyen des commandes suivantes (prévoir au moins 10 minutes).

Les packages

Installation des packages utilisés pendant la formation. Ouvrez R-Studio et tapez :

install.packages("sf")
install.packages("cartography")
install.packages("cartogram")
install.packages("readxl")
install.packages("countrycode")

Optionnel

install.packages("rmapshaper")
install.packages("osmdata")
install.packages("reshape2")
install.packages("eurostat")
install.packages("rnaturalearth")

Version des packages utilisés dans ce document

## [1] "sf (version 0.6.3), cartography (version 2.1.2), countrycode (version 1.0.0), readxl (version 1.1.0), rnaturalearth (version 0.1.0), cartogram (version 0.1.0), reshape2 (version 1.4.3), eurostat (version 3.2.9), rmapshaper (version 0.4.0), osmdata (version 0.0.7.1)"

Objectifs

• Espace d’étude : Hérault, communes
• Prise en main des fonctions spatiales de R

Téléchargement des données

Téléchargez les données ici : Télécharger

Décompressez le dossier et cliquez sur proj.Rproj. Le fichier exo1.R s’ouvre dans R-Studio.

Commandes de base

Vérifiez votre répertoire de travail courant.

getwd()

Consulter le contenu. Regarder ce qu’il y a dans le répertoire “data”

list.files()
list.files("data")

Import d’une couche SIG dans R avec le package sf. Pacakge basé sur GEOS (Geometry Engine Open Source) et GDAL (Geospatial Data Abstraction Library).

library("sf")

## Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3

communes <- st_read(dsn = "data/communes.shp", stringsAsFactors = F)

## Reading layer `communes' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/communes.shp' using driver `ESRI Shapefile'
## Simple feature collection with 343 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

Voir la table atributaire

communes

## Simple feature collection with 343 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## First 10 features:
##    INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1      34029                        BELARGA   Commune simple        415
## 2      34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3      34032                        BEZIERS   Sous-prfecture       9567
## 4      34082                    COMBAILLAUX   Commune simple        902
## 5      34108                     FRONTIGNAN Chef-lieu canton       3984
## 6      34166                      MONTBLANC   Commune simple       2716
## 7      34066                    CAZEVIEILLE   Commune simple       1619
## 8      34296                      SAUSSINES   Commune simple        630
## 9      34269        SAINT-JEAN-DE-MINERVOIS   Commune simple       3277
## 10     34012                     ARGELLIERS   Commune simple       5065
##    NOM_DEPT                       geometry
## 1   HERAULT POLYGON ((742075 6272005, 7...
## 2   HERAULT POLYGON ((770896 6289405, 7...
## 3   HERAULT POLYGON ((718759 6244261, 7...
## 4   HERAULT POLYGON ((761606 6284489, 7...
## 5   HERAULT POLYGON ((758738 6257679, 7...
## 6   HERAULT POLYGON ((728067 6248191, 7...
## 7   HERAULT POLYGON ((762395 6294296, 7...
## 8   HERAULT POLYGON ((784582 6294234, 7...
## 9   HERAULT POLYGON ((686897 6252592, 6...
## 10  HERAULT POLYGON ((756922 6287661, 7...

head(communes,3)

## Simple feature collection with 3 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##   INSEE_COM                       NOM_COMM         STATUT SUPERFICIE
## 1     34029                        BELARGA Commune simple        415
## 2     34290 SAINT-VINCENT-DE-BARBEYRARGUES Commune simple        222
## 3     34032                        BEZIERS Sous-prfecture       9567
##   NOM_DEPT                       geometry
## 1  HERAULT POLYGON ((742075 6272005, 7...
## 2  HERAULT POLYGON ((770896 6289405, 7...
## 3  HERAULT POLYGON ((718759 6244261, 7...

View(communes)

L’instruction plot permet d’afficher la couche. L’instruction st_geometry permet d’accéder à la variable définissant les géométries.

plot(st_geometry(communes))

Afficher la couche avec des parametres graphiques

plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)

Ajuster les marges

par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)

Introduction au SIG avec R / Manipuler les informations spatiales

Nous cherchons à identifier les communes qui ont une partie de leur territoire situé à moins de 20 km du centre de Sète.

Nous commençons par extraire la commune de Sète du fond de carte communal de référence.

macommune <- "SETE"
monpoly <- communes[communes$NOM_COMM == macommune,]
par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)

Calculer le centroïde de la commune de Sète.

moncentre <- st_centroid(x = monpoly)

## Warning in st_centroid.sf(x = monpoly): st_centroid assumes attributes are
## constant over geometries of x

par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)

Calculer une zone tampon

mydist <- 20000
buff <- st_buffer(x = st_geometry(moncentre), dist=mydist)
par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)

Récupérer la liste des communes

communes$buff <- st_intersects(st_geometry(communes), st_geometry(buff), sparse = FALSE)
head(communes)

## Simple feature collection with 6 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##   INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1     34029                        BELARGA   Commune simple        415
## 2     34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3     34032                        BEZIERS   Sous-prfecture       9567
## 4     34082                    COMBAILLAUX   Commune simple        902
## 5     34108                     FRONTIGNAN Chef-lieu canton       3984
## 6     34166                      MONTBLANC   Commune simple       2716
##   NOM_DEPT                       geometry  buff
## 1  HERAULT POLYGON ((742075 6272005, 7...  TRUE
## 2  HERAULT POLYGON ((770896 6289405, 7... FALSE
## 3  HERAULT POLYGON ((718759 6244261, 7... FALSE
## 4  HERAULT POLYGON ((761606 6284489, 7... FALSE
## 5  HERAULT POLYGON ((758738 6257679, 7...  TRUE
## 6  HERAULT POLYGON ((728067 6248191, 7... FALSE

Récupérer le nombre de communes situées à moins de 20 km de Sète.

nb <- dim(communes[communes$buff == T,])
nb[1]

## [1] 40

Extraction et affichage des communes

communes20km <- communes[communes$buff == TRUE,]
par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)

Nouvelle carte avec un titre (NB : La première couche affichée détermine l’emprise de la carte)

par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes20km), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
montitre <- paste0("Il y a ", nb[1], " communes situées à moins de ", mydist/1000, "km de ",tolower(macommune))
title(montitre)

Export (si besoin)

st_write(obj = communes20km, dsn = "outputs/resultat.shp")

A vous de jouer !

Reproduisez cette procédure avec la commune et la distance de votre choix.

Pour accéder au listing des communes d’Hérault

communes$NOM_COMM

Objectifs

• Thématique : Réaliser une carte des catégories socio-professionnelles.
• Espace d’étude : Occitanie, communes
• Visualisation associée : carte en figurés proportionnels (+ cartogrammes de Dorling)
• packages utilisés : sf, cartography, cartogram, readxl

Téléchargement des données

Téléchargez les données ici : Télécharger

Décompressez le dossier et cliquez sur proj.Rproj. Le fichier exo2.R s’ouvre dans R-Studio.

Préparation des données

Pour bien commencer, il est fondamental de jeter un oeil au fichier d’entrée pour comprendre son organisation et savoir comment l’importer sous R.

Ouvrez le fichier data/INSEE/pop-act2554-csp-cd-6814.xls Nous nous intéressons ici à l’onglet “COM_2014”.

library("sf")
library("readxl")

Importez correctement cette feuille du tableau Excel dans R. Les commandes de base ci-dessous permettent de visualiser le nom des colonnes, de filtrer le contenu de la table (par région) et de visualiser finalement le tableau de données dans R.

sheet <- "COM_2014"
data <- data.frame(read_excel("data/INSEE/pop-act2554-csp-cd-6814.xlsx", skip = 15, sheet = sheet)) 
colnames(data)

##  [1] "RR"                            "DR"                           
##  [3] "CR"                            "STABLE"                       
##  [5] "DR16"                          "LIBELLE"                      
##  [7] "csx_rec1taxtypac_rec1rpop2014" "csx_rec1taxtypac_rec2rpop2014"
##  [9] "csx_rec2taxtypac_rec1rpop2014" "csx_rec2taxtypac_rec2rpop2014"
## [11] "csx_rec3taxtypac_rec1rpop2014" "csx_rec3taxtypac_rec2rpop2014"
## [13] "csx_rec4taxtypac_rec1rpop2014" "csx_rec4taxtypac_rec2rpop2014"
## [15] "csx_rec5taxtypac_rec1rpop2014" "csx_rec5taxtypac_rec2rpop2014"
## [17] "csx_rec6taxtypac_rec1rpop2014" "csx_rec6taxtypac_rec2rpop2014"

data <- data[data$RR == "76",] # Ne prendre que les communes d'Occitanie (code 76)
View(data) # vOIR Si le fichier a correctement été importé. 

Pour notre analyse nous avons besoin de regrouper plusieurs colonnes. Créez de nouveaux champs en aggrégeant emploi + chômage par CSP. Renommez ces champs de façon explicite. Cela facilitera leur traitement par la suite. Agrégez également le code départemental au code communal afin de créer un identifiant (id) compatible avec le fond de carte IGN.

data$id <- paste0(data$DR, data$CR) # Creation du code commune
data[,"agr"]<- round(data[7] + data[8],0) # Actifs + chomeurs par CSP
data[,"art"] <- round(data[9] + data[10],0)
data[,"sup"] <- round(data[11] + data[12],0)
data[,"int"] <- round(data[13] + data[14],0)
data[,"emp"] <- round(data[15] + data[16],0)
data[,"ouv"] <- round(data[17] + data[18],0)
data[,"ouv_chom"] <- round(data[18],0)
data$total <- data$agr + data$art + data$sup + data$int + data$emp + data$ouv

Importez et visualisez le fond de carte (extrait IGN Geofla communes).

communes <- st_read(dsn = "data/IGN/occitanie.shp", stringsAsFactors = F)

## Reading layer `occitanie' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/IGN/occitanie.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4516 features and 17 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 428145.7 ymin: 6137116 xmax: 848019.7 ymax: 6439628
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(communes),col="white", border="black", lwd=0.1)

La Jointure est une étape cruciale qui permet de lier un fond de carte avec des données par un même identifiant. Il faut que les identifiants entre ces deux ressources soient exactement les mêmes.

Réalisez la jointure entre données et fond de carte avec l’instruction merge. Cela présuppose en amont que les identifiants géographiques entre ces deux ressources sont les mêmes (ça tombe bien, c’est le cas !). Identifiez bien le nom des champs sur lesquel doit porter la jointure.

communes <- merge(x = communes, y = data, by.x = "INSEE_COM", by.y = "id", all.x=T)

Pour finir, nettoyez le tableau. On ne garde que les colonnes utiles pour la suite…

colnames(communes)

##  [1] "INSEE_COM"                     "ID_GEOFLA"                    
##  [3] "CODE_COM"                      "NOM_COM"                      
##  [5] "STATUT"                        "X_CHF_LIEU"                   
##  [7] "Y_CHF_LIEU"                    "X_CENTROID"                   
##  [9] "Y_CENTROID"                    "Z_MOYEN"                      
## [11] "SUPERFICIE"                    "POPULATION"                   
## [13] "CODE_ARR"                      "CODE_DEPT"                    
## [15] "NOM_DEPT"                      "CODE_REG"                     
## [17] "NOM_REG"                       "RR"                           
## [19] "DR"                            "CR"                           
## [21] "STABLE"                        "DR16"                         
## [23] "LIBELLE"                       "csx_rec1taxtypac_rec1rpop2014"
## [25] "csx_rec1taxtypac_rec2rpop2014" "csx_rec2taxtypac_rec1rpop2014"
## [27] "csx_rec2taxtypac_rec2rpop2014" "csx_rec3taxtypac_rec1rpop2014"
## [29] "csx_rec3taxtypac_rec2rpop2014" "csx_rec4taxtypac_rec1rpop2014"
## [31] "csx_rec4taxtypac_rec2rpop2014" "csx_rec5taxtypac_rec1rpop2014"
## [33] "csx_rec5taxtypac_rec2rpop2014" "csx_rec6taxtypac_rec1rpop2014"
## [35] "csx_rec6taxtypac_rec2rpop2014" "agr"                          
## [37] "art"                           "sup"                          
## [39] "int"                           "emp"                          
## [41] "ouv"                           "ouv_chom"                     
## [43] "total"                         "geometry"

fields <- c("INSEE_COM", "LIBELLE", "agr", "art", "sup", "int", "emp", "ouv","ouv_chom", "total",  "geometry")
communes <- communes[,fields]
colnames(communes) <-  c("id", "name", "agr", "art", "sup", "int", "emp", "ouv","ouv_chom", "total",  "geometry")
head(communes)

## Simple feature collection with 6 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 556068.9 ymin: 6180702 xmax: 612006.7 ymax: 6219884
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##      id          name agr art sup int emp ouv ouv_chom total
## 1 09001 Aigues-Juntes   5   5   5   9   5   5        0    34
## 2 09002  Aigues-Vives   0  19  24  58  73  63       15   237
## 3 09003     Aiguillon   0   0   5  19  44  44       10   112
## 4 09004        Albiès   5  11  11  11  11  16       11    65
## 5 09005          Aleu   4   0   0   4   8   4        4    20
## 6 09006        Alliat   0   0   0  10   5   5        0    20
##                         geometry
## 1 MULTIPOLYGON (((572559.1 62...
## 2 MULTIPOLYGON (((607979.2 62...
## 3 MULTIPOLYGON (((612006.7 62...
## 4 MULTIPOLYGON (((594651 6188...
## 5 MULTIPOLYGON (((557843.5 61...
## 6 MULTIPOLYGON (((584509.3 61...

Le package cartography

library("cartography")

?cartography

Visualisation : carte en figurés proportionnels

Pour les caractères quantitatifs absolus la variable visuelle adaptée est la taille afin de préserver graphiquement les rapports de proportionnalité existant entre les différentes valeurs du tableau de données. Il n’est pas nécessaire de discrétiser les valeurs (cf exemple sur l’Europe).

# Ajuster les marges
par(mar = c(0.5,0.5,1.5,0.5)) 

# Cartographie
plot(st_geometry(communes), col="lightblue",border="white", lwd=0.2) # Fond de carte de la région Occitanie
propSymbolsLayer(x = communes, var = "total",  # Carte en figurés proportionnels 
                 symbols = "circle", col =  "red",
                 legend.pos = "right", border = "grey",
                 legend.title.txt = "Les actifs",
                 legend.style = "c")

top <- sort(data.frame(communes)[,"total"], decreasing = TRUE) 
labelLayer(x = communes[data.frame(communes)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) # Nommer les 15 premières villes d'Occitanie en nombre d'actifs
title("Région Occitanie")

Visualisation : carte combinant figurés proportionnels et typologie

Il est possible d’associer à ces figurés proportionnels la représentation d’un caractère qualitatif de différenciation.

La variable visuelle adaptée pour ce type de données est la couleur différentielle.

Exécutez le code ci-dessous, qui propose le calcul d’une typologie sur la surreprésentation d’ouvriers ou de cadres. Cette typologie peut être représentée simplement avec la fonction typoLayer.

Pour le choix des couleurs, vous pouvez utiliser color picker pour obtenir les couleurs que vous désirez, au format hexadécimal présenté ci-dessous.

# Typologie (plus d'ouvriers que de cadres ?)
communes$typo <- "Indéterminé"
communes$r <- communes$ouv / communes$sup
communes$r[is.na(communes$r)] <- 0
communes[communes$r > 1.1,"typo"] <- "Plus d'ouvriers que de cadres"
communes[communes$r < 0.91,"typo"] <- "Plus de cadres que d'ouvriers"

# Choix des couleurs
colouv <- "#dd4e44" # tonalité rouge
colsup <- "#63b269" # tonalité verte

# Ajuster les marges
par(mar = c(0.5,0.5,1.5,0.5)) 

# Cartographie
typoLayer(communes, var="typo", 
          legend.values.order = c("Plus d'ouvriers que de cadres", 
                                  "Plus de cadres que d'ouvriers",
                                  "Indéterminé" ),
          col = c(colouv, colsup, 'grey'), 
          border = NA, legend.title.txt = "Type dominant") 

On peut aussi décider d’associer cette représentation à la carte précédente en figurés proportionnels, grâce à la fonction propSymbolsTypoLayer :

# Ajuster les marges
par(mar = c(0.5,0.5,1.5,0.5)) 

# Combiner avec la carte précédente
plot(st_geometry(communes), col="lightblue",border=NA)
propSymbolsTypoLayer(x = communes, var = "total", var2 = "typo",
                     symbols = "circle",
                     inches = 0.4,
                     col = c(colouv, colsup, 'grey'),
                     border = "white", lwd = 0.1,
                     legend.var2.values.order = c("Plus d'ouvriers que de cadres", 
                                             "Plus de cadres que d'ouvriers",
                                             "Indéterminé" ),
                     legend.var.pos = "right", 
                     legend.var.title.txt = "Nombre d'actifs",
                     legend.var2.title.txt = "Type dominant")

labelLayer(x = communes[data.frame(communes)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) 
title("Cadres et ouvriers en région Occitanie")

Une astuce cartographique : le “cartogramme” de Dorling

A cette échelle géographique on remarque que les cercles se superposent les uns sur les autres. Cela ne nuit pas à la compréhension globale de l’organisation spatiale du phénomène, mais ne permet pas d’en avoir une lecture précise (à cette échelle de visualisation et pour cet espace d’étude). On peut aussi décider d’utiliser un package et une fonction qui permet d’éviter cela.

Pour ce faire, lancez la librairie cartogram et la fonction cartogram_dorling. Le cartogramme de Dorling replace la localisation des cercles en conservant leurs positions relatives et permet d’éviter leurs chevauchements.

library("cartogram")

# Ajuster les marges
par(mar = c(0.5,0.5,1.5,0.5)) 

# cartogram de Dorling
actifs <- cartogram_dorling(communes, "total", k = 1, m_weight = 1, itermax = 20)

plot(st_geometry(communes), col="lightblue",border=NA)

typoLayer(actifs, var="typo", 
          legend.values.order = c("Plus d'ouvriers que de cadres", 
                                  "Plus de cadres que d'ouvriers",
                                  "Indéterminé" ),
          col = c(colouv, colsup, 'grey'), 
          border = "white", lwd=0.2, legend.title.txt = "Type dominant",
          add = T) 
labelLayer(x = actifs[data.frame(actifs)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) 
title("Ouvriers et cadres en Occitanie, 2014")

A vous de jouer !

On souhaite réaliser une carte avec une autre catégorie socio-professionnelle : les agriculteurs (ou une autre CSP de votre choix). Pourquoi ne pas choisir des carrés à la place des cercles ?

Voici le code à compléter pour réaliser la carte :

plot(st_geometry(communes), col="",border="", lwd="")

propSymbolsLayer(...)

Objectifs

• Thématique : Réaliser des cartes de densité de population et de PIB par habitant pour les régions européennes.
• Espace d’étude : Europe, NUTS3
• Visualisation associée : cartes choroplèthes et cartes en grilles
• packages utilisés : sf, cartography, eurostat, reshape2

Téléchargement des données

Téléchargez les données ici : Télécharger

Décompressez le dossier et cliquez sur proj.Rproj. Le fichier exo3.R s’ouvre dans R-Studio.

Préparation des données

Chargez les géométries qui seront utilisées pour visualiser vos résultats. Il s’agit des géométries de référence des territoires européens adaptées à la cartographie thématique et proposées par l’UMS RIATE.

# Charger la librairie sf
library(sf)

# Fond de carte principal
regions <- st_read(dsn = "data/GREAT/nuts2013-level3.shp", stringsAsFactors = F) 

## Reading layer `nuts2013-level3' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/GREAT/nuts2013-level3.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1538 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -28.82873 ymin: 27.64158 xmax: 44.78307 ymax: 71.141
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

# Couche d'habillage
countries <- st_read("data/GREAT/world.shp",stringsAsFactors = F)

## Reading layer `world' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/GREAT/world.shp' using driver `ESRI Shapefile'
## Simple feature collection with 199 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -54.92918 xmax: 180 ymax: 83.57445
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Astuce ! Le package eurostat permet d’accéder à la base de données de référence européenne directement dans l’environnement de R. Cela permet d’avoir des données à jour et d’éviter des potentielles erreurs de manipulation de données.

Le code ci-dessous (à ne pas exécuter) est un exemple de procédure d’extraction de données (population, surface, PIB).

library("eurostat")
library("reshape2") 

# POPULATION
var <- "demo_r_pjangrp3" # Table Eurostat d'intérêt
data <- get_eurostat(var, time_format = "num") # Telecharger la table ESTAT
data <- subset(data, data$sex == "T") # Filtre des dimensions du tibble
data <- subset(data, data$age == "TOTAL") # Filtre des dimensions du tibble
data <- dcast(data, geo ~ time, value.var = "values") # Transformation en dataframe des dimensions restantes
colnames(data) <- c("geo","POP_2014","POP_2015","POP_2016","POP_2017") # Renommer les colonnes du dataframe de façon explicite

# SURFACE
var <- "reg_area3"
area <- get_eurostat(var, time_format = "num") # Telecharger la table ESTAT
area <- subset(area, area$landuse == "TOTAL") # Filtre des dimensions de la table
area <- dcast(area, geo ~ time, value.var = "values")
colnames(area) <- c("geo","AREA")

# PIB
var <- "nama_10r_3gdp"
gdp <- get_eurostat(var, time_format = "num")
gdp <- subset(gdp, gdp$unit == "MIO_EUR")
gdp <- dcast(gdp, geo ~ time, value.var = "values")
fields <- c("geo", "2014", "2015", "2016") # On ne garde que les valeurs de PIB correspondant aux valeurs de population disoponibles
gdp <- gdp[,fields]
colnames(gdp) <- c("geo","GDP_2014","GDP_2015","GDP_2016")

# Jointure avec les géométries
regions <- merge(x = regions, y = data, by.x = "id", by.y = "geo", all.x=T) # population
regions <- merge(x = regions, y = area, by.x = "id", by.y = "geo", all.x=T) # surface
regions <- merge(x = regions, y = gdp, by.x = "id", by.y = "geo", all.x=T) # surface

# Export du fichier contenant géométries et données
write_sf(obj = regions, dsn = "data/regions_data.shp")

Importez le jeu de données créé suite à cette étape de collecte de données, qui est maintenant associé aux géométries européennes de réference.

regions <- st_read("data/regions_data.shp", stringsAsFactors = F)

## Reading layer `regions_data' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/regions_data.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1538 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -28.82873 ymin: 27.64158 xmax: 44.78307 ymax: 71.141
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

par(mar = c(0.5,0.5,1.5,0.5)) 
plot(st_geometry(regions))

head(regions)

## Simple feature collection with 6 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 19.37065 ymin: 40.73386 xmax: 20.62819 ymax: 42.71386
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##      id    name POP_2014 POP_2015 POP_2016 POP_2017 AREA GDP_2014 GDP_2015
## 1 AL011   Diber   134332   131054   129056   125579   NA      315       NA
## 2 AL012  Durres   276058   277989   280205   284823   NA      923       NA
## 3 AL013   Kukes    84306    82428    81294    79559   NA      194       NA
## 4 AL014   Lezha   133748   131829   130258   129019   NA      334       NA
## 5 AL015 Shkodra   216477   213148   210168   207924   NA      526       NA
## 6 AL021 Elbasan   294938   290666   287606   283822   NA      739       NA
##   GDP_2016                       geometry
## 1       NA MULTIPOLYGON (((20.31134 41...
## 2       NA MULTIPOLYGON (((19.92766 41...
## 3       NA MULTIPOLYGON (((20.42293 42...
## 4       NA MULTIPOLYGON (((20.23201 41...
## 5       NA MULTIPOLYGON (((20.01829 42...
## 6       NA MULTIPOLYGON (((20.33483 40...

Pour choisir l’année de référence sur laquelle portera la visualisation, regardez le nombre de valeurs manquantes…

apply(regions,2,function(x) sum(is.na(x)))

##       id     name POP_2014 POP_2015 POP_2016 POP_2017     AREA GDP_2014 
##        0        0       51       51       51       51       63      162 
## GDP_2015 GDP_2016 geometry 
##      176     1219        0

… Et calculez le ratio sur lequel va porter la visualisation.

regions$DENS_2017 <- regions$POP_2017 / regions$AREA

Tout est prêt pour la visualisation !

Visualisation : cartes choroplèthes

La variable visuelle adaptée pour représenter des données quantitatives relatives est la couleur ordonnée (ou niveaux de gris). Le traitement de ce type de données implique au prélable de définir le nombre de classes, une méthode de discrétisation et une palette de couleurs sur laquelle portera la représentation.

La discrétisation consiste à découper une série statistique en classes de valeurs. Il n’existe pas de méthode optimum. Chaque méthode donnera une carte différente, plus ou moins conforme à la distribution de départ. L’agrégation de données en classes introduit de fait une erreur ou une distorsion dans la perception de cette distribution. Le choix d’une méthode de discrétisation dépend donc des propriétés de la distribution mais aussi des objectifs cartographiques que l’on s’est fixés.

var <- regions$DENS_2017

hist(var, probability = TRUE, nclass = 30)
rug(var)
moy <- mean(var)
med <- median(var)
abline(v = moy, na.rm = T, col = "red", lwd = 3)

## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...):
## "na.rm" is not a graphical parameter

abline(v = med, na.rm = T, col = "blue", lwd = 3)

## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...):
## "na.rm" is not a graphical parameter

summary(regions$DENS_2017)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
##     1.195    64.788   132.206   549.643   351.530 20779.676        63

# Choix de la palette
display.carto.all(10)

cols <- carto.pal(pal1 = "red.pal", n1 = 10)

library(cartography)

par(mfrow=c(1,1))

# Methode discretisation (quantiles)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "quantile")
hist(regions$DENS_2017, probability = TRUE, xlim = c(0,3000), breaks = breaks, col = cols,
     main = "Quantiles")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)

countries <- st_transform(countries,3035)
regions <- st_transform(regions,3035)

# Mise en page
par(mar = c(0.5,0.5,1.5,0.5)) # Ajuster les marges

plot(st_geometry(regions), col = NA, border = NA)
plot(st_geometry(countries), col = "#c6c4c4", border = NA, add =TRUE)

# Cartographie (ratio)
choroLayer(x = regions, 
           var = "DENS_2017",
           method ="quantile",
           nclass = 10,
           col = cols,
           border = NA,
           add = TRUE,
           legend.pos = "left",
           legend.title.txt = "Densité de population \n 2017 (hab/km²)",
           legend.values.rnd = 0)

plot(st_geometry(countries), col = NA, border = "#ffffff", lwd = 1, add =TRUE)

# Sources + elements d'habillage
layoutLayer(title = "Densité de population - Méthode des quantiles",
            author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
            scale = 200, tabtitle = FALSE,
            frame = FALSE,theme = "red.pal",
            south = TRUE)

A vous de jouer !

Créez une carte choroplèthe avec les paramètres suivants : * Indicateur : PIB / habitant 2014 * Discrétisation : quantiles (8 classes) * Palette de couleurs : Double palette (green.pal/red.pal)

Avant de passer à la cartographie, pensez à créer la variable PIB par habitant.

par(mfrow=c(1,1)) # On revient à un affichage pleine page
choroLayer(...)

Visualisation : cartes carroyées

Sur les cartes précédentes on a pu constater que la taille des entités territoriales est très différente d’un pays à l’autre. Pour certains pays (Allemagne, Pologne) les unités territoriales sont très petites et permettent de distinguer villes et périphéries. Dans d’autres cas, comme en France, celles-sont agrégées et rendent les résutats difficilement comparables.

Ces difficultés liées aux hétérogénéités de maillage renvoient à une série de travaux connus sous l’acronyme MAUP (Modifiable Area Unit Problem)(voir), souvent illustré de la façon suivante :

Plusieurs méthodes existent pour apporter des solutions à ces problèmes en cartographie.

La méthode du carroyage, par exemple, consiste à découper l’espace géographique en un maillage formé de carrés réguliers dans une projection donnée. La donnée est répartie sur ce quadrillage régulier au prorata de la surface représentée. Le quadrillage permet ainsi de s’affranchir des mailles administratives.

Comment faire ? La fonction getGridLayer du package cartography permet de construire ces grilles régulières. Expolorez la documentation pour vous sensibiliser aux arguments de la fonction.

La procédure est la suivante :

1. Définir le pas de grille (en mètres), son type et les variables qui seront transformées.
2. Pour les calculs de densité, convertir la surface en km².
3. Pour les autres indicateurs, combiner les données de la façon suivante : ratio = stock1/stock2.
4. Choisir sa palette, son nombre de classes et sa méthode de discrétisation.
5. Représenter la grille avec la fonction choroLayer de cartography.

# Transformation des donnees dans une grille reguliere de 50 km
mygrid50k <- getGridLayer(x = regions, cellsize = 50000 * 50000, 
                       type = "regular", var = c("POP_2017","POP_2014","GDP_2014"))

# Transformation des donnees dans une grille reguliere de 100 km
mygrid100k <- getGridLayer(x = regions, cellsize = 100000 * 100000, 
                          type = "regular", var = c("POP_2017","POP_2014","GDP_2014"))

# Conversion en km²
mygrid50k$densitykm <- mygrid50k$POP_2017 * 1000 * 1000 / mygrid50k$gridarea 
mygrid100k$densitykm <- mygrid100k$POP_2017 * 1000 * 1000 / mygrid100k$gridarea 


# Carte 1
par(mfrow=c(1,1))
par(mar = c(0.5,0.5,1.5,0.5)) # Ajuster les marges

cols <- carto.pal(pal1 = "blue.pal", pal2 = "red.pal",  n1 = 5, n2 = 5)  

plot(st_geometry(regions), col = NA, border = NA)
plot(st_geometry(countries), col = "#c6c4c4", border = NA, add =TRUE)

choroLayer(x = mygrid50k, var = "densitykm", 
           border = "grey80",col = cols, method ="quantile",
           nclass = 10, legend.pos = "topleft", legend.values.rnd = 1,
           lwd=0.5, add = T,
           legend.title.txt = "Densité de population, 2017")

layoutLayer(title = "Densité de population dans une grille régulière de 50km",
            author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
            scale = 200,
            frame = TRUE,
            col = "red",
            coltitle = "white",
            south = TRUE)

# Carte 2

plot(st_geometry(regions), col = NA, border = NA)
plot(st_geometry(countries), col = "#c6c4c4", border = NA, add =TRUE)

choroLayer(x = mygrid100k, var = "densitykm", 
           border = "grey80",col = cols, method ="quantile",
           nclass = 10, legend.pos = "topleft", legend.values.rnd = 1,
           lwd=0.5, add = T,
           legend.title.txt = "Densité de population, 2017")

layoutLayer(title = "Densité de population dans une grille régulière de 100km",
            author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
            scale = 200,
            frame = TRUE,
            col = "red",
            coltitle = "white",
            south = TRUE)

Aller plus loin…

Remarque : D’autres méthodes existent pour résoudre les problèmes du MAUP, comme les cartes lissées (ou de potentiel) qui consistent à s’affranchir complètement du maillage administratif sous-jacent et qui se basent sur la prise en compte du voisinage géographique selon une fonction mathématique (généralement gaussienne).

par(mar = c(0.5,0.5,1.5,0.5)) # Ajuster les marges
regions$GDPM_2014 <- regions$GDP_2014 * 1000000
ratio <- regions$GDPM_2014 / regions$POP_2014
summary(ratio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1907   17047   25646   26921   32685  428302     162

bks <- c(500, 5000, 10000, 15000, 20000, 26921, 30000, 35000, 40000, 428302)
cols <- c(rev(carto.pal("green.pal", 5)), carto.pal("orange.pal", 4))

plot(st_geometry(regions), border = NA, col = NA)
plot(st_geometry(countries), col = "#E3DEBF", border = NA, add = TRUE)

r <- regions[substr(regions$id,1,2) != "TR",]
r <- r[substr(r$id,1,2) != "IS",]
r <- r[substr(r$id,1,2) != "CH",]

smoothLayer(x = regions,
            var = 'GDPM_2014', var2 = 'POP_2014',
            span = 75000, beta = 2,
            resolution = 50000,
            legend.title.txt = "Potentiel de PIB par habitant, 2014\n(span = 75km)",
            col = cols, breaks = bks,
            lwd=0.2,
            mask = st_buffer(r,1),
            legend.pos = "topleft", legend.values.rnd = -2,
            add = T)

layoutLayer(title = "PIB par habitant, 2014", author = "Les serial mappers, 2018", 
            sources = "Source: Eurostat, 2018", frame = TRUE, scale = 500, north = FALSE, 
            theme = "sand.pal")

Objectifs

• Thématique : Les inégalités d’IDH et/ou d’espérance de vie.
• Espace d’étude : Le Monde, pays.
• Visualisation associée : carte chorplèthe, discontinuités.
• packages utilisés : sf, cartography, rnaturalearth, readxl, countrycode.

Téléchargement des données

Téléchargez les données ici : Télécharger

Décompressez le dossier et cliquez sur proj.Rproj. Le fichier exo4.R s’ouvre dans R-Studio.

Préparation des données

library("sf")
countries <- st_read(dsn = "data/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F)  

url <- "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip"
download.file(url = url, destfile = "data/countries.zip")
unzip("data/countries.zip", exdir = "data", unzip = "internal")
file.remove("data/countries.zip")
countries <- st_read(dsn = "data/ne_110m_admin_0_countries.shp", stringsAsFactors = F)

library("rnaturalearth")
countries <- ne_countries(scale = 50, type = "countries", continent = NULL,
                          country = NULL, geounit = NULL, sovereignty = NULL,
                          returnclass = "sf")

par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(countries))

countries <- countries[,c("adm0_a3", "admin", "geometry")]
colnames(countries) <- c("id","name","geometry")

countries <- countries[countries$id != "ATA",]

countries <- countries[!is.na(countries$id),]

par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(countries))

library("readxl")
file <-"data/worldbank/API_SP.DYN.LE00.IN_DS2_en_excel_v2_10081535.xls"
sheet <- "Data"
lifexp <- data.frame(read_excel(file, skip = 3, sheet = sheet))
lifexp <- lifexp[,c("Country.Code","Country.Name","X2016")]
colnames(lifexp) <- c("id","name","lifexp2016")
lifexp$lifexp2016 <- as.numeric(as.character(lifexp$lifexp2016))

file <-"data/hdr/Human Development Index (HDI).csv"
hdi <- read.csv2(file = file, sep = ",",skip = 1)
hdi <- hdi[,c("Country","X2016")]

head(hdi,4)

##        Country X2016
## 1  Afghanistan 0.494
## 2      Albania 0.782
## 3      Algeria 0.753
## 4      Andorra 0.856

library("countrycode")
hdi$id <- countrycode(hdi$Country, 'country.name', 'iso3c')

## Warning in countrycode(hdi$Country, "country.name", "iso3c"): Some values were not matched unambiguously:  Eswatini (Kingdom of)

hdi <- hdi[,c("id","Country","X2016")]
colnames(hdi) <- c("id","name","hdi2016")
hdi$hdi2016 <- as.numeric(as.character(hdi$hdi2016))  
head(hdi,4)

##    id         name hdi2016
## 1 AFG  Afghanistan   0.494
## 2 ALB      Albania   0.782
## 3 DZA      Algeria   0.753
## 4 AND      Andorra   0.856

Visualisation : carte choroplèthe

Chargement du package

library("cartography")

Changer la projection du fond de carte (voir le site spatialreference

countries <- st_transform(countries, "+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")
plot(st_geometry(countries))

Jointures des deux fichiers de données.

countries <- merge(countries, lifexp, by="id", all.x=TRUE)
countries <- merge(countries, hdi, by="id", all.x=TRUE)
countries <- countries[,c("id","name.x","lifexp2016","hdi2016","geometry")]
colnames(countries)[2] <- "name"

On travaille sur l’indice de développement humain

Analyse de la distribution.

var <- countries$hdi2016
hist(countries$hdi2016, probability = TRUE, nclass = 30)

On choisit une discretisation par la méthode des quantiles sans classe centrale.

breaks <- getBreaks(v = var, nclass = 6, method = "quantile")
cols <- carto.pal(pal1="blue.pal", n1 = 3, pal2 = "orange.pal", n2 = 3, middle = FALSE, transparency = FALSE)

Réalisation de la carte.

par(mar = c(0.5,0.5,1.5,0.5))
layoutLayer(extent = countries,
            bg="#EEEEEE",
            title = "Indice de développement humain en 2016",
            scale = NULL,
            sources = "Human development Report, 2018",
            author = "les serial mappers"
            )
choroLayer(x = countries, 
           var = "hdi2016",
           breaks = breaks,
           col = cols,
           border = "white",
           lwd=0.2,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "IDH, 2016",
           legend.values.rnd = 2)

Et si on ajoutait des lignes de discontinuités à cette carte?

L’instruction getBorders (du package cartography) permet d’extraire les frontières entre les unités spatiales.

borders <- getBorders(countries)
plot(st_geometry(borders), col="black", lwd=1)

Quid des discontinuités maritimes ? Par exemple entre l’Italie et la Tunisie ? Ou entre le Maroc et l’Espagne ? L’instruction getOuterBorders (package cartography) permet de les générer. Cette instruction peut prendre un peu de temps.

outer <- getOuterBorders(x = countries, res = 100000, width=500000)
plot(st_geometry(borders), col="#CCCCCC", lwd=1)
plot(st_geometry(outer), col="red", lwd=1, add=T)

On assemble les deux couches avec rbind.

b <- rbind(borders,outer)

Tout est prêt pour réaliser une carte de discontinuités.

# 1 - la carte choroplèthe
par(mar = c(0.5,0.5,1.5,0.5))
layoutLayer(extent = countries,
            bg="#EEEEEE",
            title = "Indice de développement humain en 2016",
            scale = NULL,
            sources = "Human development Report, 2018",
            author = "les serial mappers"
            )

choroLayer(x = countries, 
           var = "hdi2016",
           breaks = breaks,
           col = cols,
           border = "white",
           lwd=0.2,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "IDH, 2016",
           legend.values.rnd = 2)

# 2 - les discontinuités

discLayer(x = b, df = countries,
          var = "hdi2016", col="black", nclass=4,
          method="quantile", threshold = 0.25, sizemin = 0.2,
          sizemax = 8, type = "abs", 
          legend.title.txt = "Discontinuities\nabsolues",
          legend.pos = "bottomleft", add=TRUE)

A vous de jouer

Réaliser une carte de discontinuités sur l’espérance de vie

layoutLayer(...)
choroLayer(...)
discLayer(...

Objectifs

• Utiliser des données Open Street Map
• Généraliser un fond de carte

Téléchargement des données

Téléchargez les données ici : Télécharger

Décompressez le dossier et cliquez sur proj.Rproj. Le fichier bonus.R s’ouvre dans R-Studio.

Dé-Googliser la carto !

OpenStreetMap (OSM) est un projet de cartographie participative qui a pour but de constituer une base de données géographiques libre à l’échelle mondiale. OpenStreetMap vous permet de voir, modifier et utiliser des données géographiques dans le Monde entier. En résumé, c’est comme Google Maps, mais en mieux…

library("sf")
library("cartography")

depts <- st_read(dsn = "data/IGN/DEPARTEMENT.shp", stringsAsFactors = F)

## Reading layer `DEPARTEMENT' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/IGN/DEPARTEMENT.shp' using driver `ESRI Shapefile'
## Simple feature collection with 96 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 99217.1 ymin: 6049646 xmax: 1242417 ymax: 7110480
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

colnames(depts)

##  [1] "ID_GEOFLA"  "CODE_DEPT"  "NOM_DEPT"   "CODE_CHF"   "NOM_CHF"   
##  [6] "X_CHF_LIEU" "Y_CHF_LIEU" "X_CENTROID" "Y_CENTROID" "CODE_REG"  
## [11] "NOM_REG"    "geometry"

Paris <- depts[depts$CODE_DEPT == 75,]
par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(Paris))

Paname <- getTiles(x = Paris, type = "osm", crop = TRUE)
Paname <- getTiles(x = Paris, type ="opencycle", crop = TRUE)
tilesLayer(Paname)

plot(st_geometry(Paris), lwd=4, add=T)

library("osmdata")

## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright

prj <- st_crs(Paris)

bbox <- st_bbox(st_transform(Paris,4326))

q <- opq(bbox = bbox , timeout = 2000) %>% add_osm_feature(key = 'man_made', value = 'surveillance')
cameras <- osmdata_sf(q)$osm_points
cameras <- st_transform(cameras, prj)

cameras$ok <- st_intersects(st_geometry(cameras), st_geometry(Paris), sparse = FALSE)
cameras <- cameras[cameras$ok == T,]

Paname <- getTiles(x = Paris, type ="cartodark", crop = TRUE)
tilesLayer(Paname)
par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(Paris), lwd=1, border="white", lty=2, add=T)
plot(st_geometry(cameras), pch=20, col="red", add=T)

Généraliser un fond de carte

Toute carte est une reproduction réduite d’une portion de l’espace. Elle suppose donc une simplification des tracés, une réduction de la précision géométrique pour favoriser la lisibilité. C’est ce qu’on appelle la généralisation.

library("rmapshaper")

Import et affichage d’un fond de carte.

countries <- st_read(dsn = "data/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F)

## Reading layer `ne_110m_admin_0_countries' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/naturalearth/ne_110m_admin_0_countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 94 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(countries), col="#d1a2d3")

Généralisation avec sf (ça ne marche pas très bien… Pour l’instant ! ).

countries2 <- st_simplify(x = countries, preserveTopology = TRUE, dTolerance = 2)

## Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
## st_simplify does not correctly simplify longitude/latitude data, dTolerance
## needs to be in decimal degrees

par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(countries2), col="#d1a2d3")

Généralisation avec rmapshaper (Youpie !!!)

countries2 <- ms_simplify(countries, keep = 0.1)
par(mar = c(0.5,0.5,1.5,0.5))
plot(st_geometry(countries2), col="#d1a2d3")