Cartographie et traitement de l’information géographique avec R

R

R est un outil parmis d’autres dans la chaine de traitement géographique.
Il existe de nombreux logiciels correspondants à differents métiers, des chaînes de traitement parfois complexes.

































R permet de simplifer les chaînes de traitement en travaillant dans un environnement intégré basé sur un seul langage, un seul logiciel.

































Les fonctionnalités spatiales de R

Trois packages sont incontournables.

Le package rgdal, import/export d’objets spatiaux et gestion des projections cartographiques

rgdal est une interface entre R et les librairies GDAL (Geospatial Data Abstraction Library) et PROJ4.

Import du shapefile des régions européennes (NUTS3) :

library("rgdal")
nuts3 <- readOGR(dsn = "data", layer = "nuts3", verbose = TRUE)
OGR data source with driver: ESRI Shapefile 
Source: "data", layer: "nuts3"
with 1448 features
It has 7 fields

































Le package sp, manipulation et affichage d’objets spatiaux

sp fournit des classes et des methodes pour les données spatiales dans R.

Affichage des NUTS3 :

library("sp")
class(nuts3)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
nuts3@proj4string
CRS arguments:
 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
+ellps=GRS80 +units=m +no_defs 
nuts3@bbox
      min     max
x 2641758 7313157
y 1427614 5411284
head(nuts3@data)
     id birth_2008 death_2008 gdppps1999 gdppps2008   pop1999 pop2008
0 AT111        272        445        509        641  39148.91   37452
1 AT112       1195       1480       2262       3272 137469.46  146383
2 AT113        748       1142       1368       1783 100113.63   97350
3 AT121       2480       2119       4091       5807 238591.26  241083
4 AT122       2266       2577       4130       5774 244959.00  251838
5 AT123       1363       1392       3098       4504 144269.11  147468
plot(nuts3)

plot(nuts3, col = "#DAE3E6", border = "#8A0641", lwd = 0.5)

































Le package rgeos, géotraitements

rgeos donne accès à la librairie d’opérations spatiales GEOS (Geometry Engine - Open Source) qui permet notamment d’effectuer les géotraitements suivants :

  • Area / Perimeter
  • Distances
  • Buffer
  • Overlap / intersect / difference
  • Contains / within
  • Union
  • Dissolve

Agrégation des polygones / dissolve

library("rgeos")
europe <- gUnaryUnion(spgeom = nuts3)
plot(nuts3, lwd = 0.5)
plot(europe, lwd = 2, border = "red", add = T)

Création de zones tampons / buffer

library("rgeos")
europeBuffer <- gBuffer(spgeom = europe, width = 50000)
plot(europe, col = "#92C5D6")
plot(europeBuffer, add = T, border = "red")

































































Rappel de cartographie / sémiologie graphique

Les données géographiques sont de quatre types.
A chaque type s’applique un traitement statistique et graphique de l’information spécifique.

































Les données géographiques s’appliquent à 3 types d’objets géométriques que l’on peut classer selon leur nombre de dimension sur un plan : le point, la ligne, le polygone.

































Le type de donnée associé au type d’implantation graphique détéremine le(s) mode(s) de représentation cartographique possible(s).

« La graphique utilise les propriétés de l’image visuelle pour faire apparaitre les relations de différence, d’ordre et de proportionnalité entre les données » (J. Bertin)

































Le package cartography

Conçu comme une boite à outil dédiée à la cartographie thématique, le package cartography est développé au sein de l’UMS RIATE (CNRS, CGET, Université Paris Diderot) par Nicolas Lambert et Timothée Giraud.

































Installation

La version stable est sur le CRAN.

install.packages("cartography")

La version de développement est hébergée sur GitHub.

require(devtools)
devtools::install_github("Groupe-ElementR/cartography")

Vous pouvez nous faire remonter d’éventuels bugs ici.

































Les données

Les exemples utilisent les jeux de données fournis avec le package. Ces données portent sur les maillages régionaux européens NUTS.

# Chargement de la librairie
library(cartography)

# Import de données dans la session
data(nuts2006)

ls()
 [1] "coasts.spdf"    "countries.spdf" "frame.spdf"     "graticule.spdf"
 [5] "nuts0.df"       "nuts0.spdf"     "nuts1.df"       "nuts1.spdf"    
 [9] "nuts2.df"       "nuts2.spdf"     "nuts3.df"       "nuts3.spdf"    
[13] "twincities"     "world.spdf"    

Les objets terminant par “.spdf”" sont des objets spatiaux et ceux terminant par “.df” sont des data frames.

































Carte Choroplethes

