1 Chargement des données
2 Définition du répertoire de travail
3 Chargement et visualisation des donnees de base
4 Analyse des parcellaires
5 Création du parcellaire sans contraintes
6 Création du semis de points
7 Ajout de contraintes sur la densité des parcelles
- 7.1 Création du semis de points avec contraintes
- 7.2 Exercice pratique
8 Ajout de contraintes sur l’occupation du sol
9 Ajout de contraintes sur la dynamique d’occupation du sol
- 9.1 Application des matrices de transition (ou chaines de Markov)
- 9.2 Exercice pratique

Ce cours a pour objectifs de proposer une série d’exemple d’utilisation de R pour analyser et simuler un paysage agricole sous contraintes.

Il s’adresse à des étudiants ayant déjà de bonnes connaissances de l’environnement de travail R et souhaitant réaliser des analyses statistiques sur des données spatialisées.

1 Chargement des données

Les données nécessaires à la réalisation du TP sont disponibles à l’adresse suivante:

https://www.dropbox.com/sh/c5d0p3otid0gfo8/AAACrSz1_NCV-qDz1q4weFgsa?dl=0

2 Définition du répertoire de travail

Ici sont chargées l’ensemble des fonctions du package géomatique permettant de réaliser toutes les opérations:

library("spatstat")
library("raster")
library("maptools")
library("rgdal")
library("randomForest")

Si certains packages ne sont pas installés, utiliser la commande:

Fenêtre RStudio

Nous allons également définir le dossier dans lequel se trouvent les données, en utilisant la commande:

Fenêtre RStudio

3 Chargement et visualisation des donnees de base

Avant de commencer le projet, il est nécessaire de définir le système de projection de travail, ainsi que l’étendue de la zone d’étude. Celle-ci doit être enregistrée sous différents formats afin d’être lue par les librairies de R.

proj=CRS("+init=epsg:2154")    # projection du projet (lambert 93/RGF93)
etendue=readShapePoly(paste("datas/BV_Bourdic RGF93",sep=""))  #lecture du fichier shape
projection(etendue)<-proj              # attribution de la projection RGF 93 au fichier buffer
etendueWin=as.owin(etendue)           # transformation de l'etendue polygonale en format owin
etendueWinPSP=as.psp(etendueWin) ; marks(etendueWinPSP)=1  # transformation de l'etendue polygonale en format owin

Récupération du modèle numérique de terrain MNT qui renseigne l’altitude en tout point de l’espace:

MNT=readGDAL(paste("datas/MNT La Peyne RGF93",sep=""))  #lecture du MNT (modele numerique de terrain)

datas/MNT La Peyne RGF93 has GDAL driver GTiff 
and has 3000 rows and 3000 columns

projection(MNT)=proj        # attribution de la projection RGF 93 au MNT
MNT=raster(MNT)             # transformation en format raster (lisible par package raster)
names(MNT)="MNT"            # attribution du nom du raster
MNT=crop(MNT,etendue)       # decoupage du MNT pour qu'il corresponde a l'etendue rectangulaire du buffer

Nous allons diminuer la résolution du MNT d’un facteur 5 pour accélérer les temps de calcul:

MNT=aggregate(MNT,5)

Récupération du fichier shapefile qui renseigne l’architecture et l’occupation du sol:

ocsol_1962=readShapePoly(paste("datas/parcellaire_1962",sep=""))
ocsol_2009=readShapePoly(paste("datas/parcellaire_2009",sep=""))
projection(ocsol_1962)=proj  # Attribution de la projection RGF 93
projection(ocsol_2009)=proj  # Attribution de la projection RGF 93

L’import des fichiers shapefiles transforme automatiquement toutes les variables numériques en facteur, ce qui pose problème par la suite. Il faut donc les reconvertir en variables numériques:

ocsol_1962$OCC_SOL=as.integer(as.character(ocsol_1962$OCC_SOL))
ocsol_2009$OCC_SOL=as.integer(as.character(ocsol_2009$OCC_SOL))

Ensuite on établit la liste des occupations du sol potentielles ainsi que leur code couleur:

Correspondance=data.frame(ID=1:10,Cover=c("Route et chemins","Bois","Friche","Zone naturelle indeterminee","Vigne gobelet","Vigne palissee","Autre culture","Vigne indeterminee","Bati","Autre"),Couleur=c(rgb(167,130,116,maxColorValue=255),rgb(59,122,56,maxColorValue=255),rgb(130,222,145,maxColorValue=255),rgb(191,236,204,maxColorValue=255),rgb(255,0,255,maxColorValue=255),rgb(255,170,255,maxColorValue=255),rgb(219,181,44,maxColorValue=255),rgb(125,111,44,maxColorValue=255),rgb(170,170,170,maxColorValue=255),rgb(237,225,209,maxColorValue=255)))

