Cartographie Reproductible avec R

T. Giraud & N. Lambert

9 mars 2018 - ENS, Lyon

“La science est infaillible,
mais les savants se trompent toujours.”
Anatole France, 1889
science

La recherche reproductible

  • Les sciences se reposent sur le principe de reproductibilité.

  • La reproductibilité est un élément permettant d’évaluer la validité des résultats.

  • En sciences sociales, l’utilisation de l’informatique facilite la reproductibilité.

Le spectre de la reproductibilité

Peng, 2011 Peng, 2011

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

Cartographie reproductible ?

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

Peng, 2011

Cartographie reproductible ?

  • Garder une trace cohérente des traitements effectués
  • Partager ses méthodes
  • Etre transparent, s’exposer à la critique
  • Faciliter les mises à jour
  • Travailler collectivement
  • Automatiser des taches
  • Lier fortement l’analyse et la représentation

Cartographie reproductible ?

Drawing
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.

Cartographie reproductible ?

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

Cartographie reproductible !

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

Les fonctionnalités spatiales de R
sf gif

Les packages historiques

Le package sf

  • Les fonctionnalités de sp, rgeos et rgdal dans un package unique.

  • Manipulation plus aisée, objets plus simples

  • Auteur principal et maintainer : Edzer Pebesma (auteur de sp)

  • Compatible avec les syntaxes pipe et les opérateurs du tidyverse.

Le package sf

Format des objets spatiaux sf

format sf

Le package sf

Import de données

library(sf)
mtq <- st_read("data/martinique.shp")
Reading layer `martinique' from data source `/home/tim/Documents/prz/cartographie-reproductible/data/martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID):    32620
proj4string:    +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs

Le package sf

Affichage de données

plot(st_geometry(mtq))

Le package sf

Extraire les centroïdes

mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)

Le package sf

Construire une matrice de distances

mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]
Units: m
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]     0.000 35297.56  3091.501 12131.617 17136.310
[2,] 35297.557     0.00 38332.602 25518.913 18605.249
[3,]  3091.501 38332.60     0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702     0.000  7177.011
[5,] 17136.310 18605.25 20226.198  7177.011     0.000

Le package sf

Agréger des polygones

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")

Le package sf

Construire une zone tampon

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")

Le package sf

Réaliser une intersection

m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586), c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)

Le package sf

Réaliser une intersection

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)

Le package sf

Construire des polygones de Voronoi google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r) 

mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join=st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col='lightblue')

Le package cartography

Symboles proportionnels

library(cartography)
library(sf)
# Load data
mtq <- st_read(system.file("shape/martinique.shp", package="cartography"))
# save map
png(filename = "img/map1.png", width = 400, height = 467, res = 100)
# ser margins
par(mar=c(0,0,0,0))
# Countries plot
plot(st_geometry(mtq), col = "lightblue4",
     border = "lightblue3", bg = "lightblue1")
# Population plot on proportional symbols
propSymbolsLayer(x = mtq, 
                 var = "P13_POP", 
                 legend.title.txt = "Total\npopulation (2013)")
# Title
mtext(text = "Population en Martinique", side = 3, line = -1)
dev.off()
#                                                                                         .

Symboles proportionnels

map1

Carte choroplèthe

png(filename = "img/map2.png", width = 400, height = 467, res = 100)
# Compute the compound annual growth rate
mtq$cagr <- (((mtq$P13_POP / mtq$P08_POP)^(1/4)) - 1) * 100
# set margins
par(mar=c(0,0,0,0))
# Map plot
choroLayer(x = mtq, 
           var = "cagr", 
           breaks = c(-6.14,-2,-1,0,1,2),
           col = c("#135D89", "#4D95BA", "#96D1EA", "#FCDACA", "#EC4E49"),
           legend.title.txt = "Compound annual\ngrowth rate")
# Title
mtext(text = "Evolution de la population", side = 3, line = -1)
dev.off()


#                                                                                                         .

Carte choroplèthe

map2

Palettes de couleurs

display.carto.all(20)

Palettes de couleurs

display.carto.pal("orange.pal")

Palettes de couleurs

mypal <- carto.pal(pal1 = "wine.pal", n1 = 7, pal2 = "green.pal", n2 = 12,
                   middle = TRUE, transparency = TRUE)
k <- length(mypal)
image(1:k, 1, as.matrix(1:k), col=mypal, xlab = paste(k," classes",sep=""),
      ylab = "", xaxt = "n", yaxt = "n",bty = "n")

Discrétisations

var <- mtq$cagr
moy <- mean(var)
med <- median(var)
std <- sd(var)
# Quantile intervals
breaks <- getBreaks(v = var, nclass = 6, method = "quantile")
hist(var, probability = TRUE, breaks = breaks, main="quantiles",
     col = carto.pal(pal1 = "wine.pal",3, "green.pal", 3))
rug(var)
abline(v = med, col = "blue", lwd = 3)

Discrétisations

# Mean and standard deviation (msd)
breaks <- getBreaks(v = var, method = "msd", k = 1, middle = TRUE)
hist(var, probability = TRUE, breaks = breaks, main="moyenne / écart-type",
     col = carto.pal(pal1 = "wine.pal",3 , "green.pal", 2, middle = TRUE))
rug(var)
abline(v = moy, col = "red", lwd = 3)
abline(v = moy + 0.5 * std, col = "blue", lwd = 3)
abline(v = moy - 0.5 * std, col = "blue", lwd = 3)