# Calcul du taux de croissance annuel moyen
nuts2.df$cagr <- 100 * (((nuts2.df$pop2008/nuts2.df$pop1999)^(1/9)) - 1)

# Cartographie
choroLayer(spdf = nuts2.spdf, df = nuts2.df, var = "cagr")
title("Taux de croissance en Europe")

Après ce premier jet, il est ensuite possible de paramétrer très finement la carte : palette de couleurs, discrétisation, légende, couches d’habillage…

# Construire une palette de couleurs
cols <- carto.pal(pal1 = "green.pal", n1 = 2, 
                  pal2 = "red.pal", n2 = 4) 

# Affichage de couches d'habillage
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col  = "#E3DEBF", border=NA, add=TRUE)

# Cartographie du taux de croissance annuel moyen
choroLayer(spdf = nuts2.spdf, df = nuts2.df, var = "cagr",#
           breaks = c(-2.43,-1.0,0.0,0.5,1.0,2.0,3.1), 
           col = cols,
           border = "grey40",
           lwd = 0.5, 
           legend.pos = "right",
           legend.title.txt = "taux de croissance\nannuel moyen", 
           legend.values.rnd = 2, 
           add = TRUE) 

# Affichage de couches d'habillage
plot(nuts0.spdf,border = "grey20", lwd=0.75, add=TRUE)

# Ajout des titres, légende, sources, etc.
layoutLayer(title = "Taux de croissance en Europe", 
            author = "cartography", 
            sources = "Eurostat, 2008", frame = TRUE, col = NA, 
            scale = NULL,coltitle = "black",
            south = TRUE) 

































Cartes en symboles proportionnels

Cartographie d’un stock (la population nationale) avec des figurés proportionnels.

# Affichage de couches d'habillage
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col  = "#E3DEBF", border=NA, add=TRUE)
plot(nuts0.spdf, col = "#D1914D",border = "grey80", add=TRUE)

# Cartographie de la population des pays en cercles proportionnels
propSymbolsLayer(spdf = nuts0.spdf, df = nuts0.df,#
                 var = "pop2008", 
                 symbols = "circle", col =  "#920000",
                 legend.pos = "right",
                 legend.title.txt = "Total\npopulation (2008)",
                 legend.style = "c")

# Ajout des titres, légende, sources, etc.
layoutLayer(title = "Countries Population in Europe",
            author = "cartography",
            sources = "Eurostat, 2008",
            scale = NULL,
            south = TRUE)

































Cartes de flux

Il s’agit de représenter des données, agrégées à un niveau régional, sur les jumelages entre villes.

# Données sur les jumelages
head(twincities)
     i    j fij
1 DE14 AT11   1
2 DE21 AT11   1
3 DE23 AT11   1
4 DE26 AT11   2
5 DE91 AT11   1
6 DEB3 AT11   1
# Creation d'une couche de liens
twincities.spdf <- getLinkLayer(spdf = nuts2.spdf, df = twincities) 

# Affichage des liens créés
plot(twincities.spdf, lwd = 0.2)

# Affichage de couches d'habillage
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col  = "#E3DEBF", border=NA, add=TRUE)
plot(nuts2.spdf, col = "#D1914D",border = "grey80", add=TRUE)

# Cartographie des liens
gradLinkLayer(spdf = twincities.spdf, df = twincities,   #
              spdfids = "i", spdfide = "j", dfids = "i", dfide = "j", 
              var = "fij", 
              breaks = c(2,5,15,20,30), 
              lwd = c(0.1,1,4,10), 
              col = "#92000090",
              legend.pos = "right", legend.frame = TRUE,
              legend.title.txt = "Number of Agreements\n(regional level)",
              add = TRUE)

# Ajout des titres, légende, sources, etc.
layoutLayer(title = "International Twinning Agreements Between Cities", 
            author = "cartography", 
            sources = "Sources: Adam Ploszaj & Wikipedia, 2011",
            scale = NULL, south = TRUE, frame = TRUE, col = NA, 
            coltitle = "black")

































Discontinuités

Discontinuités de richesses entre Etats.

# Construction les polylignes des frontières inter-étatiques
nuts0.contig.spdf <- getBorders(nuts0.spdf)

plot(nuts0.spdf, bg = "#A6CAE0", col = "#D1914D",border = "grey80")
plot(nuts0.contig.spdf, col = 1:nrow(nuts0.contig.spdf), lwd = 2, add=TRUE)

head(nuts0.contig.spdf@data)
         id id1 id2
