[R] Faire une carte de discontinuités avec R

Dans l’Atlas Interactif des Régions Européennes, disponible en ligne à l’adresse http://aire.ums-riate.fr, nous proposons de représenter cartographiquement l’information géographique à l’échelle des régions européennes de différentes façons, le but étant ainsi d’appréhender la complexité des réalités géographiques par des éclairages successifs. C’est le principe de la multireprésentation théorisé par cet atlas.. Parmi ces représentantions, l’Atlas en ligne propose des cartes de discontinuités dont le but spécifique est de focaliser sur les ruptures spatiales de l’indicateur représenté.

DiscWithR

Je propose dans ce billet, la méthode pour réaliser cette carte avec le langage R. Si dans le principe, la méthodologie reste assez simple, il y quelques points techniques qui nécessitent d’être soulignés. En effet, pour extraire les frontieres entre les polygones, il faut un fond de carte sans aucune erreur topologique. Le code ci dessous propose donc un moyen (un peu couteux en temps de calcul) pour corriger le fond de carte en passant par une étape de rasterisation-vectorisation. Cette étape étant assez longue, cette partie du programme permet donc d’exporter le fond de carte corrigé au format shp, ce qui permet de rendre la partie cartographique du programme totalement indépendante.

Voici le code source
(le téléchargement est proposé en bas de page)


# pgm_discontinuities.R
# Objectif: Faire une carte de discontinuités
# auteur : Nicolas Lambert, CNRS – UMS RIATE
# version : 1.0 (janvier 2014)

library(rgeos)
library(maptools)
library(raster)
library(proj4)
library(reshape2)

# [STEP 1] IMPORT DU FOND DE CARTE
# (la projection doit être définie pour utiliser la rasterisation)

rm(list=ls())
setwd(« /home/nlambert/Documents/R/carto/ »)
fdc<-readShapeSpatial(« geom/nuts2_2006.shp », proj4string=CRS(« +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs »))
colnames(fdc@data)[1]<-« id »

# [STEP 2] SIMPLIFICATION DU FOND DE CARTE
# (vecter -> raster -> vecter)
# Cette partie du programme a pour but de corriger le fond de carte.
# A la fin, il est exporté au format shp dans le repertoire resul.
# ATTENTION, TEMPS DE CALCUL ASSEZ LONG (plusieurs minutes) !!

# VARIABLES
resolutionc<-1500
resolutionl raster
r extent(r) myr<-rasterize(fdc, r, field=fdc@data$id,fun=’first’) #plot(myr) # raster -> polys
test<-rasterToPolygons(myr, n=4, na.rm=TRUE, digits=1, dissolve=TRUE)
test test@data$id<-fdc@data$id
fdc<-test
fdc fdc@data<-data.frame(fdc@data$id)
colnames(fdc@data)<-« id »
head(fdc@data)

# Sauvagarde du shp
writeSpatialShape(fdc, « resul/fdcnew.shp », factor2char = TRUE, max_nchar=254)

# [STEP 3] EXTRACTION DES FRONTIERES
# L’extraction des frontieres s’effectue en deux étapes. Premièrement, on génère un dataframe contenant les paires de polygones contigus. Deuxièmement, pour chaque paire, on extrait par intersection les frontières entre ces polygones (préalablement transformés en polylines). Si l’intersection produit une ligne, elle sera stockée dans « borders.shp », si l’intersection produit un point, il sera stocké dans « errors.shp ».

# Création d’un fichier de contiguités
contig<-gIntersects(fdc, byid = TRUE, prepared=TRUE)
row.names(contig)<-fdc@data$id
colnames(contig)<-fdc@data$id
contig<-melt(contig,variable.name=1,value.name= »fij », na.rm=TRUE)
colnames(contig)<-c(« i », »j », »cij »)
contig<-contig[contig$cij==TRUE,]
contig<-contig[contig$i!=contig$j,]
contig$id<-paste(contig$i,contig$j,sep = « _ »)
row.names(contig)<-contig$id # Polys -> lines

