{"word_count": 1181, "character_count": 17322}
Le 6 novembre, je participai à Bruxelles à la journée de conférences et ateliers sur les logiciels libres en géomatique (be-OpenGIS-fr) organisé à l'Institut de Gestion de l'Environnement et d'Aménagement du Territoire (IGEAT/ULB) et l'OSGeo­fr. Ma présentation, intitulée « Cartographie et Analyse vectorielle avec R », préparée avec Timothée Giraud avait pour but de présenter le langage R auprès d'un public de géomaticiens.La présentation s'est articulée en 3 temps (Cf infra). Tous les codes sources et les données (shp et csv) sont téléchargeables sur cette page.1) Présentation de l'interface de Rstudio et de quelques instructions SIG de base.Notions clef : Import de données, calcul de corrélation, représentation graphique, import d'un shapefile (maptools), présentation de quelques fonctionnalités SIG (rgeos). Exemple d'étude stat/carto sur les délégations tunisiennes (cartographie des résidus).# =============================================# # OBJET : INTRODUCTION A L'ANALYSE SPATIALE AVEC R# # Objectif: Petite étude stat-SIG-carto# auteur : Nicolas Lambert, Timothée Giraud, Claude Grasland# version : 1.0 (oct 2014)## ===============================================# ----------------------------------------------------# R, UN LANGAGE POUR LES STATS# ----------------------------------------------------setwd("- votre repertoire de travail - ")donnees&#60;-read.csv( "data/tunisie_data_del_2011.csv",header=TRUE,sep=";",dec=",",encoding="latin1",)head(donnees)mean(donnees$POPTO2010)sd(donnees$POPTO2010)summary(donnees$POPTO2010)plot.new()boxplot(donnees$POPTO2010,main="Population totale en 2010", horizontal=TRUE)# ----------------------------------------------------# MAIS DANS R, ON PEUT AUSSI GERRE LES OBJETS SPATIAUX# ----------------------------------------------------# On charge deux packageslibrary(maptools) # package pour lire et gerer les objets spatiauxlibrary(rgeos) # http://trac.osgeo.org/geos/ (Geometry Engine Open Source)# Ouverture de deux couchesfdcOri&#60;-"geom/Tunisie_snuts4.shp"fdcOri2&#60;-"geom/sebkhas.shp"delegations&#60;-readShapeSpatial(fdcOri)sebkhas&#60;-readShapeSpatial(fdcOri2)# Affichageplot(delegations, col="#CCCCCC")plot(sebkhas, add=T,col="blue")# Tester le shpgetinfo.shape(fdcOri)class(delegations)head(gIsValid(delegations, byid = TRUE, reason=TRUE))# Acceder à la table attributairehead(delegations@data)# ---------------------------------------# GRACE A GEOS, DANS R, ON PEUT AUSSI FAIRE DU SIG# ---------------------------------------#extraction d'un polygone par simple requetepoly1&#60;-delegations[delegations@data$id=="TS3234",]plot(poly1,col="black", add=T)#Extraction des contours (gBoundary)b = gBoundary(poly1)plot(b, col="red",lwd=3,add=T)#buffer (gBuffer)buff&#60;-gBuffer(poly1, byid=TRUE, id=NULL, width=20000, quadsegs=5, capStyle="ROUND",joinStyle="ROUND", mitreLimit=1.0)plot(buff,add=TRUE,col="yellow")plot(poly1,col="black", add=T)#centroide (gCentroid)centres&#60;-gCentroid(poly1, byid=TRUE, id = NULL)plot(centres,col="red",add=T,lwd=4)head(centres@coords)#Aggregation des géométries avec gUnaryUnionhead(delegations@data)buff&#60;-gBuffer(delegations, byid=TRUE, id=NULL, width=1, quadsegs=5, capStyle="ROUND",joinStyle="ROUND", mitreLimit=1.0)gouvernorats&#60;-gUnaryUnion(buff,id = delegations@data$id_snuts3)regions&#60;-gUnaryUnion(buff, id = delegations@data$id_snuts2)macro&#60;-gUnaryUnion(buff, id = delegations@data$id_snuts1)country&#60;-gUnaryUnion(buff, id = delegations@data$id_snuts0)par(mfrow=c(1,5))plot(delegations)title(main="Délégations")plot(gouvernorats)title(main="Gouvernorats")plot(regions)title(main="Regions")plot(macro)title(main="Zones")plot(country)title(main="Pays")par(mfrow=c(1,1))# ---------------------------------------# EXEMPLE D'UTILISATION RAPIDE# Existe t il une relation entre l'indice de développement régional# Et la distance à la côte en Tunise# ---------------------------------------# ----------------------------------------------------------# Etape 1 : Ouverture du tableau de donnéesdonnees&#60;-read.csv( "data/tunisie_data_del_2011.csv",header=TRUE,sep=";",dec=",",encoding="latin1",)head(donnees)# simplification du tableau (extraction des variables pertinentes)donnees&#60;-donnees[,c("del","del_nom","IDRVA2011")]head(donnees)# ----------------------------------------------------------# etape 2 : créer une nouvelle variable : distance à la côte# Affichage de la cote et des centroidescentres&#60;-readShapeSpatial("geom/Tunisie_snuts4_centres.shp")coast&#60;-readShapeSpatial("geom/coast.shp")plot(coast,col="red", lwd=1.5)plot(centres,add=TRUE)# calcul distance points -&#62; ligne (gDistance)dist&#60;-gDistance(coast,centres,byid=TRUE)dist&#60;-data.frame(centres@data$id,dist)colnames(dist)&#60;-c("id","dist")head(dist)# Jointure entre les deux tableaux et selections des colonnes (match)mydata = data.frame(dist, donnees[match(dist[,"id"], donnees[,"del"]),])mydata&#60;-mydata[,c("id","del_nom","IDRVA2011","dist")]colnames(mydata)&#60;-c("id","nom","idr","dist")head(mydata)# ----------------------------------------------------------# etape 3 : Analyse statistique# Variable quantitative à expliquer (Y) : idr# Variable quantitative explicative (X) : dist (en log)mydata$logDist&#60;-log(mydata$dist)# résumé et visualisation stat des variablessummary(mydata$idr)summary(mydata$logDist)plot.new()par(mfrow=c(2,2))hist(mydata$idr,main="idr",breaks=10)hist(mydata$logDist,main="dist à la côte (log)",breaks=10)boxplot(mydata$idr,main="idr", horizontal=TRUE)boxplot(mydata$logDist,main="dist à la côte (log)", horizontal=TRUE)# Etude de la relation entre X et YX &#60;- mydata$logDistY &#60;- mydata$idrpar(mfrow=c(1,1))plot(X,Y, main="Relation entre X et Y", xlab="distance à la côte (log)",  ylab="Indice de développement régional", type="p",  pch=20,   cex=0.7)   cor(X,Y)cor.test(X,Y) # une relatio est signification si p-value &#60; 0.05 (5% d'erreur)MonModele &#60;- lm(Y~X)summary(MonModele)names(MonModele)abline(MonModele,  col="red")# Calcul des résuidusmydata$Yres&#60;-MonModele$residualsmydata$Yres_std&#60;-mydata$Yres/(sd(mydata$Yres))head(mydata)# ----------------------------------------------------------# etape 4 : Cartographielibrary(RColorBrewer) # jeu de couleursdisplay.brewer.all(n=NULL, type="all", select=NULL, exact.n=TRUE)library(classInt) # methodes de discretisation# chargement du fond de cartefdc &#60;- readShapeSpatial("geom/Tunisie_snuts4.shp")# id = 1ere colcodecarto&#60;-names(fdc@data)[1]# jointurefdc@data = data.frame(fdc@data, mydata[match(fdc@data[,codecarto], mydata[,"id"]),])# discretisationdistr&#60;-c(-1000,-2,-1,-0.5,0,0.5,1,2,1000)#distr &#60;- classIntervals(fdc@data$Yres_std,5,style="quantile")$brks# Choix des couleurscolours &#60;- brewer.pal(8,"RdBu")# Affectation des couleurs aux classescolMap &#60;- colours[(findInterval(fdc@data$Yres_std,distr,all.inside=TRUE))]# Plotplot(fdc, col=colMap,border="#000000",lwd=0.2)legend(x="bottomleft", legend=leglabs(round(distr,2),over="sup. ",under="inf. "), fill=colours, bty="n",pt.cex=1,cex=0.7,title="residus standardisés")title(main="Residus standardisés",sub="Auteur: Nicolas Lambert, CNRS, 2014",cex.sub=0.7)[télécharger le programme] [télécharger les données]2) Réalisation d'une carte de discontinuités de PIB sur la Belgique.Notions clef : Requêtes, jointures, géotraitements, cartographie.# =============================================# # OBJET : REALISATION D'UNE CARTE DE DISCONTINUITES DE PIB SUR LA BELGIQUE## Objectif: Faire une carte de discontinuités# auteur : Nicolas Lambert, CNRS - UMS RIATE# version : 1.0 (oct 2014)## ===============================================rm(list=ls())library(rgeos)library(maptools)library(reshape2)# -----------------------------------------------------------# [STEP 1] IMPORT DU FOND DE CARTE# -----------------------------------------------------------setwd("- votre repertoire de travail - ")# import du fond de cartefdc&#60;-readShapeSpatial("geom/nuts2_2010.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"plot(fdc)# extraction d'un pays par requete : 2 premiers caracteres de l'ID = BEhead(fdc@data$id)fdc &#60;- fdc[substr(fdc@data$id,0,2)=="BE",]fdc@data$id&#60;-as.factor(as.vector(fdc@data$id))# Affichageplot(fdc,col="#CCCCCC")title("Belgique\n(sans erreurs topologiques)")# -----------------------------------------------------------# [STEP 2] DETECTION DES REGIONS CONTIGUES# -----------------------------------------------------------# Création d'une matrice de contiguités (gIntersect)contig&#60;-gIntersects(fdc, byid = TRUE, prepared=TRUE)row.names(contig)&#60;-fdc@data$idcolnames(contig)&#60;-fdc@data$idhead(contig)# CONVERSION : Matrice -&#62; i,j,fij (grace à reshape2)contig&#60;-melt(contig,variable.name=1,value.name="fij", na.rm=TRUE) colnames(contig)&#60;-c("i","j","cij")head(contig)# Requete : On ne conserve que les couples de regions qui sont contiguescontig&#60;-contig[contig$cij==TRUE,]contig&#60;-contig[contig$i!=contig$j,]head(contig)# On crée un id unique (concatenation des deux codes)contig$id&#60;-paste(contig$i,contig$j,sep = "_")row.names(contig)&#60;-contig$idhead(contig)# -----------------------------------------------------------# [STEP 3] EXTRACTION DES FRONTIERES# -----------------------------------------------------------# On converti les polygones et polylignesbound&#60;-gBoundary(fdc, byid=TRUE, id = fdc@data$id)# On recupere la table attributiareattribut&#60;-fdc@datarow.names(attribut)&#60;-attribut$idbound&#60;-SpatialLinesDataFrame(bound,attribut, match.ID = TRUE)plot(bound[1], col="blue")levels(bound@data$id)&#60;-bound@data$id# BOUCLE : pour chaque contiguité, on extrait la ligne ou eventuellement un pointnb_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 &#60;- spChFIDs(tmp, as.character(contig$id[x]))    if (exists("borders")){borders&#60;-spRbind(borders,tmp)}    if (!exists("borders")){borders&#60;-tmp}    }  if (class(tmp)=="SpatialPoints")  {  nb_errors&#60;-nb_errors+1  if (exists("errors")){errors&#60;-spRbind(errors,tmp)}  if (!exists("errors")){errors&#60;-tmp}    plot(tmp,col="red",add=T)  }}# 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")# AFfichageplot(fdc,col="#00000050",border="white",lwd=0.1)plot(borders,col="red",lwd=2,add=T)head(borders@data)# Eventuellement, on sauvagarde du shp#writeSpatialShape(borders, "resul/BE_borders.shp", factor2char = TRUE, max_nchar=254)#if (exists("errors")){writeSpatialShape(errors, "resul/errors.shp", factor2char = TRUE, max_nchar=254)}# -----------------------------------------------------------# [STEP 4] CARTOGRAPHIE# -----------------------------------------------------------# chargement des packages utilslibrary(RColorBrewer)library(classInt)# import des données  donnees&#60;-read.csv( "data/N2.csv",header=TRUE,sep=";",dec=",",encoding="latin1",)donnees$var&#60;-donnees$PIBperinhhead(donnees)# 2. CARTE CHOROPLETHE  head(fdc@data)fdc@data = data.frame(fdc@data, donnees[match(fdc@data[,"id"], donnees[,"id"]),])nbclass&#60;-8distr &#60;- classIntervals(fdc@data$var,nbclass,style="quantile")$brkscolours &#60;- brewer.pal(nbclass,"Greens")colMap &#60;- colours[(findInterval(fdc$var,distr,all.inside=TRUE))]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. DISCONTINUITESMaxSize&#60;-10 # epaisseur de la plus grosse ligneLineColor&#60;-"red"seuil&#60;-1.2 # seules discontinuités relaitivse superieures à "seuil" sont retenues# Double jointure (valeurs de part et d'autre de chaque ligne)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"]),])# Calcul des discontinuités (var_i / var_j)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")# Triborders&#60;- borders[order(borders@data$disc,decreasing=TRUE),]# Calcul des epaisseursborders@data$size&#60;-(borders@data$disc/max(na.omit((borders@data$disc))))*MaxSize # On ne garde que les discontinuités superieure au seuilborders&#60;- borders[borders@data$disc&#62;seuil,]# Affichageplot(borders,col=LineColor, lwd=borders@data$size ,add=T)# LegenderLeg &#60;- quantile(borders@data$size,c(0,0.5,0.80,1),type=1)rVal &#60;- round(quantile(borders@data$disc,c(0,0.5,0.80,1),type=1),1)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 programme] [télécharger les données]3) Réalisation d'une carte animée sur les centres de rétention pour migrants en Europe.Notions clef : Boucle, enveloppe convexe, intersection, zones tampon, fonctions, animation# =============================================# # OBJET : CARTOGRAPHIE DES CENTRES DE RENTENTION POUR MIGRANTS## Objectif 1 : Exemple de constructuion cartographique# Objectif 2 : Faire une carte animée# auteur : Nicolas Lambert, CNRS - UMS RIATE# version : 1.0 (oct 2014)## ===============================================# -------------------------------------------------------------# REALISATION D'UNE CARTE STATIQUE# ------------------------------------------------------------library(maptools)library(rgeos)setwd("- votre repertoire de travail - ")camps &#60;- readShapeSpatial("geom/camps_time.shp")plot(camps)countries&#60;-readShapeSpatial("geom/countries.shp")plot(countries,add=T)# On definie l'étendue géographiquecamps&#60;-camps[camps@data$longitude&#62;-20,]camps&#60;-camps[camps@data$latitude&#62;-20,]plot(camps)extent&#60;-gEnvelope(camps, byid=FALSE, id = NULL)plot(extent, add=T)# On aggrege les payscountr&#60;-gBuffer(countries, byid=FALSE, id=NULL, width=1000)countr&#60;-gBuffer(countr, byid=FALSE, id=NULL, width=-1000)# On découpe les payscountr&#60;-gIntersection(countr, extent, byid=FALSE, id=NULL)# Affichageplot(countr,add=T)# On n'affiche les camps d'une seule année =&#62; ANNEE 1992head(camps@data)camps_year&#60;-camps[camps@data[,18]==1,]plot(camps_year,type=p,pch=16,col="blue",cex=1,add=T)# Creation d'un enveloppe convexehux&#60;-gConvexHull(camps_year, byid=FALSE, id = NULL) plot(hux,add=T,col="#CCCCCC50")# Ajout d'une zone tamponbuff&#60;-gBuffer(hux, byid=FALSE, id=NULL, width=100000, quadsegs=5, capStyle="ROUND",joinStyle="ROUND", mitreLimit=1.0)plot(buff,add=T,lwd=4)# Extraction des pays couverts par la zone tampon (intersection)area&#60;-gIntersection(countr, buff, byid=FALSE, id=NULL)plot(area,add=T,col="#FF000090",border="#FF0000")# Carteplot(countr)plot(extent,add=T)plot(area,add=T,col="#FF000090",border="#FF0000")plot(camps_year,type=p,pch=16,bg="red",cex=1,add=T)title("Les Centre de rétention pour migrants en Europe en 1992")# -------------------------------------------------------------# REALISATION D'UNE CARTE ANIMEE# ------------------------------------------------------------# Principe : on met notre process dans une fonction et on itèrelibrary(animation)# creation de la fonctionplotArea&#60;-function(date){  i&#60;-date-1976  par(omi=c(0,0,0,0), mgp=c(0,0,0),mar=c(0,0,0,0) , family = "D")  par(mfrow=c(1,1),cex=1,cex.lab = 0.75,cex.main=0.2,cex.axis=0.2)  y1&#60;-1000000  y2&#60;-3000000  x1&#60;-000000  x2&#60;-7000000  plot(countr)  plot(extent,add=T)  camps_year&#60;-camps[camps@data[,i]==1,]  hux&#60;-gConvexHull(camps_year, byid=FALSE, id = NULL)   buff&#60;-gBuffer(hux, byid=FALSE, id=NULL, width=100000, quadsegs=5, capStyle="ROUND",joinStyle="ROUND", mitreLimit=1.0)  area&#60;-gIntersection(countr, buff, byid=FALSE, id=NULL)  plot(area,add=T,col="#FF000090",border="#FF0000")  plot(camps_year,type=p,pch=16,bg="red",cex=1,add=T)  text(5910000, 4000000,"Principaux lieux d'enfermement\ndes migrants en ",cex = 1,family="arial",pos=4,col="#00000090",adj = c(0,0))  text(6510000, 3900000,date,cex = 3,family="arial",pos=4,col="#00000090",adj = c(0,0))  text(1600000, 300000,"(CC) Nicolas LAMBERT, 2014\nSources : http://closethecamps.org/",cex =  1,family="arial",pos=4,col="#00000080")  rect(1462625, 4243391, 7206122, 4393391, col = "black", border = "black")   text(1462625, 4300000,"Externalisation de la politique migratoire de l'Union européenne, 1985-2013",cex = 1.7,family="arial",pos=4,col="#FFFFFF",adj = c(0,0))  points(x=5870000,y=4030000,pch=16,add=T)}# Boucle pour appeler la fonction pour chaque année (chaque colonne)carte.ani &#60;- function() {  for (i in 1985:2013) {    plotArea(i)    i&#60;-i+1  }}saveGIF(carte.ani(), interval = 0.3, autobrowse = FALSE, movie.name = "carteAnim.gif",  outdir = paste(getwd(), "/resul", sep = ""), ani.height = 700, ani.width = 1000)[télécharger le programme] [télécharger les données]Référencesreferences