AT_CH AT_CH  AT  CH
AT_CZ AT_CZ  AT  CZ
AT_DE AT_DE  AT  DE
AT_HU AT_HU  AT  HU
AT_IT AT_IT  AT  IT
AT_LI AT_LI  AT  LI
# Calcul du PIB/habitant
nuts0.df$gdpcap <- nuts0.df$gdppps2008/nuts0.df$pop2008*1000000

# Affichage de couches d'habillage
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col  = "#E3DEBF", border=NA, add=TRUE)

# Cartographie du PIB/habitants
choroLayer(spdf = nuts0.spdf, df = nuts0.df, var = "gdpcap", #
           border = "grey80",
           col = carto.pal(pal1 = "kaki.pal", n1 = 6), 
           method = "quantile",
           nclass = 6, add=TRUE, 
           legend.pos = "right", 
           legend.values.rnd = -2,
           legend.title.txt = "GDP per Capita\n(in euros)")

# Plot discontinuities
discLayer(spdf = nuts0.contig.spdf, df = nuts0.df, var = "gdpcap", 
          type = "rel", 
          method = "equal", 
          nclass = 4, 
          threshold = 0.5, 
          sizemin = 0.5, 
          sizemax = 6, 
          col="red",
          legend.values.rnd = 1,
          legend.title.txt = "Discontinuities in \nGDP per Capita\n(relative)",
          legend.pos = "topright", 
          add=TRUE)

# Ajout des titres, légende, sources, etc.
layoutLayer(title = "Wealth Disparities in Europe",
            coltitle = "black", col = NA,
            sources = "Eurostat, 2011", 
            scale = NULL,
            author = "cartography", 
            frame = FALSE)

































Carroyages

Transformation de données dans un maillage hétérogène vers une grille régulière.

# Creation d'une grille régulière
mygrid <- getGridLayer(spdf=nuts2.spdf, cellsize = 200000)

# Affichage de la grille
plot(mygrid$spdf)

# Adaptation des données à la grille
datagrid.df <- getGridData(x = mygrid, df = nuts2.df, var = "pop2008") 
datagrid.df$densitykm <- datagrid.df$pop2008_density*1000*1000

# Affichage de couches d'habillage
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col  = "#E3DEBF", border=NA, add=TRUE)

# Cartographie de la densité de population 
choroLayer(spdf = mygrid$spdf,                         #
           df = datagrid.df,
           var = "densitykm", 
           border = "grey80", 
           col = carto.pal(pal1 = "wine.pal", n1 = 6), 
           legend.pos = "topright",
           method = "q6", add = TRUE,
           legend.title.txt = "Population Density\n(inhabitant/km²)")

# Ajout des titres, légende, sources, etc.
layoutLayer(title = "Population Density", 
            coltitle = "black", col = NA,
            sources = "Eurostat, 2011", scale = NULL, 
            author = "cartography", frame = FALSE)

































Fonds de carte OpenStreetMap

Import de tuiles issues d’OpenStreetMap

data("nuts2006")
# extract Denmark
spdf <- nuts0.spdf[nuts0.spdf$id == "DK", ]
# Download the tiles, extent = Denmark
den <- getTiles(spdf = spdf, type = "osmtransport", crop = TRUE, zoom = 7)
class(den)
[1] "RasterBrick"
attr(,"package")
[1] "raster"
# Plot the tiles
tilesLayer(den)
# Plot countries
plot(spdf, border = "red", lwd = 2, add = TRUE)
# Map tiles sources
mtext(text = "Map data © OpenStreetMap contributors, under CC BY SA.", side = 1, 
    adj = 0, cex = 0.7, font = 3, line = -2)

































Cas d’usage

Etude de la répartition des ouvriers en Ile-de-France

Import des données

# Import du fond de carte communal d'IDF (source : IGN)
library(rgdal)
library(sp)
comidf.spdf <- readOGR(dsn = "data/", layer = "idf", verbose = F)
plot(comidf.spdf)
head(comidf.spdf@data)
     ID
0 78233
1 91223
2 95445
3 95142
4 95127
5 95074
# Création d'une couche du contour de l'idf
library(rgeos)
idf.spdf <- gUnaryUnion(comidf.spdf)
plot(comidf.spdf)
plot(idf.spdf, lwd = 4, border = "red", add = TRUE)

# Import des données sur les communes (source : INSEE 2011)
comidf <- read.csv("data/idfdata.csv")
head(comidf)
     ID                    NAMES population agriculteurs artisans cadres
1 75101 Paris 1er Arrondissement      17443           19      763   5230
2 75102 Paris 2e  Arrondissement      22927            7      832   7138
3 75103 Paris 3e  Arrondissement      36120           24     1326  11154
4 75104 Paris 4e  Arrondissement      27887           21     1036   7770
5 75105 Paris 5e  Arrondissement      60800           28     1388  16527
6 75106 Paris 6e  Arrondissement      43880           61     1551  10687
  intermediaires employe ouvrier total
