Analyses cartographiques avec R

1 Objectif

L’objectif de cette note est de présenter un projet fil conducteur qui permette de dérouler un ensemble d’opérateurs SIG sous R, vecteur comme raster. Le projet fil conducteur retenu concerne la classification d’un ensemble de parcelles viticoles sur un petit bassin versant vis-à-vis du risque d’érosion ravinaire. Il s’appuie sur un schéma classique de croisement entre aléa et vulnérabilité par parcelle suivant un système de score (théorique et virtuel pour simplifier l’exercice).

\[ Score_{risque}=Score_{vulnerabilite}*Score_{alea} \]

Nous négligerons dans l’approche développée beaucoup de questions et problèmes (effet de résolution spatiale des données, significativité physique de la démarche proposée, etc). Encore une fois elle n’est proposée ici qu’à des fins de simplification pédagogique.

1.1 Aléa

L’aléa par parcelle est résumé par son score \(Score_{alea}\), prenant une valeur de 1 à 5 suivant la classe d’appartenance par quantile des valeurs de SPI maximum par parcelle. L’indice SPI (Stream Power index) étant défini pour chaque cellule d’un MNT simplement par :

\[ SPI=A.S \]

avec \(A\), surface drainée et \(S\) pente (s.u.).

1.2 Vulnérabilité

La vulnérabilité par parcelle viticole est résumée par son score \(Score_{vulnerabilite}\), prenant également une valeur de 1 à 5, moyenne de deux indices \(Score_{NDVI}\) et \(Score_{pierosite}\).

\[Score_{vulnerabilite}=\frac{Score_{NDVI}+Score_{pierosite}}{2} \]

\(Score_{NDVI}\) prend une valeur de 1 à 5 suivant la classe d’appartenance par quantile des valeurs de NDVI (Normalized Difference Vegetation Index) moyen par parcelle. On se base sur des données images infrarouge. \(Score_{Pierosite}\) prend une valeur de 1 à 5 suivant la classe d’appartenance par quintile des valeurs de pierosité moyenne par parcelle. On se base ici sur des données terrain ponctuelles.

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 et 1 m de résolution spatiale, acquis par stéréophotogrammétrie aérienne (5 m, source CD 34) et par LiDAR (1 m, source INSU). Ils ont été pour le MNT à 1 m ‘incrustés’ (stream burning) par le réseau hydrographique afin de représenter fidèlement les schémas d’écoulement des eaux de surface sur le terrain considéré.

2.3 Images aériennes (NDVI)

Une image aérienne Infra-Rouge fausse couleur à 50 cm de résolution (BD-Ortho, Source IGN) est disponible sur la zone et sera utilisée pour dériver une carte de NDVI.

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

2.5 Réseau hydrographique

Le réseau hydrographique du bassin est constitué de fossés (agricoles, bord de route) buses, routes et chemins. Il résulte d’un relevé terrain (Source : ORE Omere). A NOTER que la fonction première de ce réseau est bien sur ce bassin son rôle anti-érosif car il permet pour une parcelle ‘aval’ d’intercepter et détourner le ruissellement de surface ‘amont’ (diminuer ce que l’on appelle la longueur de pente).

2.6 Piérosité

La piérosité sur le bassin est échantillonnée ponctuellement et est décrite au travers d’une variable qui correspond à la masse d’éléments grossiers supérieurs à 6 mm contenue dans un volume de sol de 15 cm de profondeur et d’environ 19x19 cm de surface (5.415 l). Un jeu de données de 268 échantillons sur cette variable (données AgroParisTech-INRA), variant de 0 à 3 kg environ, est disponible.

2.7 Téléchargement et préparation des données

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

https://www.dropbox.com/sh/clssw0wvwwv4isa/AACVC82ZfCPhBUnH0wBVmb0-a?dl=0

2.8 Chargement et installation des librairies

library(rgdal) # lien vers la librairie GDAL
library(rgeos) # lien vers la librairie GEOS
library(gridExtra)    # transformation d'un tableau en image

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

2.10 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_RO=readOGR("datas","vETENDUE_RO RGF93")  # etendue du bassin versant de Roujan

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/vinatier/Documents/A_SCIENTIFIQUE/ENSEIGNEMENTS/2019-03_Cartographie avec R (SILAT)/Cartes et R/datas", layer: "vETENDUE_RO RGF93"
## with 1 features
## It has 4 fields
## Integer64 fields read as strings:  CONT7III_ CONT7III_I

vETENDUE_BO=readOGR("datas","vETENDUE_BO RGF93")  # etendue du bassin versant du Bourdic

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/vinatier/Documents/A_SCIENTIFIQUE/ENSEIGNEMENTS/2019-03_Cartographie avec R (SILAT)/Cartes et R/datas", layer: "vETENDUE_BO RGF93"
## with 1 features
## It has 2 fields

vETENDUE_PE=readOGR("datas","vETENDUE_PE RGF93")  # etendue du bassin versant de la Peyne

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/vinatier/Documents/A_SCIENTIFIQUE/ENSEIGNEMENTS/2019-03_Cartographie avec R (SILAT)/Cartes et R/datas", layer: "vETENDUE_PE RGF93"
## with 1 features
## It has 4 fields

3 Import et analyse de données raster

Objectif général: Charger un fichier raster et réaliser des opérations simples.

3.1 Chargement et installation des librairies

library(raster)       # analyses sur données Raster
library(RSAGA)        # interfaçage avec SAGA (données Raster)

3.2 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

3.3 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_RO=crop(rMNT_PE,vETENDUE_RO)
rMNT_BO=crop(rMNT_PE,vETENDUE_BO)

par(mfcol=c(1,2))         # Division de la fenêtre graphique en deux cases
plot(rMNT_PE,main="MNT complet")  # Plot du MNT complet sur la première case
plot(vETENDUE_RO,border="red",add=T,lwd=3) # Ajout de l'étendue de découpe
plot(rMNT_RO,main="MNT découpé")  # Plot du MNT découpé sur la seconde case
plot(vETENDUE_RO,border="red",add=T,lwd=3) # Ajout de l'étendue de découpe

