Ce cours a pour objectif d’utiliser un modèle de ruissellement superficiel en interaction avec les pratiques agricoles.

Il s’adresse à des étudiants ayant déjà de bonnes connaissances de l’environnement de travail alt text et souhaitant connaître la programmation orientée-objet.

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/etdt52r7acgtw1h/AAB90YgOALXo8FCJEJZSTEBRa?dl=0

2 Contexte

La mise en place de la directive cadre sur l’eau en 2000 a permis d’établir le cadre pour une politique globale communautaire dans le domaine de l’eau.

Cette directive vise à prévenir et réduire la pollution de l’eau, promouvoir son utilisation durable, protéger l’environnement, améliorer l’état des écosystèmes aquatiques et atténuer les effets des inondations et des sécheresses.

L’usage des produits phytosanitaires a provoqué une pollution des eaux superficielles et souterraines à des teneurs variables, qu’il a fallu endiguer par la mise en place de bonnes pratiques agricoles de réduction du risque de pollution des eaux.

Parmis ces bonnes pratiques agricoles (ou BPA), définies en application de la directive européenne 91/676/CEE, il existe des mesures d’aménagement du territoire pour réduire le risque de pollution des eaux, telles que la mise en place de dispositifs enherbés.

En effet, ce principe de phytoremédiation a comme principaux avantages:

  1. de jouer le rôle de filtre lorsque l’eau ruisselle,
  2. de permettre la sédimentation des particules de terre et des résidus qui s’y sont fixés,
  3. de favoriser, via la zone racinaire, la dégradation des substances actives,
  4. d’éloigner le pulvérisateur du cours d’eau, ce qui limite les risques de contamination directe.

Parmi les éléments du paysage susceptibles d’être touchés par la directive cadre sur l’eau, les fossés agricoles sont de bons cas d’études car soumis a des pratiques agricoles variées (fauche, curage, herbicides, ou absence de pratiques) et permettent, en région méditerranéenne principalement, le drainage des eaux de pluie dans les bassins versants agricoles.

3 Objectifs

L’objectif général de ce TP est d’estimer l’impact de l’organisation spatiale des pratiques agricoles à l’échelle d’un bassin versant sur la pollution de l’eau à l’exutoire du bassin.

Il portera sur un terrain d’étude principalement agricole, avec une dominante viticole, drainé par un réseau dense de fossés agricoles. Il fera appel à un système multi-agents pour simuler les comportements individuels d’un ensemble d’agriculteurs. Les flux d’eaux seront également considérés comme des agents dans le système, afin de conserver la mémoire des charges polluantes.

4 Etapes du TP sous alt text

4.1 Création des couches raster à inclure dans le modèle

La première étape consiste à définir le répertoire de travail dans lequel les données seront rangées:

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

Fenêtre RStudio

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

Fenêtre RStudio

Fenêtre RStudio

4.2 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 alt text.

