1 Objectif
2 Description succincte des données
3 Logiciel utilisé
4 Chargement des données
5 Import et analyse de données raster
6 Chargement des données vecteur | readOGR
- 6.1 Cartes avec données vecteur | plot
7 Analyse des parcellaires
8 Création du parcellaire sans contraintes
9 Création du semis de points
10 Ajout de contraintes sur la densité des parcelles
- 10.1 Création du semis de points avec contraintes
- 10.2 Exercice pratique
11 Ajout de contraintes sur l’occupation du sol
12 Ajout de contraintes sur la dynamique d’occupation du sol
- 12.1 Application des matrices de transition (ou chaines de Markov)
- 12.2 Exercice pratique

1 Objectif

L’objectif de ce cours interactif est de vous présenter les différents concepts permettant d’analyser et de simuler l’architecture (i.e. les limites de parcelles) et la composition (i.e. les différentes catégories d’occupation du sol) d’un parcellaire agricole.

Il portera sur différentes librairies liées au logiciel R, telles que sp, spatstat et raster.

2 Description succincte des données

2.1 Site support : Les bassins versants de Roujan, du Bourdic et de la Peyne

Le domaine spatial étudié est le bassin versant de Roujan (43.300 N, 3.190 E) d’une surface de 0.91 \(km^2\). Le bassin de Roujan est un observatoire de recherche à long-terme (ORE OMERE) équipé par l’INRA depuis 1992. La pluviométrie totale annuelle varie entre 500 et 1400 mm. Le bassin est couvert principalement de vignes avec 237 parcelles cultivées d’une surface variant de 0.03 à 2.2 ha. Ce bassin accueille également un réseau hydrographique dense (11 km). Le domaine étudié débordera cependant des strictes limites du bassin, pour limiter les effets de bord.

Situation des différents bassins versants

2.2 MNTs

Les modèles numériques de terrain (MNTs) disponibles pour ce projet sont à 5 m de résolution spatiale, acquis par stéréophotogrammétrie aérienne (5 m, source CD 34).

2.3 Images aériennes (NDVI)

Une image aérienne RGB à 50 cm de résolution (BD-Ortho, Source IGN) est disponible sur la zone et sera utilisée pour illustrer certaines parties du projet.

2.4 Occupation du sol

Une carte d’occupation du sol par parcelle a été réalisée par photo-interprétation et relevé terrain (Source : ORE Omere).

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.

3 Logiciel utilisé

4 Chargement des données

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

https://www.dropbox.com/sh/3kgclyo75h77yvr/AACmnJH7wI8atel3HqZwk9Mta?dl=0

Dossier datas à télécharger (clic droit sur le dossier) et dézipper.

4.1 Chargement et installation des librairies

library(sp)
library(spatstat)
library(raster)
library(maptools)
library(rgdal)
library(randomForest)
library(gridExtra)
library(scales)
library(fields)
library(doBy)
library(rgeos)

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

4.2 Table récapitulative des données calculées

4.3 Projection et étendue du projet

proj_RGF93=CRS("+init=epsg:2154")      # projection du projet en RGF 93
# lecture du fichier shape délimitant l'étendue du projet
vETENDUE_BO=readOGR("datas","vETENDUE_BO RGF93")  # etendue du bassin versant du Bourdic

## OGR data source with driver: ESRI Shapefile 
## Source: "datas", layer: "vETENDUE_BO RGF93"
## with 1 features
## It has 2 fields

5 Import et analyse de données raster

5.1 Chargement du fichier Raster | raster

Objectif: Charger un fichier raster de type MNT

# lecture du MNT (modele numerique de terrain)
rMNT_PE=raster("datas/rMNT_PE RGF93.tif")
projection(rMNT_PE)=proj_RGF93  # attribution de la projection RGF 93 au MNT
names(rMNT_PE)="MNT"      # Nommage du raster

5.2 Découpage d’un raster | crop

Objectif: Découper un MNT selon une étendue donnée

# découpage de la zone à l'aide de la fonction crop
rMNT_BO=crop(rMNT_PE,vETENDUE_BO)

5.3 Agrégation d’un raster | aggregate

Objectif: Changer la résolution d’un raster