Visualisation des cartes:

par(mfcol=c(1,1))
plot(MNT,col=grey(100:1/100))

plot of chunk unnamed-chunk-8

par(mfcol=c(1,1))
tmp=ocsol_1962
tmp$OCC_SOL=factor(ocsol_1962$OCC_SOL,levels=Correspondance$ID,labels=Correspondance$Cover)
spplot(tmp,col.regions=as.character(Correspondance$Couleur),axes=F,zcol="OCC_SOL",main="Année 1962")

plot of chunk unnamed-chunk-9

par(mfcol=c(1,1))
tmp=ocsol_2009
tmp$OCC_SOL=factor(ocsol_2009$OCC_SOL,levels=Correspondance$ID,labels=Correspondance$Cover)
spplot(tmp,col.regions=as.character(Correspondance$Couleur),axes=F,zcol="OCC_SOL",main="Année 2009")

plot of chunk unnamed-chunk-10

4 Analyse des parcellaires

4.1 Calcul des facteurs environnementaux

Nous allons d’abord calculer plusieurs caractéristiques géo-morphologiques sur la base du modèle numérique de terrain.

Pente=terrain(MNT,"slope")
Orientation=terrain(MNT,"aspect")
Ombrage=hillShade(Pente,Orientation) ; names(Ombrage)="ombrage"
FACTEURS=stack(MNT,Pente,Orientation,Ombrage)

par(mfcol=c(1,4),mar=rep(4,4))
plot(MNT,main="MNT",col=grey(100:1/100))
plot(Pente,main="Pente")
plot(Orientation,main="Orientation")
plot(Ombrage,col=bpy.colors(),main="Ombrage")

plot of chunk unnamed-chunk-12

4.2 Analyse des densités de parcelles

L’analyse des densités des parcelles consiste d’abord a établir la carte de chaleur des positions des centres de chaque parcelle. Une carte de chaleur correspond à un raster avec des valeurs plus élevées lorsque les points sont proches les uns des autres.

Centres=as.ppp(coordinates(ocsol_2009),etendueWin)
spatstat.options(npixel=c(ncol(MNT),nrow(MNT))) # definition de la grille de base pour les densites
CentresChaleur=raster(density(Centres))

Ensuite on représente cette carte de chaleur:

par(mfcol=c(1,2))
plot(Centres,main="",pch=19)
plot(CentresChaleur,col=bpy.colors(),main="",legend=F,axes=F,box=F)

plot of chunk unnamed-chunk-14

Enfin, on calcule l’importance de chaque caractéristique géo-morphologique sur la densité des parcelles, à l’aide d’un modèle regressif multiple:

CentresChaleur=resample(CentresChaleur,FACTEURS)  # il est nécessaire d'aligner les résolutions des cartes
datas=as.data.frame(stack(CentresChaleur,FACTEURS))
model=lm(datas)
summary(model)


Call:
lm(formula = datas)

Residuals:
      Min        1Q    Median        3Q       Max 
-6.44e-05 -1.63e-05 -2.10e-06  1.20e-05  1.04e-04 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.66e-05   6.05e-06   15.96   <2e-16 ***
MNT         -5.19e-07   1.22e-08  -42.64   <2e-16 ***
slope        1.59e-04   4.93e-06   32.23   <2e-16 ***
aspect       1.26e-06   1.27e-07    9.91   <2e-16 ***
ombrage      9.74e-05   8.15e-06   11.96   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.34e-05 on 12003 degrees of freedom
  (21481 observations deleted due to missingness)
Multiple R-squared:  0.212, Adjusted R-squared:  0.212 
F-statistic:  809 on 4 and 12003 DF,  p-value: <2e-16

4.3 Analyse de l’occupation du sol

L’analyse de l’occupation du sol consiste à comprendre comment sont réparties les cultures dans l’espace. Pour y parvenir, la méthode la plus simple consiste à rasteriser les fichiers vecteur des occupations du sol, i.e. caler une grille sur la zone et attribuer pour chaque case de la grille la valeur d’occupation du sol qui lui correspond.

ocsol_2009RAST=rasterize(ocsol_2009,MNT,field="OCC_SOL",progress="text")

Ensuite on réalise les graphiques des cartes d’occupation du sol rasterisée ou non:

par(mfcol=c(1,1),mar=rep(4,4))
tmp=ocsol_2009
tmp$OCC_SOL=factor(ocsol_2009$OCC_SOL,levels=Correspondance$ID,labels=Correspondance$Cover)
spplot(tmp,col.regions=as.character(Correspondance$Couleur),axes=F,zcol="OCC_SOL",main="Vecteur",ribbon=F)

plot of chunk unnamed-chunk-17

par(mfcol=c(1,1),mar=rep(4,4))
plot(ocsol_2009RAST,col=as.character(Correspondance$Couleur),axes=F,main="Raster",legend=F)

plot of chunk unnamed-chunk-18

Enfin, on calcule l’importance de chaque caractéristique géo-morphologique sur l’occupation du sol, à l’aide d’un modèle de type randomForest.

datasRF=na.omit(as.data.frame(stack(ocsol_2009RAST,FACTEURS)))
datasRF[,1]=factor(datasRF[,1])

On calcule le modèle randomForest:

modelRF=randomForest(layer~.,data=datasRF,importance=T)
modelRF


Call:
 randomForest(formula = layer ~ ., data = datasRF, importance = T) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 2

        OOB estimate of  error rate: 31.61%
Confusion matrix:
   1   2   3  4 5    6   7  8   9 10 class.error
1  0  11  19  2 0  161  15  2  10  1      1.0000
2  4 274  85  1 0  188  18  0  27  1      0.5418
3  1  95 420  7 0  525  58  3 103  2      0.6540
4  1   2  13 18 0   56   8  2   8  0      0.8333
5  0   0   0  0 0    0   2  0   0  0      1.0000
6  8  91 181  8 0 6087 201 27 197  0      0.1049
7  3  15  62  2 0  763 414  5  78  0      0.6915
8  1   5  12  0 0   92   4 76  14  0      0.6275
9  1  42  83  2 0  301  66  6 772  0      0.3936
10 0   2   4  0 0   21   1  0   0  6      0.8235

On peut ainsi déterminer l’importance de chacune des variables dans l’explication des occupations du sol:

par(mfcol=c(1,1),mar=rep(4,4))
varImpPlot(modelRF,type=2,main="",pch=19)

plot of chunk unnamed-chunk-21

4.4 Analyse de l’évolution temporelle des occupations du sols

On rasterise la deuxième carte datant de 1962:

ocsol_1962RAST=rasterize(ocsol_1962,MNT,field="OCC_SOL",progress="text")

On peut ainsi évaluer la matrice de transition d’une culture à une autre entre 1962 et 2009:

Trans=na.omit(crosstab(stack(ocsol_1962RAST,ocsol_2009RAST)))
vecCovers=1:10
TransMat=array(dim=c(10,10))
for(i in vecCovers)
{
   tmp=Trans[Trans[,1]==i,2:3]
   tmp=merge(vecCovers,tmp,by=1,all.x=T)
   tmp[is.na(tmp$Freq),"Freq"]=0
   tmp[,"Freq"]=tmp[,"Freq"]/sum(tmp[,"Freq"])
   TransMat[i,]=tmp[,"Freq"]
}
rownames(TransMat)=colnames(TransMat)=Correspondance$Cover

knitr::kable(round(TransMat,2))

	Route et chemins	Bois	Friche	Zone naturelle indeterminee	Vigne palissee	Autre culture	Vigne indeterminee	Bati
Route et chemins	0.72	0.01	0.00	0.00	0.03	0.00	0.00	0.22
Bois	0.01	0.91	0.04	0.00	0.01	0.00	0.00	0.03
Friche	0.00	0.34	0.31	0.02	0.17	0.03	0.00	0.12
Zone naturelle indeterminee	0.00	0.00	0.00	0.00	0.00	0.00	0.00	1.00
Vigne gobelet	0.00	0.02	0.10	0.01	0.66	0.12	0.02	0.07
Vigne palissee	0.00	0.03	0.07	0.01	0.66	0.14	0.00	0.09
Autre culture	0.00	0.00	0.18	0.01	0.47	0.19	0.05	0.09
Vigne indeterminee	0.00	0.00	0.11	0.00	0.61	0.18	0.00	0.10
Bati	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.99
Autre	0.00	0.00	0.00	0.00	0.33	0.00	0.00	0.67

4.5 Exercice pratique

Analyser la densité de parcelle et l’occupation du sol pour l’année 1962.
Calculer la fréquence de chaque occupation du sol et décrire les modifications.\
Utiliser les fonctions ?stack et ?raster::freq