1           2134    1778     438 10362
2           3052    2248    1169 14446
3           5003    3418    1160 22085
4           3578    2696     638 15739
5           6857    4671    1282 30753
6           4065    3414     897 20675

Répartition des ouvriers (valeur absolue)

library(cartography)
plot(idf.spdf, col = "grey80")
propSymbolsLayer(spdf = comidf.spdf, df = comidf,     #
                 var = "ouvrier", 
                 col = "#920000", legend.pos = "topright",inches = 0.15, 
                 legend.title.txt = "Nombre d'ouvriers")
layoutLayer(title = "Répartition des ouvriers en Ile-de-France", 
            col = "#080E5B",
            sources = "Insee, 2011", author = "T. Giraud", 
            frame = T)

Part des ouvriers dans la population

# Calcul de la part des ouvriers dans la population active occupée
comidf$partouvrier <- comidf$ouvrier/comidf$total * 100

# discrétisation en quantiles de la distribution de la part des ouvriers
bks <- quantile(comidf$partouvrier, seq(0,1,length.out = 9))

# création d'une palette de couleurs
cols <- carto.pal(pal1 = "turquoise.pal", n1 = length(bks)-1)

# Affichage de la carte
choroLayer(spdf = comidf.spdf, df = comidf,             #
           col = cols, border = "grey50",
           legend.title.txt = "Part des ouvriers\n(en %)", 
           legend.pos = "topright",
           legend.title.cex = 0.7, 
           var = "partouvrier", breaks = bks)
layoutLayer(title = "Répartition des ouvriers en IDF",col = "#080E5B",
            sources = "Insee, 2011", author = "T. Giraud", 
            frame = T)

Combinaison du nombre et de la part des ouvriers

# Affichage de la carte
plot(idf.spdf, col = "grey80")
propSymbolsChoroLayer(spdf = comidf.spdf, df = comidf,                       #
                      var = "ouvrier", var2 = "partouvrier",
                      col = cols, border = "grey50", inches = 0.15,
                      breaks = bks, legend.var2.values.rnd = 0,
                      legend.var.title.txt = "Nombre\nd'ouvriers", 
                      legend.var2.title.txt = "Part des ouvriers",
                      legend.var.pos = "bottomright",
                      add=TRUE)

Part potentiel d’ouvrier dans un voisinnage

Des modèles d’interaction spatiale permettant de saisir l’influence exercée par un lieu sur tous les autres.

Le modèle de Stewart (accessibilité potentielle, potentiel gravitationnel ou accessibilité gravitationnelle) :

\[ A_i = \sum_{j=1}^n O_j f(d_{ij}) \]

  • \(A_i\) potentiel en \(i\)
  • \(O_j\) stock de population sur \(j\)
  • \(f(d_{ij})\) fonction négative de la distance entre \(i\) et \(j\)

Potentiel dans un maillage

library(SpatialPosition)
# Jointure fond de carte / données
comidf.spdf@data <- data.frame(comidf.spdf@data, 
                               comidf[match(comidf.spdf@data$ID,comidf$ID), -1])
row.names(comidf.spdf) <- as.character(comidf.spdf@data$ID)
head(comidf.spdf@data)
         ID             NAMES population agriculteurs artisans cadres
78233 78233      Feucherolles       2847           16       88    624
91223 91223           Etampes      24013           34      309   1222
95445 95445 Nerville-la-Foret        642            0       35     89
95142 95142             Chars       1862            4       42    147
95127 95127             Cergy      58341           12      706   5829
95074 95074         Boisemont        746            0       40    120
      intermediaires employe ouvrier total partouvrier
78233            264     144      72  1208    5.960265
91223           2768    3971    3046 11350   26.837004
95445             82      47      23   276    8.333333
95142            293     280     172   938   18.336887
95127           8521   10100    5254 30422   17.270396
95074            104      96      36   396    9.090909
# Potentiel d'ouvriers dans un voisinnage de 7.5 kilomètres
ouv <- stewart(knownpts = comidf.spdf, unknownpts = comidf.spdf,          #
               typefct = "exponential",
               varname = "ouvrier", span = 7500, beta = 3)
# Potentiel de population active dans un voisinnage de 7.5 kilomètres
tot <- stewart(knownpts = comidf.spdf, unknownpts = comidf.spdf,               #
               typefct = "exponential",
               varname = "total", span = 7500, beta = 3)