# Utilisation de la fonction aggregate
rMNT_BO<-aggregate(rMNT_BO,fact=5)

5.4 Opérateurs sur données Raster | terrain

Objectif: Calculer des variables dérivées du MNT et compilation dans un même raster.

rSLOPE_BO=terrain(rMNT_BO,opt="slope") ; names(rSLOPE_BO)="Slope"
rASP_BO=terrain(rMNT_BO,opt="aspect") ; names(rASP_BO)="Aspect"
rHILL_BO=hillShade(terrain(rMNT_BO,opt="slope"),terrain(rMNT_BO,opt="aspect"))
names(rHILL_BO)="Hillshade"
rDERIVMNT_BO=stack(rSLOPE_BO,rASP_BO,rHILL_BO)

plot(rDERIVMNT_BO)

5.5 Cartes avec données raster | plot alpha

La transparence avec la fonction alpha permet de superposer des cartes, et on peut également ajouter des barres de couleur à divers endroits du graphique:

par(mfcol=c(1,1))
plot(rMNT_BO,col=bpy.colors(),main="Altitude avec ombrage en transparence",axes=F)
plot(rHILL_BO,col=alpha(grey(1:100/100),0.5),add=T,legend=F)
image.plot(legend.only=TRUE,zlim=range(na.omit(rHILL_BO[])),col=grey(1:100/100),horizontal=T)

6 Chargement des données vecteur | readOGR

Objectif: Charger des données vecteur de type Polygone, Ligne et Point et une table de correspondance associée.

# Chargement de la carte d'occupation du sol (emprise Bourdic)
vOCSOL_BO_2012=readOGR("datas","vOCSOL_BO_2012 RGF93")
vOCSOL_BO_1962=readOGR("datas","vOCSOL_BO_1962 RGF93")
# corrections des erreurs de topologie
vOCSOL_BO_2012=gBuffer(vOCSOL_BO_2012,width=0,byid=TRUE)
vOCSOL_BO_1962=gBuffer(vOCSOL_BO_1962,width=0,byid=TRUE)
# Decoupe des parcellaires selon l'étendue du bassin versant
vOCSOL_BO_2012=crop(vOCSOL_BO_2012,vETENDUE_BO)
vOCSOL_BO_1962=crop(vOCSOL_BO_1962,vETENDUE_BO)
# Chargement des données de pédologie (emprise Peyne)
vPEDO_PE=readOGR("datas","vPEDO_PE RGF93")

Ici le tableau de correspondance entre les occupations du sol et leurs identifiants:

Correspondance=
  data.frame(ID=c(1220,3100,3220,2211,2212,2100,2210,1110,2220,1100),
             Cover=c("Réseaux routiers et chemins (tous)","Bois",
                     "Landes et broussailles","Vignoble en gobelet","Vignoble palissé",
                     "Autres terres arables","Vignoble en transition",
                     "Tissu urbain continu","Arboriculture","Zones urbanisées"),
             CoverEng=c("Road and rail networks","Forests","Sclerophyllous vegetation",
                        "Goblet vineyard","Aligned vineyard","Other arable lands",
                        "Transitional vineyard","Continuous urban area","Orchards",
                        "Discontinuous urban area"),
             Color=c("#8B0000","#4DFF00","#A6FF80","#B414F1","#FFAAFF","#E0EC05",
                     "#FF00FF","#00004D","#DAA520","#CD0000"))
Correspondance=unique(orderBy(~ID,Correspondance))  # Tri de la table par ordre croissant de ID

On peut afficher le tableau de correspondance comme une image:

plot.new()
grid.table(Correspondance)

6.1 Cartes avec données vecteur | plot

Il convient au départ pour des variables factorielles, ou qualitatives, de créer un vecteur de couleur correspondant aux attributs à représenter:

fOCSOL_BO_2012=factor(vOCSOL_BO_2012$Classvalue,
                      levels=Correspondance$ID,
                      labels=Correspondance$Color)
fOCSOL_BO_2012=as.character(fOCSOL_BO_2012)

par(mfrow=c(1,2)) # Division de la fenêtre graphique
plot(vETENDUE_BO,border="transparent")
plot(vOCSOL_BO_2012,col=alpha(fOCSOL_BO_2012,0.7),add=T,lwd=0.1) # plot de l'occupation du sol
plot(vETENDUE_BO,add=T,lwd=3,lty=3,border="red") # plot de l'étendue du projet
plot.new()
legend("center",legend=Correspondance$Cover,bty="n",fill=as.character(Correspondance$Color))

7 Analyse des parcellaires

7.1 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.

# Recuperation des centres des parcelles
vCENTRES_BO=as.ppp(coordinates(vOCSOL_BO_2012),vETENDUE_BO)
# definition de la grille de base pour les densites
spatstat.options(npixel=c(ncol(rMNT_BO),nrow(rMNT_BO)))
# Calcul de la carte de chaleur
rCENTRES_BO=raster(density(vCENTRES_BO))
names(rCENTRES_BO)="Centres"

Ensuite on représente cette carte de chaleur:

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

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:

# Realignement des raster pour correspondance parfaite
rCENTRES_BO=resample(rCENTRES_BO,rDERIVMNT_BO)
# Creation du tableau complet
datas=as.data.frame(stack(rCENTRES_BO,rMNT_BO,rDERIVMNT_BO))
# Suppression des valeurs manquantes
datas=na.omit(datas)
# Creation du modèle
model=lm(Centres~.,datas)
summary(model)


Call:
lm(formula = Centres ~ ., data = datas)

Residuals:
       Min         1Q     Median         3Q        Max 
-6.046e-05 -1.662e-05 -2.392e-06  9.700e-06  1.092e-04 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.868e-05  6.095e-06   16.19   <2e-16 ***
MNT         -4.396e-07  1.225e-08  -35.88   <2e-16 ***
Slope        1.669e-04  4.961e-06   33.64   <2e-16 ***
Aspect       1.370e-06  1.281e-07   10.70   <2e-16 ***
Hillshade    8.882e-05  8.202e-06   10.83   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.354e-05 on 12003 degrees of freedom
Multiple R-squared:  0.1907,    Adjusted R-squared:  0.1905 
F-statistic: 707.3 on 4 and 12003 DF,  p-value: < 2.2e-16

7.2 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.

rOCSOL_BO_2012=rasterize(vOCSOL_BO_2012,rMNT_BO,field="Classvalue",progress="text")
projection(rOCSOL_BO_2012)=proj_RGF93

par(mfrow=c(1,3)) # Division de la fenêtre graphique
plot(vETENDUE_BO,border="transparent",main="Vecteur")
plot(vOCSOL_BO_2012,col=alpha(fOCSOL_BO_2012,0.7),add=T,lwd=0.1) # plot de l'occupation du sol
plot(vETENDUE_BO,add=T,lwd=3,lty=3,border="red") # plot de l'étendue du projet
plot(rOCSOL_BO_2012,col=as.character(Correspondance$Color),axes=F,main="Raster",legend=F)
plot(vETENDUE_BO,add=T,lwd=3,lty=3,border="red") # plot de l'étendue du projet
plot.new()
legend("center",legend=Correspondance$Cover,bty="n",fill=as.character(Correspondance$Color))

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(rOCSOL_BO_2012,rMNT_BO,rDERIVMNT_BO)))
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: 33.41%
Confusion matrix:
     1100 1110 1220 2100 2210 2211 2212 2220 3100 3220 class.error
1100   40    9    0   18    0    0   46    3    3    7   0.6825397
1110   11  745    0   93    6    0  295    5   43   33   0.3948010
1220    4   11    0   22    1    0  157    4   15    8   1.0000000
2100   10   88    4  782    1    0 1011   16   51   47   0.6109453
2210    0   14    0    1   51    0   48    0    1    0   0.5565217
2211    0    0    0    1    0    0    1    0    0    0   1.0000000
2212   16  189   12  402   21    0 5774   29  111   59   0.1268713
2220    1   13    3   21    0    0  104   33   17   13   0.8390244
3100    4   25    5   48    0    0  297    9  328   32   0.5614973
3220    7   64    0   77    0    0  215    8   61  121   0.7811935

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)

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

On rasterise la deuxième carte datant de 1962:

rOCSOL_BO_1962=rasterize(vOCSOL_BO_1962,rMNT_BO,field="Classvalue",progress="text")
projection(rOCSOL_BO_1962)=proj_RGF93

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

dTRANS_BO=crosstab(stack(rOCSOL_BO_1962,rOCSOL_BO_2012))
dTRANS_BO=xtabs(Freq~Var1+Var2,dTRANS_BO)

plot.new()
grid.table(dTRANS_BO)

7.4 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

8 Création du parcellaire sans contraintes

9 Création du semis de points

vPLOT_BO=runifpoint(100,win=vETENDUE_BO)

par(mfcol=c(1,2))
plot(vETENDUE_BO,main="")
plot(vPLOT_BO,pch=19,main="")

9.1 Tesselation du semis de points

vTESS_BO=dirichlet(vPLOT_BO)

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

9.2 Ajout des occupations du sol

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

vOCSOL_BO_sim=data.frame(ID=1:length(tiles(vTESS_BO)),
                         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:

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

par(mfrow=c(1,5)) # Division de la fenêtre graphique
for(i in paste("Annee",1:5,sep=""))
{
  fOCSOL_BO=factor(vOCSOL_BO_sim[,i],
                   levels=Correspondance$ID,
                   labels=Correspondance$Color)
  fOCSOL_BO=as.character(fOCSOL_BO)
plot(vETENDUE_BO,border="transparent",main=i)
for(j in 1:length(tiles(vTESS_BO)))  plot(as.owin(vTESS_BO[j]),col=fOCSOL_BO[j],add=T)
}

9.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.

10 Ajout de contraintes sur la densité des parcelles

10.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=rMNT_BO)
{
  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 également créer une autre fonction de normalisation des données:

Rescale=function(x){
tmp<-(x-minValue(x))
tmp<-tmp/maxValue(tmp)
tmp
}

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

rMNT_BO_norm=Rescale(rMNT_BO)*0.0002
rMNT_BO_norm=as.rasterToim(rMNT_BO_norm)

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

vPLOT_BO_mnt=rpoispp(lambda=rMNT_BO_norm)

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

vTESS_BO_mnt=dirichlet(vPLOT_BO_mnt)

Enfin, le graphique du résultat

par(mfcol=c(1,3),mar=rep(1,4))
plot(rMNT_BO_norm,ribbon=F,main="")
plot(rMNT_BO_norm,ribbon=F,col=grey(100:1/100),main="")
plot(vPLOT_BO_mnt,pch=19,col="blue",add=T)
plot(vETENDUE_BO,add=T,border="red")
plot(rMNT_BO_norm,ribbon=F,col=grey(100:1/100))
plot(vPLOT_BO_mnt,pch=19,col="blue",add=T)
plot(vTESS_BO_mnt,add=T,col="darkblue")

10.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.

11 Ajout de contraintes sur l’occupation du sol

11.1 Création du parcellaire

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

vPLOT_BO_distmin=rSSI(r=200,win=vETENDUE_BO)
vTESS_BO_distmin=dirichlet(vPLOT_BO_distmin)

11.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

11.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(stack(rMNT_BO,rDERIVMNT_BO),data.frame(vPLOT_BO_distmin))
vec=predict(modelRF,dataSIM,type="response")
vec=as.integer(as.character(vec))

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

11.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.

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

12.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:

vOCSOL_BO_sim_distmin=data.frame(ID=1:length(tiles(vTESS_BO_distmin)),
                                 Annee1962=NA,Annee2009=NA,Annee2056=NA)
vOCSOL_BO_sim_distmin[,"Annee1962"]=extract(rOCSOL_BO_1962,data.frame(vPLOT_BO_distmin))

for(year in 3:4)
  for(i in 1:dim(vOCSOL_BO_sim_distmin)[1])
{index=as.character(vOCSOL_BO_sim_distmin[i,year-1])
  if(!is.na(index)) vOCSOL_BO_sim_distmin[i,year]=sample(x=Correspondance$ID,
                                                         size=1,
                                                         prob=as.vector(dTRANS_BO[index,]))
}

On visualise ensuite la dynamique des parcellaires:

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

12.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/12/2015