3.4 Agrégation d’un raster | aggregate

Objectif: Changer la résolution d’un raster

# Utilisation de la fonction aggregate
rMNT_RO_LowRes<-aggregate(rMNT_RO,fact=20) 

par(mfcol=c(1,2))   # Division de la fenêtre graphique en deux cases
plot(rMNT_RO,col=bpy.colors(),main="Resolution de base")
plot(rMNT_RO_LowRes,col=bpy.colors(),main="Resolution de base divisée 20 fois")

3.5 Adaptation des résolutions entre deux rasters | resample

Objectif: Adapter deux rasters avec une resolution différente

rMNT_RO_resampled=resample(rMNT_RO,rMNT_RO_LowRes)

3.6 Reprojection d’un raster | projectRaster

Objectif: Charger un nouveau système de projection et reprojeter un raster

# création d'un objet avec le système de projection WGS 84
proj_WGS84 <- CRS("+init=epsg:4326")
# reprojection du fichier raster
rMNT_RO_wgs84=projectRaster(from=rMNT_RO,crs=proj_WGS84)  

par(mfcol=c(1,2)) # Division de la fenêtre graphique en deux cases
plot(rMNT_RO,main="Système RGF 93")
plot(rMNT_RO_wgs84,main="Système WGS 84")

3.7 Enregistrement d’un raster | writeRaster

Objectif: Enregistrer un fichier raster pour qu’il soit réutilisé avec d’autres logiciels (type Qgis ou SAGA…)

writeRaster(rMNT_RO_wgs84,"MNT_wgs84.tif",overwrite=TRUE)

3.8 Opérateurs sur données Raster | terrain

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

rSLOPE_RO=terrain(rMNT_RO,opt="slope") ; names(rSLOPE_RO)="Slope"
rASP_RO=terrain(rMNT_RO,opt="aspect") ; names(rASP_RO)="Aspect"
rHILL_RO=hillShade(terrain(rMNT_RO,opt="slope"),terrain(rMNT_RO,opt="aspect"))
names(rHILL_RO)="Hillshade"
rDERIVMNT_RO=stack(rSLOPE_RO,rASP_RO,rHILL_RO)

plot(rDERIVMNT_RO)

3.9 Opérateurs sur données Raster | focal

Objectif: Lisser des données raster

FenetreLissage=matrix(1,15,15) # Création de la fenêtre de lissage
rSLOPE_RO_med<- focal(rSLOPE_RO,fun=median,w=FenetreLissage)  # Lissage avec filtre médian

par(mfcol=c(1,2)) # Division de la fenêtre graphique en deux cases
plot(rSLOPE_RO,col=bpy.colors(),main="Pente")
plot(rSLOPE_RO_med,col=bpy.colors(),main="Pente mediane")

3.10 Import et analyse de données raster multicanaux

Objectif général: Calculer, à partir d’une image multicanal de type Infrarouge couleur (IRC), un indice de végétation de type Normalized Difference Vegetation Index (NDVI).

3.11 Chargement des données multicanaux RGB et IRC | merge

Objectif: Charger des données raster multicanaux enregistrées en dalles jointives et les aligner selon une nouvelle emprise.

Ici, on charge l’orthophoto de la zone.

# identification de la liste des fichiers d'orthophoto
listOrthos=dir("datas/ORTHO/")
# ajout du chemin relatif vers les fichiers
listOrthos=paste("datas/ORTHO/",listOrthos,sep="")

rORTHOS=list() # creation d'une liste vide
# Travail sur la première bande
rORTHOS[[1]]=lapply(listOrthos,raster,band=1) # chargement des rasters
rORTHOS[[1]]=lapply(rORTHOS[[1]],aggregate,fact=5) # aggregation des rasters
rORTHOS[[1]]=lapply(rORTHOS[[1]],round)       # arrondi des valeurs
rORTHOS[[1]]=do.call(merge,rORTHOS[[1]])      # fusion des rasters

# Travail sur la deuxième bande
rORTHOS[[2]]=lapply(listOrthos,raster,band=2)
rORTHOS[[2]]=lapply(rORTHOS[[2]],aggregate,fact=5)
rORTHOS[[2]]=lapply(rORTHOS[[2]],round)
rORTHOS[[2]]=do.call(merge,rORTHOS[[2]])

# Travail sur la troisième bande
rORTHOS[[3]]=lapply(listOrthos,raster,band=3)
rORTHOS[[3]]=lapply(rORTHOS[[3]],aggregate,fact=5)
rORTHOS[[3]]=lapply(rORTHOS[[3]],round)
rORTHOS[[3]]=do.call(merge,rORTHOS[[3]])

# Empilement des rasters entre eux
rORTHO_BO=stack(rORTHOS) # creation d'un RasterStack
names(rORTHO_BO)=c("R","G","B") # nommage du RasterStack
rORTHO_RO=crop(rORTHO_BO,vETENDUE_RO) # decoupage du rasterStack

Ici, on charge l’orthophoto IRC de la zone pour calculer le NDVI:

# Creation d'une fonction pour réaliser le calcul
getDALLES=function(listFiles=listOrthos,band=1,fact=5)
                  {tmp=lapply(listFiles,raster,band=band)
                   tmp=lapply(tmp,aggregate,fact=fact)
                   tmp=lapply(tmp,round)
                   tmp=do.call(merge,tmp)
                   tmp}
# Création de la liste des fichiers IRC
listIRCs=dir("datas/IRC/")
listIRCs=paste("datas/IRC/",listIRCs,sep="")
# Application de la fonction à toutes les bandes
rIRC_BO=lapply(1:3,function(x){getDALLES(listFiles=listIRCs,band=x,fact=5)})
# Empilement des rasters entre eux et decoupage
rIRC_BO=stack(rIRC_BO)
names(rIRC_BO)=c("V","R","PIR")
rIRC_RO=crop(rIRC_BO,vETENDUE_RO)

rNDVI_RO=(rIRC_RO[["PIR"]]-rIRC_RO[["R"]])/(rIRC_RO[["PIR"]]+rIRC_RO[["R"]])
rNDVI_RO=resample(rNDVI_RO,rMNT_RO)
rNDVI_RO=focal(rNDVI_RO,fun=median,w=matrix(1,9,9))

par(mfcol=c(1,2),mar=rep(4,4))
plot(vETENDUE_RO,main="Orthophoto en RGB",border="transparent")
plotRGB(rORTHO_RO,add=T)
plot(vETENDUE_RO,main="NDVI",border="transparent")
plot(rNDVI_RO,add=T,horizontal=T)

3.12 Import de données via des serveurs extérieurs | getData

Objectif: Récupérer les données altimétriques (alt) pour différents pays. Changer la résolution et le système de projection pour qu’il corresponde au projet.

getData('ISO3') # Recuperation des codes de chaque pays

dir.create("tmp", showWarnings = FALSE)
dir.create("tmp/raster", showWarnings = FALSE)

rALT_FR=getData(name="alt",country="FRA",path = "tmp/raster/")
rALT_FR=projectRaster(from=rALT_FR,crs=proj_RGF93)
rALT_PE=crop(rALT_FR,extent(rMNT_PE))
rALT_PE=resample(rALT_PE,rMNT_PE)
diffRast=rALT_PE-rMNT_PE

par(mfcol=c(1,3))
plot(rMNT_PE)
plot(rALT_PE)
plot(diffRast)

listPAYS=c("URY","ESP","GRC","GUF")
# Creation d'une fonction pour récupérer les données d'un pays et calculer l'ombrage
calcHillshade=function(x){alt=getData(name="alt",country=x,path = "tmp/raster/")
                          hillShade(slope=terrain(alt,"slope"),
                                    aspect=terrain(alt,"aspect"))}
# Application de la fonction
rHILL_pays=lapply(listPAYS,calcHillshade)

par(mfcol=c(2,2),mar=rep(4,4))
lapply(rHILL_pays,plot)

[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

3.13 Utilisation de plateforme tierce pour l’analyse raster | RSAGA

Objectif général: Faire appel, via R et une librairie appropriée, ici RSAGA, à un logiciel de calcul de données Raster, ici SAGA (ce dernier doit être installé sur l’ordinateur au préalable http://www.saga-gis.org/en/index.html). Réaliser la détection des plateaux et des vallées sur le bassin en calculant l’index MRVBF (multiresolution index of valley bottom flatness).

3.13.1 Création du repertoire de travail et export du fichier temporaire | dir.create

# creation d'un répertoire temporaire pour stocker les couches SAGA
dir.create("tmp", showWarnings = FALSE)
dir.create("tmp/SAGA", showWarnings = FALSE)
# ecriture du raster sous forme de fichier temporaire dans l'environnement de travail
tmpSGRD<-writeRaster(x=rMNT_RO,
                     filename="tmp/SAGA/tmpSGRD",
                     format="SAGA",overwrite=T)

3.13.2 Pointage vers le logiciel SAGA | rsaga.env

if(Sys.info()['sysname']=="Linux")   
  myenv_saga=rsaga.env()

if(Sys.info()['sysname']=="Windows") 
  myenv_saga=rsaga.env(path="C:/Program Files (x86)/SAGA-GIS",
                       modules="C:/Program Files (x86)/SAGA-GIS/modules")

3.13.3 Choix de la librairie SAGA | rsaga.get.libraries

rsaga.get.libraries(path=myenv_saga$modules)
rsaga.get.lib.modules("ta_morphometry", env = myenv_saga, interactive = FALSE)
rsaga.get.usage("ta_morphometry", module=8, env = myenv_saga)

3.13.4 Calcul et interfaçage avec SAGA | rsaga.geoprocessor

rsaga.geoprocessor("ta_morphometry",8,env=myenv_saga,
                   list(DEM="tmp/SAGA/tmpSGRD.sgrd",
                        MRVBF="tmp/SAGA/tmpSAGA",
                        CLASSIFY=T,T_SLOPE=100),
                   warn=F,
                   ignore.stdout=T,
                   check.module.exists = F)

3.13.5 Récupération du résultat

resultSAGA<-raster("tmp/SAGA/tmpSAGA.sdat")

par(mfcol=c(1,2)) # Division de la fenêtre graphique en 2 cases
plot(rMNT_RO,main="MNT")
plot(resultSAGA,main="MRVBF")

4 Import et analyse de données vecteur

4.1 Chargement et installation des librairies

Objectif général:Cette partie introductive sur les opérateurs de type vecteur a pour vocation de vous familiariser avec les possibilités offertes par les librairies sp, spdep, igraph ainsi que l’utilisation d’une plateforme tierce telle que GRASS via la librairie rgrass7.

library(sp)           # analyses sur données vecteur
library(spdep)        # géostatistiques sur données vecteur
library(igraph)       # analyses de réseaux
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
library(rgrass7)      # interfaçage avec Grass 7.0 (données vecteur)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
library(doBy)         # ordre et tri de tableaux
library(gridExtra)    # transformation d'un tableau en image

4.2 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=readOGR("datas","vOCSOL_BO RGF93")
# Chargement des données de pédologie (emprise Peyne)
vPEDO_PE=readOGR("datas","vPEDO_PE RGF93")
# Chargement du réseau hydrographique (emprise Bourdic)
vHYDRO_BO=readOGR("datas","vHYDRO_BO RGF93")
# Chargement des données de pierosite (emprise Roujan)
vPIERO_RO=readOGR("datas","vPIERO_RO 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)

4.3 Calcul des aires de polygones | gArea

dOCSOL_BO_area=gArea(vOCSOL_BO,byid=T)
hist(dOCSOL_BO_area,breaks=100,xlab="Aire de chaque polygone (en m2)",main="")

4.4 Découpe d’un Spatial Object crop […,]

Objectif: Decouper des objets vecteurs selon une emprise donnée.

vHYDRO_RO=crop(vHYDRO_BO,vETENDUE_RO)
vOCSOL_RO=crop(vOCSOL_BO,vETENDUE_RO)
vPEDO_RO=crop(vPEDO_PE,vETENDUE_RO)

A NOTER: L’opérateur […,] peut également être utilisé et est équivalent à la fonction crop.

4.5 Calcul des centres de polygones coordinates

vOCSOL_BO_centres=coordinates(vOCSOL_BO) # Récupération des centres de chacun des polygones

par(mfcol=c(1,2)) # Division de la fenêtre graphique en deux cases
plot(vOCSOL_BO,axes=F,main="")
plot(vOCSOL_BO,axes=F,main="Centres de chaque parcelle")
points(vOCSOL_BO_centres,pch=".",col="red",cex=3)

4.6 Accès aux données slot des vecteurs | slot

listPolygons=slot(vOCSOL_BO, "polygons") # liste de polygones du fichier vecteur
slot(listPolygons[[1]],"area")  # aire du premier polygone
slot(listPolygons[[1]],"ID")    # identifiant du premier polygone
vecArea=sapply(listPolygons, slot, "area") # recuperation de toutes les aires des polygones

A NOTER: Les aires calculées via les slots sont mises à jour lors des opérations de découpe par la fonction crop, au contraire des données de la table attributaire.

dOCSOL_RO_area1=sapply(slot(vOCSOL_RO, "polygons"), slot, "area")
dOCSOL_RO_area2=gArea(vOCSOL_RO,byid=T)
range(dOCSOL_RO_area1-dOCSOL_RO_area2)

4.7 Manipulation de la table attributaire de données vectorielles | $

Objectif: Ajouter des colonnes supplémentaires dans la table attributaire d’un shapefile. Ici on calcule l’aire, le périmètre et le rapport Périmètre/Aire de chaque polygone.

vOCSOL_BO$Area=vecArea  # Ajout d'une colonne avec les aires de chaque polygone en m2
vOCSOL_BO$AreaHa=vecArea/10000 # Ajout des aires en Hectares

# Calcul des périmètres de chaque polygone
vOCSOL_BO_lines=as(vOCSOL_BO,"SpatialLines")# Transformation de chaque polygone en SpatialLines
vecPerimeters=SpatialLinesLengths(vOCSOL_BO_lines) # calcul des longueurs totales
vOCSOL_BO$Perimeter=vecPerimeters
vOCSOL_BO$Ratio=vecPerimeters/vecArea
head(vOCSOL_BO)

  SP_ID  Id Classvalue         Area       AreaHa Perimeter        Ratio
0     3   3       3220 5.584522e-07 5.584522e-11  93.53468 1.674891e+08
1     4   4       2212 1.080526e+03 1.080526e-01 567.37843 5.250948e-01
2    18  18       2100 5.187976e+03 5.187976e-01 437.53198 8.433578e-02
3    50  50       3100 2.153310e+01 2.153310e-03  18.72176 8.694413e-01
4   844 844       2212 7.281970e+01 7.281970e-03  48.01518 6.593708e-01
5   845 845       2212 8.643284e+03 8.643284e-01 371.66466 4.300040e-02

4.8 Reprojection de données vectorielles | spTransform

Objectif: Changer le système de projection d’un shapefile.

vOCSOL_BO_wgs84=spTransform(vOCSOL_BO,proj_WGS84)   # transformation du shapefile initial
par(mfcol=c(1,2))                       # Division de la fenêtre graphique
plot(vOCSOL_BO,main="Système Lambert 93",axes=T)
plot(vOCSOL_BO_wgs84,main="Système WGS 84",axes=T)

4.9 Analyses vecteur: overlays | gUnion gDifference + -

Objectif: Réaliser des intersections, unions et différenciations symétriques entre deux shapefiles.

vPEDOuOCSOL=gUnion(vPEDO_RO,vOCSOL_RO,byid=T,drop_lower_td=T)
vOCSOL_RO_BO_diff=gDifference(vOCSOL_BO,vETENDUE_RO,byid=T)

A NOTER: Les opérateurs arithmétiques + et - fonctionnent également pour remplacer les fonctions gUnion et gDifference, respectivement. A conseiller car dans ce cas on conserve les données des tables attributaires.

vPEDOuOCSOL=vPEDO_RO+vOCSOL_RO
vOCSOL_RO_BO_diff=vOCSOL_BO-vETENDUE_RO

  par(mfcol=c(1,4)) # Division de la fenêtre graphique
  plot(vPEDO_RO,main="Couches pedologiques",col="yellow")
  plot(vOCSOL_RO,main="Couches d'occupation du sol", col="red") 
  plot(vPEDOuOCSOL,main="Intersection OCSOL et PEDO",col="orange")
  plot(vOCSOL_RO_BO_diff,main="Différence entre OCSOL et etendue")
  plot(vETENDUE_RO,add=T,border="red",lwd=3)

A NOTER: Les fonctions gIntersection, gUnion et gDifference sont très sensibles aux erreurs de topologie dans les shapefiles. Pour y remédier, on peut opérer un nettoyage topologique rapide des parcelles selon la méthode du Buffer0:

  vOCSOL_RO=gBuffer(vOCSOL_RO,width=0,byid=TRUE)    # corrections des erreurs de topologie
  vPEDO_RO=gBuffer(vPEDO_RO,width=0,byid=TRUE)      # corrections des erreurs de topologie

4.10 Agrégation de données vectorielles | gUnaryUnion

Objectif: Agréger les polygones d’un shapefile selon les valeurs communes de sa table attributaire. Ici on agrège les occupations du sol d’un parcellaire.

IDs=vOCSOL_RO$Classvalue
IDs=factor(IDs)
vOCSOL_RO_fus=gUnaryUnion(vOCSOL_RO,id=IDs)

par(mfcol=c(1,2)) # Division de la fenêtre graphique
plot(vOCSOL_RO,main="Occupation du sol initiale")
# ajout des identifiants d'occupation du sol pour chaque parcelle
text(coordinates(vOCSOL_RO),labels=vOCSOL_RO$Classvalue,font=2,cex=0.5)
plot(vOCSOL_RO_fus,col=1:nlevels(IDs),main="Occupation du sol fusionnée")

IDs=vOCSOL_RO$Classvalue  # recuperation des occupations du sol de chaque polygone
IDs=substr(IDs,1,1)   # selection du premier niveau de Corine Land Cover
IDs=factor(IDs)      # Factorisation
vOCSOL_RO_fus2=gUnaryUnion(vOCSOL_RO,id=IDs)
par(mfcol=c(1,2)) # Division de la fenêtre graphique
plot(vOCSOL_RO,main="Occupation du sol initiale")
# ajout des identifiants d'occupation du sol pour chaque parcelle
text(coordinates(vOCSOL_RO),labels=vOCSOL_RO$Classvalue,font=2,cex=0.5)
plot(vOCSOL_RO_fus2,col=1:nlevels(IDs),main="Occupation du sol fusionnée")

4.11 Création d’une zone tampon | gBuffer

Objectif: Réaliser un polygone correspondant à la surface située à une distance maximale d’un SpatialLinesDataFrame.

Ici on applique la fonction buffer au réseau hydrographique.

vHYDRO_RO_buffer=gBuffer(vHYDRO_RO,width=30)

par(mfcol=c(1,2))
plot(vHYDRO_RO,col="blue",lwd=2)
plot(vHYDRO_RO_buffer,col="blue")

4.12 Distances entre données vectorielles | gDistance

Objectif: Calculer des distances entre deux couches vecteurs. Ici on étudiera les distances entre les polygones de l’occupation du sol et du réseau hydrographique. A NOTER: Les distances sont calculées de centroide à centroide lorsqu’il s’agit de polygones.

matDist_OCSOL=gDistance(vOCSOL_RO,byid=T)
matDist_OCSOLHYDRO=gDistance(vOCSOL_RO,vHYDRO_RO,byid=T)

Ici on sélectionne les parcelles viticoles situées entre 50 et 100 mètres du réseau hydrographique.

SelViti=vOCSOL_RO$Classvalue%in%c(2210,2211,2212)

Calcul sur la base des matrices de distances centroide a centroide:

matDist_OCSOLHYDRO_seuil=(matDist_OCSOLHYDRO>50 & matDist_OCSOLHYDRO<100)
SelSeuilDist1=apply(matDist_OCSOLHYDRO_seuil,2,any)
vOCSOL_RO_sel1=vOCSOL_RO[SelSeuilDist1 & SelViti,]

Calcul avec buffers:

vBUF_50=gBuffer(vHYDRO_RO,width=50)
vBUF_100=gBuffer(vHYDRO_RO,width=100)
vOCSOL_RO_sel2=vOCSOL_RO[SelViti,]
vOCSOL_RO_sel2=vOCSOL_RO_sel2[vBUF_100-vBUF_50,]

par(mfcol=c(1,3)) # Division de la fenêtre graphique
plot(vOCSOL_RO,main="Occupation du sol initiale")
plot(vHYDRO_RO,col="blue",add=T,lwd=2)
# ajout des identifiants d'occupation du sol pour chaque parcelle
text(coordinates(vOCSOL_RO),labels=vOCSOL_RO$Classvalue,font=2,cex=0.5)
plot(vOCSOL_RO,main="Parcelles sélectionnées par Matrice",border="grey",lwd=0.5)
plot(vOCSOL_RO_sel1,add=T,col="darkgreen")
plot(vHYDRO_RO,col="blue",add=T,lwd=2)
plot(vOCSOL_RO,main="Parcelles sélectionnées par Buffer",border="grey",lwd=0.5)
plot(vOCSOL_RO_sel2,add=T,col="darkgreen")
plot(vHYDRO_RO,col="blue",add=T,lwd=2)

4.13 Calcul de l’aire de chaque cellule d’un MNT | surfaceArea

Objectif: Calculer la surface réelle de chaque cellule d’un MNT en tenant compte de la pente de cette surface.

rMNT_RO_sp=as(rMNT_RO,"SpatialGridDataFrame")  # transformation du raster en SpatialGridDataFrame
rSA_RO=surfaceArea(rMNT_RO_sp,byCell=T)     # Calcul de l'aire de chaque cellule
rSA_RO=raster(rSA_RO)              # re-transformation du résultat en raster

par(mfcol=c(1,1))
plot(rSA_RO,main="Aire de chaque cellule en fonction du MNT")

4.14 Connectivité et topologies des données vectorielles | poly2nb

Objectif: Calculer les connectivités des parcelles entre elles. Ici on prend l’exemple des connectivités entre polygones de l’occupation du sol.

connec_nb=poly2nb(vOCSOL_RO)     #calcul des connections topologiques entre polygones
connec_centres=coordinates(vOCSOL_RO) # récupération des centres des polygones

Ici, on sélectionne les parcelles viticoles connectées aux parcelles de forêt.

IDviti=vOCSOL_RO$Classvalue%in%c(2210,2211,2212)
IDforet=(1:length(IDviti))[vOCSOL_RO$Classvalue==3100]
IDsel=sapply(connec_nb,function(x){sum(x%in%IDforet)})

par(mfcol=c(1,2),mar=rep(4,4))
plot(vOCSOL_RO,border="grey",main="Connections topologiques \n entre tous les polygones")
plot(connec_nb,connec_centres,add=T,pch=".",col="black")
points(connec_centres,col="red",pch=19)
plot(vOCSOL_RO,main="Parcelles viticoles \n connectées aux forêts")
plot(vOCSOL_RO[(IDsel & IDviti),],col="purple",add=T)
plot(vOCSOL_RO[IDforet,],col="darkgreen",add=T)
text(coordinates(vOCSOL_RO),labels=vOCSOL_RO$Classvalue,font=2,cex=0.5) # ajout des identifiants d'occupation du sol pour chaque parcelle

Ici, on compare, pour la parcelle centrale, les différences entre les distances à vol d’oiseau et les ordres de connectivité.

# Selection du polygone central
Poly_central=which.min(gDistance(vOCSOL_RO,SpatialPoints(vETENDUE_RO),byid=T))
# Calcul des distances à vol d'oiseau
DistPoly_central=gDistance(vOCSOL_RO[Poly_central,],vOCSOL_RO,byid=T)
vOCSOL_RO$Dist=as.vector(DistPoly_central)
# Calcul des ordres de connectivité
listNB=nblag(connec_nb,20)
listNB_Poly_central=sapply(listNB,function(x){x[[Poly_central]]})
vecID=sapply(listNB_Poly_central,length)
vecRepID=rep(1:length(vecID),times=vecID)
Ordre_Poly_Central=data.frame(ID=c(listNB_Poly_central,recursive=T),Ordre=vecRepID)
Ordre_Poly_Central=orderBy(~ID,Ordre_Poly_Central)
Ordre=merge(data.frame(ID=1:length(vOCSOL_RO)),Ordre_Poly_Central,all.x=T)
Ordre[is.na(Ordre[,2]),2]=0
vOCSOL_RO$Ordre=Ordre[,2]

par(mfcol=c(1,2))
plot(vOCSOL_RO,col=grey(vOCSOL_RO$Dist/max(vOCSOL_RO$Dist)),main="Distance à vol d'oiseau")
plot(vOCSOL_RO[Poly_central,],add=T,col="red")
plot(vOCSOL_RO,col=grey(vOCSOL_RO$Ordre/max(vOCSOL_RO$Ordre)),main="Distance topologique")
plot(vOCSOL_RO[Poly_central,],add=T,col="red")

4.15 Utilisation de plateforme tierce pour l’analyse vectorielle |rgrass7

4.15.1 Utilisation de la librairie rgrass7

La librairie rgrass7 offre la possibilité d’accéder à toutes les fonctions de l’environnement de travail GRASS version 7.X. Il faut également créer l’environnement de travail dans GRASS pour s’y référer par la suite. Une fois installée, on peut faire appel aux fonctions de rgrass7 de la manière suivante (on prendra l’exemple ici de la vérification et correction des topologies d’un shapefile):

4.15.2 Création du repertoire de travail | dir.create

dir.create("tmp", showWarnings = FALSE)
dir.create("tmp/GRASS", showWarnings = FALSE)

4.15.3 Pointage vers le logiciel GRASS | rsaga.env

if(Sys.info()['sysname']=="Linux") myenv_grass="/usr/lib/grass74"
if(Sys.info()['sysname']=="Windows")myenv_grass="C:/Program Files/GRASS GIS 7.0.3"

4.15.4 Définition de l’environnement de travail et export du fichier temporaire | initGRASS

#lancement du logiciel grass et affiliation au bon repertoire de travail
initGRASS(myenv_grass,home="tmp/GRASS",gisDbase="tmp/GRASS",override=T)

#reduction de la zone d'etude a celle de l'étendue de la zone
region=extent(vOCSOL_RO)
execGRASS("g.region",flags=c("p"),parameters=list(n=as.character(ymax(region)),
                                                  s=as.character(ymin(region)),
                                                  e=as.character(xmax(region)),
                                                  w=as.character(xmin(region))))
# Enregistrement du shapefile dans un fichier temporaire
writeVECT(vOCSOL_RO,"vec1",v.in.ogr_flags=c("o", "overwrite"))

4.15.5 Calcul de la fonction et dépôt du résultat dans l’environnement de travail

cpt=0
cpt=cpt+1 ; execGRASS("v.clean",parameters=list(input=paste("vec",cpt,sep=""),
                                                output=paste("vec",cpt+1,sep=""),
                                                threshold=c(0,1,5),
                                                type="line",
                                                tool="break,snap,rmdangle"),
                      flags=c("overwrite"))
cpt=cpt+1 ; execGRASS("v.build.polylines",parameters=list(input=paste("vec",cpt,sep=""),
                                                          output=paste("vec",cpt+1,sep=""),
                                                          cats="first"),
                      flags=c("overwrite"))      
cpt=cpt+1 ; execGRASS("v.clean",parameters=list(input=paste("vec",cpt,sep=""),
                                                output=paste("vec",cpt+1,sep=""),
                                                threshold=c(0,0,1,0,-1),
                                                type="line",
                                                tool="break,rmdupl,snap,break,rmdangle"),
                      flags=c("overwrite"))
cpt=cpt+1 ; execGRASS("v.generalize",parameters=list(input=paste("vec",cpt,sep=""),
                                                     output=paste("vec",cpt+1,sep=""),
                                                     method="snakes",threshold=1),flags=c("overwrite"))
cpt=cpt+1 ; execGRASS("v.type",parameters=list(input=paste("vec",cpt,sep=""),
                                               output=paste("vec",cpt+1,sep=""),
                                               from_type="line",
                                               to_type="boundary"),
                      flags=c("overwrite"))
cpt=cpt+1 ; execGRASS("v.centroids",parameters=list(input=paste("vec",cpt,sep=""),
                                                    output=paste("vec",cpt+1,sep="")),flags=c("overwrite"))

4.15.6 Récupération du résultat de calcul

OCSOL_clean<-readVECT(paste("vec",cpt+1,sep="")) # recuperation du fichier grass

5 Intégration raster-vecteur

5.1 Chargement et installation des librairies

Objectif général: L’objectif de cette partie est de vous familiariser avec les extractions de données vecteur et raster, ainsi que les fonctions de polygonisation-rasterisation disponibles.

library(akima)        # interpolation spatiale
library(maptools)

5.2 Extraction de données vecteur polygone-points | over

Objectif: Appliquer une fonction sur les données attributaires d’un SpatialPointsDataFrame séparées par les limites d’un SpatialPolygonsDataFrame.

Ici on applique la fonction over pour extraire puis faire la moyenne des collections des relevés de piérosité correspondant à chaque polygone de l’occupation du sol.

dPIERO_RO=over(vOCSOL_RO,vPIERO_RO[,"M"],fn=mean)
names(dPIERO_RO)="PIERO"

5.3 Rasterisation d’un fichier vecteur | rasterize

Objectif: Rasteriser un fichier shapefile selon une grille donnée par une couche raster. Les données présentes dans la couche raster n’ont pas d’importance.

Ici, on réalise la rasterisation des occupations du sol, l’identifiant des parcelles et du réseau hydrographique selon la résolution du MNT:

rOCSOL_RO=rasterize(vOCSOL_RO,rMNT_RO,FUN="last",field="Classvalue",progress="text")
rID_RO=rasterize(vOCSOL_RO,rMNT_RO,FUN="last",field="Id",progress="text")
rHYDRO_RO=rasterize(vHYDRO_RO,rMNT_RO,field=1,progress="text",background=NA)

par(mfcol=c(1,3))
plot(vOCSOL_RO)
plot(vHYDRO_RO,col="blue",lwd=2,add=T)
text(coordinates(vOCSOL_RO),labels=vOCSOL_RO$Classvalue,font=2,cex=0.5)
plot(rOCSOL_RO,axes=F,main="Occupation du sol rasterisée",legend=F)
plot(rHYDRO_RO,add=T,col="blue",legend=F)
plot(rID_RO,axes=F,main="Identifiants polygones rasterisés",legend=F)

5.4 Polygonisation d’un raster | rasterToPolygons

Objectif: Transformer un raster en SpatialPolygonsDataFrame

Ici, on transorme le réseau hydrographique:

vHYDRO_RO_trans=rasterToPolygons(rHYDRO_RO>0,dissolve=T)

par(mfcol=c(1,3))
plot(vHYDRO_RO,main="Réseau vecteur SpatialLines")
plot(rHYDRO_RO,main="Réseau raster")
plot(vHYDRO_RO_trans,main="Réseau vecteur SpatialPolygons")

Ici, on distingue trois zones de pentes (faible,moyen et forte):

# split des valeurs en trois modalités (faible,moyen et fort)
rSLOPE_RO_med_cut=cut(rSLOPE_RO_med,3)
# Transformation des limites entre modalités en SpatialLinesDataFrame
vSLOPE_RO_med_cut=rasterToPolygons(rSLOPE_RO_med_cut,dissolve=T)

par(mfcol=c(1,3))
plot(rSLOPE_RO_med,main="Pente médiane")
plot(rSLOPE_RO_med_cut,main="Pente médiane en trois modalités")
plot(vSLOPE_RO_med_cut,col=c("lightgrey","yellow","green"),main="Pente polygonisée")

5.5 Simplification des limites d’un SpatialPolygonsDataFrame | thinnedSpatialPoly

Objectif: Simplifier la géométrie d’un SpatialPolygonsDataFrame, en particulier après une polygonisation issue d’un raster, pour supprimer les contours en marche d’escalier.

vSLOPE_RO_med_cut_thin=thinnedSpatialPoly(vSLOPE_RO_med_cut,tolerance=20,topologyPreserve = T)

par(mfcol=c(1,2))
plot(vSLOPE_RO_med_cut,main="Polygone en marche d'escaliers")
plot(vSLOPE_RO_med_cut_thin,main="Polygone simplifié")

5.6 Extraction de données raster-polygone rasterisé | zonal

Objectif: Extraire la collection de cellules d’un raster correspondant aux identifiants d’un autre raster, puis d’appliquer une fonction sur chacune des collections. On obtient un data.frame avec le résultat de la fonction, rangé dans l’ordre d’apparition des polygones dans le shapefile.

Ici, on réalise le calcul de la pente moyenne par occupation du sol:

dSLOPE_RO=zonal(rSLOPE_RO,rOCSOL_RO,fn=mean)

plot.new()
grid.table(round(dSLOPE_RO,2))

5.7 Extraction de données raster-vecteur | extract

Objectif: Extraire les valeurs de raster correspondantes aux positions d’un SpatialPoints, SpatialPolygons ou SpatialLines.

Ici, on extrait les valeurs des dérivées du MNT correspondantes aux centres des parcelles.

extract(rDERIVMNT_RO,coordinates(vOCSOL_RO))

Ici, on extrait les valeurs des dérivées du MNT pour un polygone.

extract(rDERIVMNT_RO,vOCSOL_RO[3,])

5.8 Création d’une zone tampon | buffer

Objectif: sélectionner un ensemble de cellules d’un raster situées à une distance fixée.

Ici, on applique la fonction buffer au réseau hydrographique.

# Rasterisation du réseau hydrographique
rHYDRO_RO=rHYDRO_RO>0 # Binarisation du réseau hydrographique
rHYDRO_RO_buffer=buffer(rHYDRO_RO,width=30)  # creation d'un buffer a une distance de 30 mètres
rHYDRO_RO[is.na(rHYDRO_RO)]=0

par(mfcol=c(1,2)) # Division de la fenêtre graphique
plot(rHYDRO_RO_buffer,main="buffer à 30 mètres (RASTER)")
plot(vHYDRO_RO_buffer,main="buffer à 30 mètres (VECTEUR)")

5.9 Interpolation spatiale | interpp

Objectif: Interpolation de données ponctuelles (type SpatialPointsDataFrame) aux noeuds de la grille d’un raster suivant deux méthodes: une interpolation bilinéaire (dite interpolation TIN) et une interpolation par spline bi-cubique.

Ici, on applique l’interpolation aux données de piérosité:

vGRID_RO=rasterToPoints(rMNT_RO,spatial=T) # transformation du raster en SpatialPoints

# interpolation par spline bi-cubique des données de piérosité
vPIERO_RO_bicub=interpp(x=vPIERO_RO,z="M",xo=vGRID_RO,linear=F) 

# Rasterisation des donnees sur la grille de reference
rPIERO_RO_bicub=rasterize(vPIERO_RO_bicub,rMNT_RO,field="M") # Rasterisation des points

par(mfcol=c(1,1))
plot(rPIERO_RO_bicub,main="interpolation bicubique") ; plot(vPIERO_RO,add=T)

5.10 Calcul de la surface drainée et du Stream Power Index

Objectif général: On souhaite calculer la surface drainée par cellule d’un MNT à 5 m sur le bassin de Roujan mais en forçant les écoulements à suivre un réseau hydrographique (fossés, routes, chemins, buses. . .). Pour ce faire on va effectuer successivement 4 opérations classiques allant du “Stream Burning” (incrustation du réseau hydro dans le MNT) au calcul de surfaces drainées. Trois de ces opérations seront sous-traitées aux fonctions de SAGA appelées à l’aide de la librairie RSAGA (sous réserve que le logiciel SAGA soit installé): 1. Incrustation du réseau dans le MNT 2. Bouchage des trous du MNT incrusté 3. Calcul des directions de drainage et surfaces de drainage par cellule du MNT

5.10.1 Incrustation du réseau hydrographique dans le MNT

Objectif: Incruster le réseau hydrographique dans le MNT en abaissant la valeur du réseau rasterisé à -50 mètres:

rMNT_RO_incrust=rMNT_RO-50*rHYDRO_RO  # Incrustation du réseau

5.10.2 Bouchage des trous du MNT incrusté | rsaga.geoprocessor

Objectif: Ensuite, sur le MNT incrusté, on réalise une classique opération dite de “bouche-trous” (fill sinks) qui permet d’assurer un écoulement continu amont-aval dans la MNT. Pour cela on sous-traite cette opération via SAGA par l’intermédiaire du package RSAGA:

5.10.3 Création du repertoire de travail et export du fichier temporaire | dir.create

# creation d'un répertoire temporaire pour stocker les couches SAGA
dir.create("tmp", showWarnings = FALSE)
dir.create("tmp/SAGA", showWarnings = FALSE)
# ecriture du raster sous forme de fichier temporaire dans l'environnement de travail
tmpSGRD<-writeRaster(x=rMNT_RO_incrust,
                     filename="tmp/SAGA/tmpSGRD",
                     format="SAGA",overwrite=T)

5.10.4 Pointage vers le logiciel SAGA | rsaga.env

if(Sys.info()['sysname']=="Linux")   
  myenv_saga=rsaga.env()

if(Sys.info()['sysname']=="Windows") 
  myenv_saga=rsaga.env(path="C:/Program Files (x86)/SAGA-GIS",
                       modules="C:/Program Files (x86)/SAGA-GIS/modules")

5.10.5 Choix de la librairie SAGA | rsaga.get.libraries

rsaga.get.libraries(path=myenv_saga$modules)
rsaga.get.lib.modules("ta_preprocessor", env = myenv_saga, interactive = FALSE)
rsaga.get.usage("ta_preprocessor", module=5, env = myenv_saga)

rsaga.geoprocessor("ta_preprocessor",5,env=myenv_saga,
                   list(ELEV="tmp/SAGA/tmpSGRD.sgrd",
                        FILLED="tmp/SAGA/tmpSGRD.sgrd",
                        MINSLOPE=0.05),
                   warn=F,ignore.stdout=T,check.module.exists = FALSE)

5.10.6 Calcul de la surface drainée | rsaga.geoprocessor

Objectif: Calcul de la surface drainée par cellule en appliquant l’algorithme classique (D8) du module Catchment Area (Top-Down) de la librairie ta_hydrology.

rsaga.get.usage("ta_hydrology",0,env=myenv_saga)

rsaga.geoprocessor("ta_hydrology",0,env=myenv_saga,
                   list(ELEVATION="tmp/SAGA/tmpSGRD.sgrd",
                        METHOD=0,
                        FLOW="tmp/SAGA/tmpSGRD.sgrd"),
                   warn=F,ignore.stdout=T,check.module.exists = FALSE)

rSD_RO=raster("tmp/SAGA/tmpSGRD.sdat")

Illustration du calcul:

par(mfcol=c(1,2))
plot(rMNT_RO_incrust,main="Modèle numérique de terrain")
plot(log(rSD_RO+1),main="Surfaces drainées")

Calcul du SPI:

rSPI_RO=rSD_RO*rSLOPE_RO

par(mfcol=c(1,3),mar=rep(4,4))
plot(log(rSD_RO+1),main="Surfaces drainées")
plot(rSLOPE_RO,main="Pente")
plot(rSPI_RO,main="Stream Power Index")

5.11 Evaluation du risque

Objectif: Estimer le risque par parcelle en fonction de l’aléa et des vulnérabilités calculées tout au long du TD.

dSPI_RO=zonal(rSPI_RO,rID_RO,FUN=max)    # SPI maximal par parcelle
dNDVI_RO=zonal(rNDVI_RO,rID_RO,FUN=mean) # NDVI moyen par parcelle

# Creation de la table complète avec l'aléa et les vulnérabilités par parcelle
tabDonnees=data.frame(vOCSOL_RO$Classvalue,vOCSOL_RO$Id,dPIERO_RO)
tabDonnees=merge(tabDonnees,dSPI_RO,by.x="vOCSOL_RO.Id",by.y="zone",all.x=T)
tabDonnees=merge(tabDonnees,dNDVI_RO,by.x="vOCSOL_RO.Id",by.y="zone",all.x=T)
names(tabDonnees)=c("Id","Classvalue","PIERO","SPI","NDVI")

# Fonction de calcul des quantiles
calcQuantile=function(x){x=cut(x,breaks=quantile(x,probs=seq(0,1,0.2),na.rm=T),
                               include.lowest = T)
                         x=as.numeric(x)
                         x}
# Application de la fonction aux colonnes de la table de données
tabQuantiles=apply(tabDonnees[,-(1:2)],2,calcQuantile)

# Calcul du risque
Risque=tabQuantiles[,"SPI"]*(tabQuantiles[,"NDVI"]+tabQuantiles[,"PIERO"])/2
Risque[is.na(Risque)]=0
vOCSOL_RO$Risque=Risque

Illustration du risque par parcelle:

par(mfcol=c(1,1))
plot(vOCSOL_RO,col=grey(vOCSOL_RO$Risque/max(vOCSOL_RO$Risque)),main="Risque par parcelle")
plot(vOCSOL_RO[Risque==0,],col="red",add=T)