pot <- data.frame(id = tot$ID, partouvrierpot = ouv$OUTPUT/tot$OUTPUT * 100)

# cartographie
choroLayer(spdf = comidf.spdf, df = pot,              #
           border = NA,
           legend.title.txt = "Part des\nouvriers\n(en %)*", 
           legend.pos = "left",
           legend.title.cex = 0.7, 
           var = "partouvrierpot", 
           breaks = bks, add=F)
titre <- "Répartition des ouvriers en IDF - Voisinnage de 7.5 kilomètres"
layoutLayer(title = titre,
            col = "#080E5B",
            sources = "* Fonction d'interaction spatiale : exponentielle
  Facteur d'impédence (beta) : 3
  Portée : 7.5 kilomètres
  Sources : Insee, 2011", author = "T. Giraud", 
            frame = T, south = T)

Potentiel sur une surface continue

pot.spdf <- quickStewart(spdf = comidf.spdf, df = comidf.spdf@data,                #
                         var = "ouvrier", var2 = "total", breaks = bks/100,
                         span = 7500, beta = 3, mask = idf.spdf)
breaks <- sort(c(unique(pot.spdf$min), max(pot.spdf$max)), decreasing = FALSE)
cartography::choroLayer(spdf = pot.spdf, df = pot.spdf@data,              #
                        var = "center", breaks = breaks, border = NA, 
                        legend.pos = "left", legend.values.rnd = 2,
                        legend.title.txt = "Part des\nouvriers\n(en %)*")

  titre <- "Répartition des ouvriers en IDF - Voisinnage de 7.5 kilomètres"
layoutLayer(title = titre,
            col = "#080E5B",
            sources = "* Fonction d'interaction spatiale : exponentielle
  Facteur d'impédence (beta) : 3
  Portée : 7.5 kilomètres
  Sources : Insee, 2011", author = "T. Giraud", 
            frame = T, south = T)
plot(idf.spdf, add=T)

Discontinuités dans la répartition des ouvriers

# extraction de la zone Paris et petite couronne
comppc.spdf <- comidf.spdf[substr(comidf.spdf$ID,1,2) %in% c("75", "92", "93", "94"), ]

# extraction de frontières de polygones
ppccomcontig <- getBorders(spdf = comppc.spdf)

# téléchargement d'une tuile osm
osm <- getTiles(spdf = comppc.spdf, type = "mapquestosm", crop = TRUE, zoom = 11)
tilesLayer(osm)

# affichage de la part des ouvriers
choroLayer(spdf = comppc.spdf, df = comidf, var = "partouvrier",  #
           border = NA,
           col = paste0(carto.pal(pal1 = "wine.pal", n1 = 6),"CC"), 
           method = "quantile",
           nclass = 6,
           legend.pos = "topleft", legend.frame = TRUE, add=TRUE,
           legend.title.txt = "Part des ouvriers")

# affichage des discontinuités
discLayer(spdf = ppccomcontig, df = comidf, var = "partouvrier", 
          type = "rel", 
          method = "equal", 
          nclass = 4, 
          threshold = 0.4, 
          sizemin = 1, 
          sizemax = 6, 
          col="blue",
          legend.values.rnd = 1, legend.frame = TRUE,
          legend.title.txt = "Discontinuités dans\nla part des ouvriers",
          legend.pos = "topright", 
          add=TRUE)

Un peu d’interactivité?

library(leaflet)
# initialiser une carte
m <- leaflet()
# carte avec emprise mondiale
m <- addTiles(map = m)
# centrer sur un point
m <- setView(map = m, lng = 2.519182, lat = 48.896768, zoom = 12)
# ajout d'une couche de poly
m <- addPolygons(map = m,                                                  #
                 data = spTransform(comppc.spdf, "+init=epsg:4326"), 
                 opacity = 0.5,
                 color = "blue", stroke = TRUE,
                 weight = 2, popup = NULL,
                 options = list(clickable = FALSE),
                 fill = T, fillColor = "#B3C4B3",
                 fillOpacity = 0.2)

# ajout d'un pop-up
m <- addPopups(map = m, lng = 2.519182, lat = 48.896768, 
               popup = '<a href="https://fr.wikipedia.org/wiki/Le_Raincy">
               "Le Petit Neuilly du 93"</a>')


# affichage de la carte
m

































Quelques ressources en ligne

































Merci de votre attention!

Le package sur GitHub : https://github.com/Groupe-ElementR/cartography/

Carnet de recherche R Géomatique : http://rgeomatic.hypotheses.org/

































Timothée Giraud - CNRS
e-mail : mjsh GitHub : Blog :

R à l’usage des sciences sociales - 2016