Ressources numériques en sciences humaines et sociales OpenEdition Nos plateformes OpenEdition Books OpenEdition Journals Hypothèses Calenda Bibliothèques OpenEdition Freemium Suivez-nous

[R] be-OpenGIS-fr

par · Publié · Mis à jour

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.

pictR

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).
tun

# =============================================
# 
# 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<-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 packages

library(maptools) # package pour lire et gerer les objets spatiaux
library(rgeos) # http://trac.osgeo.org/geos/ (Geometry Engine Open Source)

# Ouverture de deux couches

fdcOri<-"geom/Tunisie_snuts4.shp"
fdcOri2<-"geom/sebkhas.shp"
delegations<-readShapeSpatial(fdcOri)
sebkhas<-readShapeSpatial(fdcOri2)

# Affichage

plot(delegations, col="#CCCCCC")
plot(sebkhas, add=T,col="blue")

# Tester le shp
getinfo.shape(fdcOri)
class(delegations)
head(gIsValid(delegations, byid = TRUE, reason=TRUE))

# Acceder à la table attributaire
head(delegations@data)

# ---------------------------------------
# GRACE A GEOS, DANS R, ON PEUT AUSSI FAIRE DU SIG
# ---------------------------------------

#extraction d'un polygone par simple requete
poly1<-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<-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<-gCentroid(poly1, byid=TRUE, id = NULL)
plot(centres,col="red",add=T,lwd=4)
head(centres@coords)

#Aggregation des géométries avec gUnaryUnion
head(delegations@data)
buff<-gBuffer(delegations, byid=TRUE, id=NULL, width=1, quadsegs=5, capStyle="ROUND",joinStyle="ROUND", mitreLimit=1.0)
gouvernorats<-gUnaryUnion(buff,id = delegations@data$id_snuts3)
regions<-gUnaryUnion(buff, id = delegations@data$id_snuts2)
macro<-gUnaryUnion(buff, id = delegations@data$id_snuts1)
country<-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ées

donnees<-read.csv( "data/tunisie_data_del_2011.csv",header=TRUE,sep=";",dec=",",encoding="latin1",)
head(donnees)
# simplification du tableau (extraction des variables pertinentes)
donnees<-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 centroides
centres<-readShapeSpatial("geom/Tunisie_snuts4_centres.shp")
coast<-readShapeSpatial("geom/coast.shp")
plot(coast,col="red", lwd=1.5)
plot(centres,add=TRUE)

# calcul distance points -> ligne (gDistance)
dist<-gDistance(coast,centres,byid=TRUE)
dist<-data.frame(centres@data$id,dist)
colnames(dist)<-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<-mydata[,c("id","del_nom","IDRVA2011","dist")]
colnames(mydata)<-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<-log(mydata$dist)

# résumé et visualisation stat des variables
summary(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 Y
X <- mydata$logDist
Y <- mydata$idr
par(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 < 0.05 (5% d'erreur)
MonModele <- lm(Y~X)
summary(MonModele)
names(MonModele)
abline(MonModele,  col="red")

# Calcul des résuidus
mydata$Yres<-MonModele$residuals
mydata$Yres_std<-mydata$Yres/(sd(mydata$Yres))
head(mydata)

# ----------------------------------------------------------
# etape 4 : Cartographie

library(RColorBrewer) # jeu de couleurs
display.brewer.all(n=NULL, type="all", select=NULL, exact.n=TRUE)
library(classInt) # methodes de discretisation

# chargement du fond de carte
fdc <- readShapeSpatial("geom/Tunisie_snuts4.shp")
# id = 1ere col
codecarto<-names(fdc@data)[1]
# jointure
fdc@data = data.frame(fdc@data, mydata[match(fdc@data[,codecarto], mydata[,"id"]),])
# discretisation
distr<-c(-1000,-2,-1,-0.5,0,0.5,1,2,1000)
#distr <- classIntervals(fdc@data$Yres_std,5,style="quantile")$brks
# Choix des couleurs
colours <- brewer.pal(8,"RdBu")
# Affectation des couleurs aux classes
colMap <- colours[(findInterval(fdc@data$Yres_std,distr,all.inside=TRUE))]
# Plot
plot(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)

2) Réalisation d’une carte de discontinuités de PIB sur la Belgique.
Notions clef : Requêtes, jointures, géotraitements, cartographie.

bel

# =============================================
# 
# 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 carte
fdc<-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]<-"id"
plot(fdc)

# extraction d'un pays par requete : 2 premiers caracteres de l'ID = BE
head(fdc@data$id)
fdc <- fdc[substr(fdc@data$id,0,2)=="BE",]
fdc@data$id<-as.factor(as.vector(fdc@data$id))
# Affichage
plot(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<-gIntersects(fdc, byid = TRUE, prepared=TRUE)
row.names(contig)<-fdc@data$id
colnames(contig)<-fdc@data$id
head(contig)

# CONVERSION : Matrice -> i,j,fij (grace à reshape2)
contig<-melt(contig,variable.name=1,value.name="fij", na.rm=TRUE) 
colnames(contig)<-c("i","j","cij")
head(contig)

# Requete : On ne conserve que les couples de regions qui sont contigues
contig<-contig[contig$cij==TRUE,]
contig<-contig[contig$i!=contig$j,]
head(contig)

# On crée un id unique (concatenation des deux codes)
contig$id<-paste(contig$i,contig$j,sep = "_")
row.names(contig)<-contig$id
head(contig)

# -----------------------------------------------------------
# [STEP 3] EXTRACTION DES FRONTIERES
# -----------------------------------------------------------

# On converti les polygones et polylignes
bound<-gBoundary(fdc, byid=TRUE, id = fdc@data$id)

# On recupere la table attributiare
attribut<-fdc@data
row.names(attribut)<-attribut$id
bound<-SpatialLinesDataFrame(bound,attribut, match.ID = TRUE)
plot(bound[1], col="blue")
levels(bound@data$id)<-bound@data$id

# BOUCLE : pour chaque contiguité, on extrait la ligne ou eventuellement un point

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 <- spChFIDs(tmp, as.character(contig$id[x]))
    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)
  }
}