5 Création du parcellaire sans contraintes

6 Création du semis de points

semis=runifpoint(100,win=etendueWin)

Graphique du résultat:

par(mfcol=c(1,2))
plot(etendueWin,main="")
plot(semis,pch=19,main="")

plot of chunk unnamed-chunk-26

6.1 Tesselation du semis de points

Tesselation=dirichlet(semis)

Graphique du résultat:

par(mfcol=c(1,2))
plot(semis,pch=19,main="")
plot(Tesselation,main="")
plot(semis,pch=19,add=T,col="darkgrey")

plot of chunk unnamed-chunk-28

6.2 Ajout des occupations du sol

Il s’agit tout d’abord de créer la table permettant d’ajouter les occupations du sol:

OccSol=data.frame(ID=1:length(tiles(Tesselation)),Annee1=NA,Annee2=NA,Annee3=NA,Annee4=NA,Annee5=NA)

Enfin on attribue de manière aléatoire les occupations du sol dans la table:

OccSol[,"Annee1"]=sample(Correspondance$Cover,size=dim(OccSol)[1],replace=T)
OccSol[,"Annee2"]=sample(Correspondance$Cover,size=dim(OccSol)[1],replace=T)
OccSol[,"Annee3"]=sample(Correspondance$Cover,size=dim(OccSol)[1],replace=T)
OccSol[,"Annee4"]=sample(Correspondance$Cover,size=dim(OccSol)[1],replace=T)
OccSol[,"Annee5"]=sample(Correspondance$Cover,size=dim(OccSol)[1],replace=T)

Graphique du résultat:

par(mfcol=c(1,5),mar=rep(0,4))
for(i in names(OccSol)[-1])
{
  plot(etendueWin,main=i)
    for(j in 1:dim(OccSol)[1])
        plot(as.owin(Tesselation[j]),col=as.character(Correspondance[Correspondance$Cover==OccSol[j,i],]$Couleur),add=T)
}

plot of chunk unnamed-chunk-31

6.3 Exercice pratique

En vous aidant de l’aide de la librairie spatstat:

Vous devrez réaliser un parcellaire représentant un pavage carré, puis un pavage rectangulaire de la zone.
Ensuite, il faudra appliquer une légère déformation à toutes les parcelles pour enlever la régularité parfaite du pavage.

7 Ajout de contraintes sur la densité des parcelles

7.1 Création du semis de points avec contraintes

Nous allons faire l’hypothèse que la densité des parcelles est dépendante de l’altitude. Plus l’altitude est élevée, et plus les parcelles sont concentrées. Il s’agit donc de réaliser un semis de point dépendant de l’altitude de la zone.

Cela nécessitera de créer une fonction permettant de transformer un objet raster en objet image:

as.rasterToim<-function(datas=MNT)
{
  tmp<-as.matrix(datas)
  win<-as.owin(c(xmin(extent(datas)),xmax(extent(datas)),ymin(extent(datas)),ymax(extent(datas))))
  tmp[dim(tmp)[1]:1,]<-tmp
  as.im(tmp,win)
}

Il faut tout d’abord transformer le mnt en une image dont les valeurs sont comprises entre 0 et 0.0001:

Rescale=function(x){
tmp<-(x-minValue(x))
tmp<-tmp/maxValue(tmp)
tmp
}
mntNorm=Rescale(MNT)*0.0001
mntNorm=as.rasterToim(mntNorm)

Ensuite on réalise un semis selon un loi de poisson contrainte par les valeurs du mnt:

semis=rpoispp(lambda=mntNorm)
semis=semis[etendueWin]

On réalise la tesselation de ce semis de points:

Tesselation=dirichlet(semis)

Enfin, le graphique du résultat

par(mfcol=c(1,3),mar=rep(1,4))
plot(mntNorm,ribbon=F,main="")
plot(mntNorm,ribbon=F,col=grey(100:1/100),main="")
plot(semis,pch=19,col="blue",add=T)
plot(etendueWin,add=T,border="red")
plot(mntNorm,ribbon=F,col=grey(100:1/100))
plot(semis,pch=19,col="blue",add=T)
plot(Tesselation,add=T,col="darkblue")

plot of chunk unnamed-chunk-37

7.2 Exercice pratique

En vous aidant de l’aide de la librairie raster et en trouvant la (ou les) fonctions adéquates:

?raster

