{"word_count": 538, "character_count": 6758}
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é.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&#60;-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]&#60;-"id"# [STEP 2] SIMPLIFICATION DU FOND DE CARTE# (vecter -&#62; raster -&#62; 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) !!# VARIABLESresolutionc&#60;-1500resolutionl rasterr extent(r) myr&#60;-rasterize(fdc, r, field=fdc@data$id,fun='first') #plot(myr) # raster -&#62; polystest&#60;-rasterToPolygons(myr, n=4, na.rm=TRUE, digits=1, dissolve=TRUE)test test@data$id&#60;-fdc@data$idfdc&#60;-testfdc fdc@data&#60;-data.frame(fdc@data$id)colnames(fdc@data)&#60;-"id"head(fdc@data)# Sauvagarde du shpwriteSpatialShape(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éscontig&#60;-gIntersects(fdc, byid = TRUE, prepared=TRUE)row.names(contig)&#60;-fdc@data$idcolnames(contig)&#60;-fdc@data$idcontig&#60;-melt(contig,variable.name=1,value.name="fij", na.rm=TRUE)colnames(contig)&#60;-c("i","j","cij")contig&#60;-contig[contig$cij==TRUE,]contig&#60;-contig[contig$i!=contig$j,]contig$id&#60;-paste(contig$i,contig$j,sep = "_")row.names(contig)&#60;-contig$id # Polys -&#62; linesbound&#60;-gBoundary(fdc, byid=TRUE, id = fdc@data$id)attribut&#60;-fdc@datarow.names(attribut)&#60;-attribut$idbound&#60;-SpatialLinesDataFrame(bound,attribut, match.ID = TRUE) # Lines -&#62; bordersnb_errors&#60;-0for(x in 1:nrow(contig)){polyline1&#60;-bound[bound@data$id==contig$i[x],]polyline2&#60;-bound[bound@data$id==contig$j[x],]tmp&#60;-gIntersection(polyline1, polyline2)if (class(tmp)=="SpatialLines"){tmp if (exists("borders")){borders&#60;-spRbind(borders,tmp)}if (!exists("borders")){borders&#60;-tmp}}if (class(tmp)=="SpatialPoints"){nb_errors&#60;-nb_errors+1if (exists("errors")){errors&#60;-spRbind(errors,tmp)}if (!exists("errors")){errors&#60;-tmp}# plot(tmp,col="red",add=T)}}coords&#60;-data.frame(coordinates(errors)[,1],coordinates(errors)[,2])addAttrToGeom(errors, coords, match.ID=F)errors&#60;-SpatialPointsDataFrame(errors,coords, match.ID = FALSE)# Recuparation des attributsborders&#60;-SpatialLinesDataFrame(borders,contig, match.ID = TRUE)borders@data&#60;-data.frame(borders@data$id,borders@data$i,borders@data$j)colnames(borders@data)&#60;-c("id","i","j")# simplification#borders&#60;-gSimplify(borders, tol=100000, topologyPreserve=TRUE)# Sauvagarde du shpwriteSpatialShape(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&#60;-"red"seuil# 1. IMPORTSfdc names(fdc@data)[1]&#60;-"id"borders donnees&#60;-read.csv( "data/N2.csv",header=TRUE,sep=";",dec=",",encoding="latin1",)donnees$var&#60;-donnees$PIBperinh# 2. CARTE CHOROPLETHEhead(fdc@data)fdc@data = data.frame(fdc@data, donnees[match(fdc@data[,"id"], donnees[,"id"]),])nbclass&#60;-8distr 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. DISCONTINUITESborders@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&#60;-pmax(borders@data$var.1/borders@data$var,borders@data$var/borders@data$var.1)borders@data&#60;-data.frame(borders@data$id,borders@data$disc)names(borders)&#60;-c("id","disc")borders borders@data$size&#60;-(borders@data$disc/max(na.omit((borders@data$disc))))*MaxSizebordersseuil,]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)&#160;Télécharger le code source et les données