# 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")

# AFfichage
plot(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 utils
library(RColorBrewer)
library(classInt)

# import des données  
donnees<-read.csv( "data/N2.csv",header=TRUE,sep=";",dec=",",encoding="latin1",)
donnees$var<-donnees$PIBperinh
head(donnees)

# 2. CARTE CHOROPLETHE  
head(fdc@data)
fdc@data = data.frame(fdc@data, donnees[match(fdc@data[,"id"], donnees[,"id"]),])
nbclass<-8
distr <- classIntervals(fdc@data$var,nbclass,style="quantile")$brks
colours <- brewer.pal(nbclass,"Greens")
colMap <- 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. DISCONTINUITES
MaxSize<-10 # epaisseur de la plus grosse ligne
LineColor<-"red"
seuil<-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<-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")
# Tri
borders<- borders[order(borders@data$disc,decreasing=TRUE),]
# Calcul des epaisseurs
borders@data$size<-(borders@data$disc/max(na.omit((borders@data$disc))))*MaxSize 
# On ne garde que les discontinuités superieure au seuil
borders<- borders[borders@data$disc>seuil,]
# Affichage
plot(borders,col=LineColor, lwd=borders@data$size ,add=T)
# Legende
rLeg <- quantile(borders@data$size,c(0,0.5,0.80,1),type=1)
rVal <- 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)

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

carteAnim

# =============================================
# 
# 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 <- readShapeSpatial("geom/camps_time.shp")
plot(camps)
countries<-readShapeSpatial("geom/countries.shp")
plot(countries,add=T)

# On definie l'étendue géographique
camps<-camps[camps@data$longitude>-20,]
camps<-camps[camps@data$latitude>-20,]
plot(camps)
extent<-gEnvelope(camps, byid=FALSE, id = NULL)
plot(extent, add=T)

# On aggrege les pays
countr<-gBuffer(countries, byid=FALSE, id=NULL, width=1000)
countr<-gBuffer(countr, byid=FALSE, id=NULL, width=-1000)
# On découpe les pays
countr<-gIntersection(countr, extent, byid=FALSE, id=NULL)
# Affichage
plot(countr,add=T)

# On n'affiche les camps d'une seule année => ANNEE 1992
head(camps@data)
camps_year<-camps[camps@data[,18]==1,]
plot(camps_year,type=p,pch=16,col="blue",cex=1,add=T)

# Creation d'un enveloppe convexe
hux<-gConvexHull(camps_year, byid=FALSE, id = NULL) 
plot(hux,add=T,col="#CCCCCC50")

# Ajout d'une zone tampon
buff<-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<-gIntersection(countr, buff, byid=FALSE, id=NULL)
plot(area,add=T,col="#FF000090",border="#FF0000")

# Carte
plot(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ère

library(animation)

# creation de la fonction
plotArea<-function(date){  i<-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<-1000000
  y2<-3000000
  x1<-000000
  x2<-7000000
  plot(countr)
  plot(extent,add=T)

  camps_year<-camps[camps@data[,i]==1,]
  hux<-gConvexHull(camps_year, byid=FALSE, id = NULL) 
  buff<-gBuffer(hux, byid=FALSE, id=NULL, width=100000, quadsegs=5, capStyle="ROUND",joinStyle="ROUND", mitreLimit=1.0)
  area<-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 <- function() {
  for (i in 1985:2013) {
    plotArea(i)
    i<-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)

Références

references

OpenEdition vous propose de citer ce billet de la manière suivante :
Nicolas Lambert (12 novembre 2014). [R] be-OpenGIS-fr. Carnet (neo)cartographique. Consulté le 7 octobre 2024 à l’adresse https://doi.org/10.58079/rres

Nicolas Lambert

Ingénieur de recherche CNRS en sciences de l'information géographique. Membre de l'UMS RIATE et du réseau MIGREUROP / CNRS research engineer in geographical information sciences. Member of UMS RIATE and the MIGREUROP network.

Vous aimerez aussi...

2 réponses

  1. Juliette dit :
    18/04/2017 à 16:30

    Travail sympa, vous pourriez mettre le fichier des les limites administratives (uniquement des gouvernorats) digitalisées en .shp ? Si vous avez fait ce travail

    Répondre

Laisser un commentaire

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

Ce site utilise Akismet pour réduire les indésirables. En savoir plus sur comment les données de vos commentaires sont utilisées.