Vous devrez réaliser un parcellaire dont les parcelles situées sur l’ubac sont plus concentrées que celles de l’adret. Nota: L’ubac désigne les versants d’une vallée de montagne qui bénéficient de la plus courte exposition au soleil. Le versant opposé s’appelle l’adret.

8 Ajout de contraintes sur l’occupation du sol

8.1 Création du parcellaire

Nous allons nous baser cette fois-ci sur un parcellaire conditionné par une distance minimale entre chaque parcelle:

semis=rSSI(r=200,win=etendueWin)
Tesselation=dirichlet(semis)

8.2 Calcul des facteurs de contraintes

A présent, nous allons considérer que la probabilité d’avoir une culture donnée dépend des caractéristiques géo-morphologiques de la zone dans laquelle elle se trouve. Les caractéristiques géo-morphologiques ont été calculées dans la partie Analyse des parcellaires

8.3 Règles de décision pour l’allocation des cultures

Ensuite on définit des règles de décision pour l’allocation de chaque culture en fonction des critères géomorphologiques. On peut se baser sur le modèle randomForest qui a été calculé dans la partie Analyse des parcellaires

On réalise la même manipulation pour le parcellaire simulé:

dataSIM=extract(FACTEURS,data.frame(semis))
vec=predict(modelRF,dataSIM,type="response")
vec=as.integer(as.character(vec))

par(mfcol=c(1,3))
plot(semis,main="Position des parcelles",pch=19)
plot(Tesselation,main="Forme des parcelles")
plot(etendueWin,main="Occupation des parcelles") ; cpt=0
for(i in vec) {cpt=cpt+1 ; plot(as.owin(Tesselation[cpt]),col=as.character(Correspondance[Correspondance$ID==i,]$Couleur),add=T)}

plot of chunk unnamed-chunk-41

8.4 Exercice pratique

Appliquer le modèle randomForest sur le fichier raster, et non plus le fichier vecteur et comparer les résultats. Vous pouvez vous aider en tapant ?raster::predict Refaites ensuite l’application du modèle pour deux résolutions différentes: 5 m et 50 m.

9 Ajout de contraintes sur la dynamique d’occupation du sol

9.1 Application des matrices de transition (ou chaines de Markov)

On se base sur la matrice de transition calculée dans la partie Analyse des parcellaires:

OccSol=data.frame(ID=1:length(tiles(Tesselation)),Annee1962=NA,Annee2009=NA,Annee2056=NA)
OccSol[,"Annee1962"]=extract(ocsol_1962RAST,data.frame(semis))

for(year in 3:4)
for(i in 1:dim(OccSol)[1])
{index=OccSol[i,year-1]
  if(!is.na(index)) OccSol[i,year]=sample(x=1:10,size=1,prob=as.vector(TransMat[index,]))
}

On visualise ensuite la dynamique des parcellaires:

par(mfcol=c(1,3))
for(year in 2:4)
{
plot(etendueWin,main=paste("Année",(year-2)*47+1962)) ; cpt=0
vec=OccSol[,year]
for(i in vec) {cpt=cpt+1 ; plot(as.owin(Tesselation[cpt]),col=as.character(Correspondance[Correspondance$ID==i,]$Couleur),add=T)}
}

plot of chunk unnamed-chunk-43

9.2 Exercice pratique

Plutôt que considérer les probabilités de passage parcelle par parcelle, on peut les appliquer directement à toutes les cellules des cartes raster. La fonction ?getValues permet d’accéder à toutes les valeurs du raster

Analyse et simulation de parcellaires agricoles

Fabrice Vinatier

01/10/2014

1 Chargement des données

2 Définition du répertoire de travail

3 Chargement et visualisation des donnees de base

4 Analyse des parcellaires

4.1 Calcul des facteurs environnementaux

4.2 Analyse des densités de parcelles

4.3 Analyse de l’occupation du sol

4.4 Analyse de l’évolution temporelle des occupations du sols

4.5 Exercice pratique

5 Création du parcellaire sans contraintes

6 Création du semis de points

6.1 Tesselation du semis de points

6.2 Ajout des occupations du sol

6.3 Exercice pratique

7 Ajout de contraintes sur la densité des parcelles

7.1 Création du semis de points avec contraintes

7.2 Exercice pratique

8 Ajout de contraintes sur l’occupation du sol

8.1 Création du parcellaire

8.2 Calcul des facteurs de contraintes

8.3 Règles de décision pour l’allocation des cultures

8.4 Exercice pratique

9 Ajout de contraintes sur la dynamique d’occupation du sol

9.1 Application des matrices de transition (ou chaines de Markov)

9.2 Exercice pratique