Combinaisons

png(filename = "img/map3.png", width = 400, height = 467, res = 100)
# set margins
par(mar=c(0,0,0,0))
# Map plot
# Countries plot
plot(st_geometry(mtq), col = "lightblue4",
     border = "lightblue3", bg = "lightblue1")
propSymbolsChoroLayer(x = mtq, var= "P13_POP", 
                 legend.var.title.txt = "Total\npopulation (2013)",
           var2 = "cagr", legend.var.pos = "bottomleft", 
           breaks = c(-6.14,-2,-1,0,1,2),
           col = c("#135D89", "#4D95BA", "#96D1EA", "#FCDACA", "#EC4E49"),
           legend.var2.title.txt = "Compound annual\ngrowth rate")
# Title
mtext(text = "Evolution de la population", side = 3, line = -1)
dev.off()


#                                                                                                         .

Combinaisons

map3

Labels

png(filename = "img/map4.png", width = 400, height = 467, res = 100)
# set margins
par(mar=c(0,0,0,0))
plot(st_geometry(mtq), col = "darkseagreen3", 
     border = "darkseagreen4", 
     bg = "#A6CAE0")
labelLayer(x = mtq, txt = "LIBGEO", col= "black", cex = 0.7, font = 4,
           halo = TRUE, bg = "white", r = 0.1, 
           overlap = FALSE, show.lines = FALSE)
dev.off()

Labels

map3

Tuiles OpenStreetMap

png(filename = "img/map4bis.png", width = 400, height = 467, res = 100)
# set margins
par(mar=c(0,0,0,0))
tiles <- getTiles(x = mtq, type = "osm", crop=T)
tilesLayer(tiles)
plot(st_geometry(mtq), add=T)
dev.off()

Tuiles OpenStreetMap

map3

Carte carroyée

library(cartography)
data(nuts2006)
nuts3.spdf@data = nuts3.df
# Create a grid layer
mygrid <- getGridLayer(x = nuts3.spdf, cellsize = 100000 * 100000, 
                       type = "hexagonal", var = c("pop1999","pop2008"))
# 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")

# save the map in png format
png(filename = "img/map5.png", width = 600, height = 535, res = 100)
par(mar = c(0, 0, 1.2, 0))
# 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 = "Package cartography v2.0.2", 
            sources = "Source: Eurostat, 2011", frame = TRUE, 
            scale = 500, north = TRUE, theme = "taupe.pal")
dev.off()



#                                                                                                         .

Carte carroyée

map5

Carte de discontinuités

# Get a SpatialLinesDataFrame of countries borders
nuts2.contig <- getBorders(x = nuts2.spdf)
nuts2.df$gdpcap <- nuts2.df$gdppps2008/nuts2.df$pop2008 * 1e+06

# save the map in png format
png(filename = "img/map6.png", width = 600, height = 535, res = 100)
par(mar = c(0, 0, 1.2, 0))

# 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 = "Package cartography v2.0.2", 
            sources = "Source: Eurostat, 2011", frame = TRUE, 
            scale = 500, north = FALSE, theme = "grey.pal")
north(pos="topleft")
dev.off()







#                                                                                                         .

Carte de discontinuités

map6

Lissage

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

# transform gdp in euros
mygrid$gdppps2008 <- mygrid$gdppps2008 * 1e+06

# list of breaks
v <- c(2920, 5000, 10000, 15000, 20000, 23500, 30000, 35000, 40000, 42720)

# save the map to a png file
# save the map in png format
png(filename = "img/map7.png", width = 600, height = 535, res = 100)
par(mar = c(0, 0, 1.2, 0))

# Plot a background layers
plot(nuts0.spdf, border = NA, col = NA, bg = "#A6CAE0")
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 = "gdppps2008", var2 = "pop2008", breaks = v, 
            span = 1e+05, beta = 2, mask = nuts0.spdf, 
            resolution = 50000, 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 = "Package 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 = 3600000, 
     labels = "Distance function:\n- type = exponential\n- beta = 2\n- span = 100 km", 
     cex = 0.8, adj = 0, font = 3)

dev.off()











#                                                                                                         .

Lissage

map7

Merci

Merci

Code & Présentation
https://github.com/riatelab/cartographie-reproductible

Télécharger les données https://github.com/riatelab/cartographie-reproductible/archive/master.zip

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

Blogs
neocarto: neocarto.hypotheses.org
rgeomatic: rgeomatic.hypotheses.org

sessionInfo()

sessionInfo()
R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

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

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

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

other attached packages:
[1] cartography_2.0.2.9000 sf_0.6-0               knitr_1.20            

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.15     magrittr_1.5     units_0.5-1      lattice_0.20-35  rlang_0.2.0     
 [6] stringr_1.3.0    udunits2_0.13    tools_3.4.3      revealjs_0.9     grid_3.4.3      
[11] e1071_1.6-8      DBI_0.7          htmltools_0.3.6  rgeos_0.3-26     class_7.3-14    
[16] yaml_2.1.17      rprojroot_1.3-2  digest_0.6.15    codetools_0.2-15 evaluate_0.10.1 
[21] rmarkdown_1.8    sp_1.2-7         stringi_1.1.6    compiler_3.4.3   pillar_1.2.1    
[26] backports_1.1.2  classInt_0.1-24