bound<-gBoundary(fdc, byid=TRUE, id = fdc@data$id)
attribut<-fdc@data
row.names(attribut)<-attribut$id
bound<-SpatialLinesDataFrame(bound,attribut, match.ID = TRUE) # Lines -> borders

nb_errors<-0
for(x in 1:nrow(contig))
{
polyline1<-bound[bound@data$id==contig$i[x],]
polyline2<-bound[bound@data$id==contig$j[x],]
tmp<-gIntersection(polyline1, polyline2)
if (class(tmp)== »SpatialLines »)
{
tmp if (exists(« borders »)){borders<-spRbind(borders,tmp)}
if (!exists(« borders »)){borders<-tmp}
}

if (class(tmp)== »SpatialPoints »)
{
nb_errors<-nb_errors+1
if (exists(« errors »)){errors<-spRbind(errors,tmp)}
if (!exists(« errors »)){errors<-tmp}
# plot(tmp,col= »red »,add=T)
}
}
coords<-data.frame(coordinates(errors)[,1],coordinates(errors)[,2])
addAttrToGeom(errors, coords, match.ID=F)
errors<-SpatialPointsDataFrame(errors,coords, match.ID = FALSE)

# Recuparation des attributs
borders<-SpatialLinesDataFrame(borders,contig, match.ID = TRUE)
borders@data<-data.frame(borders@data$id,borders@data$i,borders@data$j)
colnames(borders@data)<-c(« id », »i », »j »)

# simplification
#borders<-gSimplify(borders, tol=100000, topologyPreserve=TRUE)

# Sauvagarde du shp
writeSpatialShape(borders, « resul/borders.shp », factor2char = TRUE, max_nchar=254)
if (exists(« errors »)){writeSpatialShape(errors, « resul/errors.shp », factor2char = TRUE, max_nchar=254)}

# [STEP 4] CARTOGRAPHIE
# Cette partie du programme permet de realiser une carte de discontinuités.
# Elle est indépendante.
# Si vous disposez déjà d’un shp contenant les frontières entre les polygones, cette partie dui programme se suffit à elle même.

library(maptools)
library(RColorBrewer)
library(classInt)

rm(list=ls())
setwd(« /home/nlambert/Documents/R/carto/ »)

#——————————————-
# CARTOGRAPHIE
#——————————————-

MaxSize LineColor<-« red »
seuil

# 1. IMPORTS

fdc names(fdc@data)[1]<-« id »
borders donnees<-read.csv( « data/N2.csv »,header=TRUE,sep= »; »,dec= », »,encoding= »latin1″,)
donnees$var<-donnees$PIBperinh

# 2. CARTE CHOROPLETHE
head(fdc@data)
fdc@data = data.frame(fdc@data, donnees[match(fdc@data[, »id »], donnees[, »id »]),])
nbclass<-8
distr colours colMap plot(fdc, col=colMap,border= »black »,lwd=0.1)
legend(x= »topright », legend=leglabs(round(distr,1),over= »plus de »,under= »moins de »), fill=colours, bty= »n »,pt.cex=0.3,cex=0.4,title= »GDP/inh »)

# 3. DISCONTINUITES
borders@data = data.frame(borders@data, donnees[match(borders@data[, »i »], donnees[, »id »]),])
borders@data = data.frame(borders@data, donnees[match(borders@data[, »j »], donnees[, »id »]),])
borders@data$disc<-pmax(borders@data$var.1/borders@data$var,borders@data$var/borders@data$var.1)
borders@data<-data.frame(borders@data$id,borders@data$disc)
names(borders)<-c(« id », »disc »)
borders borders@data$size<-(borders@data$disc/max(na.omit((borders@data$disc))))*MaxSize
bordersseuil,]
plot(borders,col=LineColor, lwd=borders@data$size ,add=T)

rLeg rVal legend(« bottomright »,legend=rVal, lwd = rLeg, col= »red »,bty= »n »,title= »Discontinuités »,cex=0.4,pt.cex=1)
title(main= »Discontinuités de PIB/hab en 2007″,sub= »Auteur: Nicolas LAMBERT, UMS RIATE, 2014″,cex.sub=0.6)

 

Télécharger le code source et les données


Vous aimerez aussi...

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *