Cartographie avec R

Logiciels pour la statistique spatiale
RTP CNRS “Réseau Interdisciplinaire autour de la Statistique” | Jeudi 23 Novembre 2017

Timothée Giraud
UMS RIATE
http://rgeomatic.hypotheses.org/

R

1

R est un langage et un environnement permettant de réaliser une variété de traitements statistiques et de représentations graphiques.

R est un logiciel libre sous license GNU General Public License.

R est multiplateforme (GNU/Linux, Windows, OS X…).

2

Un écosystème de packages

Les extensions (packages) sont mis à disposition sur le CRAN (Comprehensive R Archive Network).

11851 packages le 17 novembre 2017


http://blog.revolutionanalytics.com/2017/01/cran-10000.html

Les packages sont organisés en Task Views

La Task View “Spatial” : https://CRAN.R-project.org/view=Spatial

3

Un environnement de développement intégré (IDE).

RStudio

4

Et aussi

Des solutions de literate programming (Markdown, R Mardown)

Literate programming : Une explication de la logique du programme en langage naturel, entremêlée de morceaux de code source.

Des solutions de gestion de version (git, Github)

Pour diffuser et partager ses résultats, travailler en équipe ou faire appel à des contributeurs.

R pour la cartographie et l’analyse spatiale ?

Recherche reproductible

1

Accompagner les publications scientifiques des jeux de données et codes sources pour permettre aux collègues de reproduire les résultats.

Peng, 2011 Peng, 2011

2

Les cartes, comme les autres production graphiques ou statistiques sont des éléments à part entière des études scientifique.

mapspectrum

Chaînes de traitements

1

La grande majorité des cartes produites dans un contexte académique sont issues de processus complexes. Elles sont donc souvent produites en utilisant une grande variété de logiciels et de formats.

chain1

2

Cette variété de formats et de logiciels rend difficile la reproduction des cartes.

chain2

3

Simplifier les chaines de traitement pour couvrir les différentes étapes de la construction cartographique.

chain2

Les fonctionnalités spatiales de R

  • Import / Export
  • Manipulation et affichage
  • Géotraitements
  • Statistiques spatiales

Les packages “historiques”

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

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

  • 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
  • Dissolve
  • Buffer
  • Overlap / intersect / difference
  • Contains / within
  • Union

L’actualité et le futur : le package sf

  • Première release : 31 octobre 2016
  • Auteur principal et maintainer : Edzer Pebesma
  • Financement :

Import / Export

library(sf)
nuts3 <- st_read(dsn = "data/nuts3.shp")
Reading layer `nuts3' from data source `/home/tg/Documents/prz/Intro_a_la_carto_avec_R/data/nuts3.shp' using driver `ESRI Shapefile'
Simple feature collection with 1448 features and 7 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 2641758 ymin: 1427614 xmax: 7313157 ymax: 5411284
epsg (SRID):    NA
proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

Manipulation et affichage

plot(st_geometry(nuts3))

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

Géotraitements

Agrégation des polygones / dissolve

europe <- st_union(x = nuts3)
plot(st_geometry(nuts3), lwd = 0.5)
plot(europe, lwd = 2, border = "red", add=T)

Création de zones tampons / buffer

europeBuffer <- st_buffer(x = europe, dist = 50000)
plot(st_geometry(europe), col = "#92C5D6")
plot(europeBuffer, add = T, border = "red")

Ressources

La page GitHub du package sf

C’est dans ce dépôt GitHub que se déroule le développement du package et que se tiennent les discussions à son sujet.

Les ressources publiées par Edzer Pebesma

Le créateur et maintainer de sf a publié un certain de nombre de vignettes et de billets autour du package.

Billets de blogs / tutoriels

Le package cartography

Create and integrate maps in your R workflow. This package allows various cartographic representations such as proportional symbols, chroropleth, typology, flows or discontinuities maps. It also offers several features enhancing the graphic presentation of maps like cartographic palettes, layout elements (scale, north arrow, title…), labels, legends or access to some cartographic APIs.

Installation

  • Version stable (CRAN)
install.packages("cartography")



  • Version de développement (Github)
devtools::install_github(repo = "riatelab/cartography")

Cartes choroplèthes

library(cartography)

# chargement de données
data(nuts2006)

# 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 v2.0.2", 
            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 =  "seagreen4",
                 legend.pos = "right", inches = 0.35,
                 legend.title.txt = "Total\npopulation (2008)",
                 legend.style = "c")

# Ajout de labels
dflab <- nuts0.df[order(nuts0.df$pop2008, decreasing = TRUE),][1:8,]
dflab$lab <- paste(dflab$id, "\n", round(dflab$pop2008/1000000,0), "M", sep ="")

# Label plot of the 8 most populated countries
labelLayer(spdf = nuts0.spdf, 
           df = dflab, 
           txt = "lab", 
           col = "#690409", 
           cex = 0.8, 
           font = 2) 

# Ajout des titres, légende, sources, etc.
layoutLayer(title = "Countries Population in Europe",
            theme = 'green.pal',
            frame = FALSE,
            author = "cartography v2.0.2",
            sources = "Eurostat, 2008",
            scale = NULL,
            south = TRUE)

Cartes en symboles proportionnels colorés

# Load data
data(nuts2006)

# Compute the compound annual growth rate
nuts2.df$cagr <- (((nuts2.df$pop2008 / nuts2.df$pop1999)^(1/9)) - 1) * 100

# Plot a layer with the extent of the EU28 countries with only a background color
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
# Plot non european space
plot(world.spdf, col  = "#E3DEBF", border = NA, add = TRUE)
# Plot Nuts2 regions
plot(nuts2.spdf, col = "grey60",border = "white", lwd = 0.4, add = TRUE)

# Set a custom color palette
cols <- carto.pal(pal1 = "blue.pal", n1 = 2, pal2 = "red.pal", n2 = 4)

# Plot symbols with choropleth coloration
propSymbolsChoroLayer(spdf = nuts2.spdf, df = nuts2.df, 
                      var = "pop2008", inches = 0.1, 
                      var2 = "cagr", col = cols, 
                      breaks = c(-2.43,-1,0,0.5,1,2,3.1), 
                      border = "grey50", lwd = 0.75, 
                      legend.var.pos = "topright", 
                      legend.var.values.rnd = -3,
                      legend.var.title.txt = "Total Population", 
                      legend.var.style = "e", 
                      legend.var2.pos = "right", 
                      legend.var2.title.txt = "Compound Annual\nGrowth Rate") 

# layout
layoutLayer(title = "Demographic trends, 1999-2008", coltitle = "black",
            sources = "Eurostat, 2011", scale = NULL,
            author = "cartography v2.0.2", frame ="", col = NA)

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.df)
i j fij
DE14 AT11 1
DE21 AT11 1
DE23 AT11 1
DE26 AT11 2
DE91 AT11 1
DEB3 AT11 1
# Creation d'une couche de liens
twincities.sf <- getLinkLayer(x = nuts2.spdf, df = twincities.df[,1:2])

# Affichage des liens créés
plot(st_geometry(twincities.sf), 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(x = twincities.sf, df = twincities.df,   
              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 v2.0.2", 
            sources = "Sources: Adam Ploszaj & Wikipedia, 2011",
            scale = NULL, south = TRUE, frame = TRUE, col = NA, 
            coltitle = "black")

Discontinuités

library(cartography)

# Load data
data(nuts2006)

# Get a layer of borders between regions of countries borders
nuts2.contig <- getBorders(nuts2.spdf)

plot(nuts2.spdf, col = "grey", border = NA)
plot(st_geometry(nuts2.contig), col = 1:5, add=T, lwd = 3)

# get gdp/inhabitants
nuts2.df$gdpcap <- nuts2.df$gdppps2008/nuts2.df$pop2008 * 1e+06

# Plot a background layers
plot(nuts2.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col = "#E3DEBF", border = NA, add = TRUE)

# Plot GDP per capita with a choropleth layer
choroLayer(spdf = nuts2.spdf, df = nuts2.df, var = "gdpcap", 
           border = "grey20", lwd = 0.2, 
           col = carto.pal(pal1 = "green.pal", n1 = 3, "sand.pal", 3), 
           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(x = nuts2.contig, df = nuts2.df, 
          var = "gdpcap", type = "rel", 
          method = "equal", nclass = 3, threshold = 0.4, 
          sizemin = 0.7, sizemax = 6, col = "red", 
          legend.values.rnd = 1, legend.pos = "topright", add = TRUE,
          legend.title.txt = 
            "Discontinuities in \nGDP per Capita\n(relative)")

# Layout
layoutLayer(title = "Wealth Disparities in Europe, 2008", 
            author = "cartography v2.0.2", 
            sources = "Source: Eurostat, 2011", frame = TRUE, 
            scale = 500, north = FALSE, theme = "grey.pal")

Carroyages

library(cartography)
# Load data
data(nuts2006)

nuts3.spdf@data <- nuts3.df

# Create a grid layer
mygrid <- getGridLayer(x = nuts3.spdf, 
                       cellsize = 100000 * 100000, 
                       var = c("pop2008", "pop1999"), 
                       type = "hexagonal")

# Compute the compound annual growth rate
mygrid$cagr <- (((mygrid$pop2008/mygrid$pop1999)^(1/9)) - 1) * 100
v <- getBreaks(v = mygrid$cagr, method = "quantile", nclass = 10)
v[5] <- 0

# set a color palette
cols <- c("#f18b61", "#f7b48c", "#f3d9b7", "#f1eccd", 
          "#c0dec2", "#91caa4", "#63b285", "#329966", 
          "#26734d", "#1a4c33")

# Plot a background layers
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
plot(world.spdf, col = "#CCCCCC", border = NA, add = TRUE)

# plot the choropleth grid
choroLayer(x = mygrid, var = "cagr", add = TRUE, 
           col = cols, lwd = 0.6, border = "#FFFFFF60", 
           legend.pos = "right", breaks = v, legend.values.rnd = 2, 
           legend.title.txt = "Compound Annual\nGrowth Rate")

# plot countries boundaries
plot(nuts0.spdf, add = T, col = NA, border = "#56514c", lwd = 0.7)

# Plot a layout
layoutLayer(title = "Demographic Trends, 1999-2008", 
            author = "cartography v2.0.2", 
            sources = "Source: Eurostat, 2011", frame = TRUE, 
            scale = 500, north = TRUE, theme = "taupe.pal")

Cartes lissées

# Load data
data(nuts2006)

nuts3.spdf@data = nuts3.df
# Create a grid layer
mygrid <- getGridLayer(x = sf::st_as_sf(nuts3.spdf), 
                       cellsize = 50000 * 50000, 
                       type = "regular", 
                       var = c("pop2008", "gdppps2008"))


# Compute data for the grid layer
mygrid$gdp <- mygrid$gdppps2008*1000000

# list of breaks
v <- c(2920, 5000, 10000, 15000, 20000, 23500, 30000, 35000, 40000, 42720)
# Plot a layer with the extent of the EU28 countries with only a background
# color
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")

# Plot non european space
plot(world.spdf, col = "#E3DEBF", border = NA, add = TRUE)

# set a color palette
cols <- c(rev(carto.pal("green.pal", 5)), carto.pal("orange.pal", 4))

# compute & display the potential map
smoothLayer(x = mygrid, var = "gdp", var2 = "pop2008", 
            breaks = v, span = 1e+05, beta = 2, mask = nuts0.spdf, 
            resolution = 49000, col = cols, 
            legend.title.txt = "Potential\nGDP per capita\n(in euros)", 
            legend.values.rnd = -2, 
            border = "grey80", lwd = 0.5, add = T, 
            legend.pos = "topright")

# plot Europe contour
plot(rgeos::gBuffer(nuts0.spdf, FALSE, 1), 
     add = T, col = NA, border = "grey50")

# plot a layout
layoutLayer(title = "Wealth Inequalities in Europe, 2008", 
            author = "cartography v2.0.2", 
            sources = "Source: Eurostat, 2011", 
            frame = TRUE, scale = 500, north = FALSE, 
            theme = "sand.pal")

# plot a text on the map
text(x = 6271272, y = 3743765, labels = "Distance function:\n- type = exponential\n- beta = 2\n- span = 100 km", 
     cex = 0.8, adj = 0, font = 3)

La présentation est accessible à cette adresse :
https://rcarto.github.io/Intro_a_la_carto_avec_R/index.html

Son code source est accessible sur GitHub.

cartography
Github: https://github.com/riatelab/cartography
CRAN: https://cran.r-project.org/web/packages/cartography/

rgeomatic
https://rgeomatic.hypotheses.org

sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.3 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rmdformats_0.3.4 knitr_1.17      

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.14    bookdown_0.5    digest_0.6.12   rprojroot_1.2  
 [5] mime_0.5        R6_2.2.2        xtable_1.8-2    backports_1.1.1
 [9] magrittr_1.5    evaluate_0.10.1 highr_0.6       stringi_1.1.6  
[13] rstudioapi_0.7  miniUI_0.1.1    rmarkdown_1.8   tools_3.4.2    
[17] stringr_1.2.0   questionr_0.6.2 shiny_1.0.5     httpuv_1.3.5   
[21] yaml_2.1.14     compiler_3.4.2  htmltools_0.3.6