proj=CRS("+init=epsg:2154")    # projection du projet (lambert 93/RGF93)
etendue=readShapePoly(paste("datas/BV_Roujan 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

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

On applique une grille des mêmes résolution et emprise que le MNT sur l’occupation du sol (rasterisation):

ocsolRast_1962 = rasterize(x = ocsol_1962, y = MNT,field = "OCC_SOL", progress = "text")
ocsolRast_2009 = rasterize(x = ocsol_2009, y = MNT,field = "OCC_SOL", progress = "text")

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

Recuperation et rasterisation des lineaires de fosses:

fosses<-readShapeLines(paste("datas/fosses_Roujan RGF93",sep=""))
projection(fosses)<-proj
fossesRast<-rasterize(fosses, MNT, field="Id",progress="text",background=0)

4.3 Préparation des données

Calcul des distances à chaque élément du paysage:

# Distance a l'exutoire
tmp=(fossesRast==466|fossesRast==4522)
tmp[tmp[]==0]=NA
Distance_Exutoire=raster::distance(tmp)
# Distance aux routes
tmp=(ocsolRast_2009==1)
tmp[tmp[]==0]=NA
Distance_Routes=raster::distance(tmp)
# Distance aux vignes
tmp=(ocsolRast_2009==7)
tmp[tmp[]==0]=NA
Distance_Vignes=raster::distance(tmp)

Incrustation des fossés dans le modèle numérique de terrain. On considère la profondeur moyenne des fossés aux alentours de 2 mètres:

mntIncrust=MNT-((fossesRast>0)*2)

Calcul de l’ombrage:

hillShade=hillShade(terrain(mntIncrust,"slope"),terrain(mntIncrust,"aspect"))
par(mfcol=c(1,2))
plot(ocsolRast_2009,col=as.character(Correspondance$Couleur),axes=F,main="Occupation du sol",legend=F)
plot(fosses,add=T,col="blue",lwd=3)
plot(hillShade,col=grey(1:100/100),axes=F,main="Ombrage",legend=F)
plot of chunk unnamed-chunk-12

plot of chunk unnamed-chunk-12

par(mfcol=c(1,3))
plot(Distance_Exutoire,legend=F,axes=F,col=grey(100:1/100),main="Distance a l'exutoire")
plot(Distance_Vignes,legend=F,axes=F,col=grey(100:1/100),main="Distance aux vignes")
plot(Distance_Routes,legend=F,axes=F,col=grey(100:1/100),main="Distance aux routes")
plot of chunk unnamed-chunk-13

plot of chunk unnamed-chunk-13

4.4 Enregistrement des données afin d’être lues par le modèle

writeRaster(mntIncrust,paste("datas/mnt",sep=""),format="ascii",overwrite=T)
class       : RasterLayer 
dimensions  : 213, 267, 56871  (nrow, ncol, ncell)
resolution  : 5, 5  (x, y)
extent      : 725300, 726635, 6264840, 6265905  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/vinatier/Documents/A-SCIENTIFIQUE/ENSEIGNEMENTS/2014-UE Simulation paysagere/TP_Pratiques_Ruissellement/datas/mnt.asc 
names       : mnt 
writeRaster(Distance_Exutoire,paste("datas/Distance_Exutoire",sep=""),format="ascii",overwrite=T)
class       : RasterLayer 
dimensions  : 213, 267, 56871  (nrow, ncol, ncell)
resolution  : 5, 5  (x, y)
extent      : 725300, 726635, 6264840, 6265905  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/vinatier/Documents/A-SCIENTIFIQUE/ENSEIGNEMENTS/2014-UE Simulation paysagere/TP_Pratiques_Ruissellement/datas/Distance_Exutoire.asc 
names       : Distance_Exutoire 
writeRaster(Distance_Vignes,paste("datas/Distance_Vignes",sep=""),format="ascii",overwrite=T)
class       : RasterLayer 
dimensions  : 213, 267, 56871  (nrow, ncol, ncell)
resolution  : 5, 5  (x, y)
extent      : 725300, 726635, 6264840, 6265905  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/vinatier/Documents/A-SCIENTIFIQUE/ENSEIGNEMENTS/2014-UE Simulation paysagere/TP_Pratiques_Ruissellement/datas/Distance_Vignes.asc 
names       : Distance_Vignes 
writeRaster(Distance_Routes,paste("datas/Distance_Routes",sep=""),format="ascii",overwrite=T)
class       : RasterLayer 
dimensions  : 213, 267, 56871  (nrow, ncol, ncell)
resolution  : 5, 5  (x, y)
extent      : 725300, 726635, 6264840, 6265905  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/vinatier/Documents/A-SCIENTIFIQUE/ENSEIGNEMENTS/2014-UE Simulation paysagere/TP_Pratiques_Ruissellement/datas/Distance_Routes.asc 
names       : Distance_Routes 
writeRaster(hillShade,paste("datas/hillShade",sep=""),format="ascii",overwrite=T)
class       : RasterLayer 
dimensions  : 213, 267, 56871  (nrow, ncol, ncell)
resolution  : 5, 5  (x, y)
extent      : 725300, 726635, 6264840, 6265905  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/vinatier/Documents/A-SCIENTIFIQUE/ENSEIGNEMENTS/2014-UE Simulation paysagere/TP_Pratiques_Ruissellement/datas/hillShade.asc 
names       : hillShade

6 Cahier des charges du TD

6.1 Comportement du modèle aux limites

Lancer le modèle sans modifier les paramètres initiaux et observer le niveau de pollution au bout de 50, 100 et 150 pas de temps. Comparer les résultats avec vos voisins.

Pourquoi y a t-il des différences?

Lancer différentes simulations en jouant sur les boutons suivants:

Paramètre Description
rain-rate L’importance des pluies.
pluie La concentration spatiale de la pluie.
agriculteurs-nb Le nombre d’agriculteurs.
Prop-polluteur La proportion d’agriculteurs utilisant des herbicides par rapport aux agriculteurs bio.
radius-agriculteur Le rayon d’action des agriculteurs.
jours_avant_traitement Le nombre de jours avant le traitement des fossés par les agriculteurs.

Le comportement du modèle semble-t-il réaliste? Lister au moins quatre points soulignant le caractère irréaliste du modèle.

6.2 Optimisation des pratiques agricoles

Selon vous, pour une proportion d’agriculteurs donnée, quelle serait la meilleure répartition des agriculteurs sur la zone pour minimiser la pollution à l’exutoire du bassin versant?

Déterminer, par la simulation, la réponse du niveau de pollution à l’exutoire selon la répartition des agriculteurs face à un gradient donné. Le gradient peut être la distance au vignoble, aux routes ou à l’exutoire. Etablir les abaques de cette réponse pour les trois gradients donnés. Optimiser la répartition spatiale des agriculteurs en fonction de ces abaques, et tester ensuite cette répartition avec le modèle.

On considère que deux structures de prescription doivent s’installer sur la zone: l’une préconise l’enherbement tandis que l’autre l’usage d’herbicides sur la zone. Elles possèdent toutes les deux le même rayon d’action, de l’ordre de 500 mètres. Où doivent se situer ces deux structures pour un minimum de pollution à l’exutoire?