Optimisation des pratiques agricoles avec R et NETLOGO
Fabrice Vinatier
03/12/2015
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 
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/3kgclyo75h77yvr/AACmnJH7wI8atel3HqZwk9Mta?dl=0
Dossier datas à télécharger (clic droit sur le dossier) et dézipper.
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:
- de jouer le rôle de filtre lorsque l’eau ruisselle,
- de permettre la sédimentation des particules de terre et des résidus qui s’y sont fixés,
- de favoriser, via la zone racinaire, la dégradation des substances actives,
- 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 R
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("rgeos")
library("doBy")
library("gridExtra")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
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 du Bourdic## OGR data source with driver: ESRI Shapefile
## Source: "datas", layer: "vETENDUE_RO RGF93"
## with 1 features
## It has 4 fields# 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
# découpage de la zone à l'aide de la fonction crop
rMNT_RO=crop(rMNT_PE,vETENDUE_RO)
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"vOCSOL_RO=readOGR("datas","vOCSOL_BO_2012 RGF93") # Chargement de la carte d'occupation du sol
vOCSOL_RO=gBuffer(vOCSOL_RO,width=0,byid=TRUE) # corrections des erreurs de topologie
vOCSOL_RO=crop(vOCSOL_RO,vETENDUE_RO) # Decoupe des parcellaires selon l'étendue du bassin versant
vHYDRO_RO=readOGR("datas","vHYDRO_BO RGF93") # Chargement du reseau hydrographique
vHYDRO_RO=crop(vHYDRO_RO,vETENDUE_RO) # Decoupe des parcellaires selon l'étendue du bassin versantCorrespondance=
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
# Rasterisation de l'occupation du sol
rOCSOL_RO=rasterize(vOCSOL_RO,rMNT_RO,field="Classvalue",progress="text")
projection(rOCSOL_RO)=proj_RGF93
# Rasterisation du reseau hydrographique
rHYDRO_RO=rasterize(vHYDRO_RO,rMNT_RO,field="Id",progress="text",background=0)
projection(rHYDRO_RO)=proj_RGF93# Rasterisation des identifiants des parcelles
rOCSOL_RO_ID=rasterize(vOCSOL_RO,rMNT_RO,field="Id",progress="text")
projection(rOCSOL_RO_ID)=proj_RGF93
# Rasterisation des identifiants des fosses
rHYDRO_RO_ID=rasterize(vHYDRO_RO,rMNT_RO,field="Id",progress="text",background=NA)
projection(rHYDRO_RO_ID)=proj_RGF934.2 Préparation des données
Calcul des distances à chaque élément du paysage:
# Distance a l'exutoire
tmp=(rHYDRO_RO==466|rHYDRO_RO==4522)
tmp[tmp[]==0]=NA
rDIST_RO_exutoire=raster::distance(tmp)
# Distance aux routes
tmp=(rOCSOL_RO==1220)
tmp[tmp[]==0]=NA
rDIST_RO_routes=raster::distance(tmp)
# Distance aux vignes
tmp=(rOCSOL_RO==2212)
tmp[tmp[]==0]=NA
rDIST_RO_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:
rMNT_RO_incrust=rMNT_RO-((rHYDRO_RO>0)*2)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"4.3 Illustration des données préparées
par(mfrow=c(1,2)) # Division de la fenêtre graphique
plot(rOCSOL_RO,col=as.character(Correspondance$Color),legend=F,axes=F)
plot(rHYDRO_RO,legend=F,axes=F)
par(mfcol=c(1,4))
plot(rMNT_RO,col=bpy.colors(),axes=F)
plot(rSLOPE_RO,col=bpy.colors(),axes=F)
plot(rASP_RO,col=bpy.colors(),axes=F)
plot(rHILL_RO,col=bpy.colors(),axes=F)
par(mfcol=c(1,3))
plot(rDIST_RO_exutoire,legend=F,axes=F,col=grey(100:1/100),main="Distance a l'exutoire")
plot(rDIST_RO_vignes,legend=F,axes=F,col=grey(100:1/100),main="Distance aux vignes")
plot(rDIST_RO_routes,legend=F,axes=F,col=grey(100:1/100),main="Distance aux routes")
4.4 Enregistrement des données afin d’être lues par le modèle
writeRaster(rOCSOL_RO,paste("datas/ocsol",sep=""),format="ascii",overwrite=T)
writeRaster(rMNT_RO_incrust,paste("datas/mnt",sep=""),format="ascii",overwrite=T)
writeRaster(rDIST_RO_exutoire,paste("datas/Distance_Exutoire",sep=""),format="ascii",overwrite=T)
writeRaster(rDIST_RO_vignes,paste("datas/Distance_Vignes",sep=""),format="ascii",overwrite=T)
writeRaster(rDIST_RO_routes,paste("datas/Distance_Routes",sep=""),format="ascii",overwrite=T)
writeRaster(rHILL_RO,paste("datas/hillShade",sep=""),format="ascii",overwrite=T)
writeRaster(rOCSOL_RO_ID,paste("datas/IDplots",sep=""),format="ascii",overwrite=T)
writeRaster(rHYDRO_RO_ID,paste("datas/IDfosses",sep=""),format="ascii",overwrite=T)5 Etapes du TP sous NETLOGO
5.1 Présentation de la plateforme
Logo et les tortues. Ce langage a été créé en 1967 par W. Feurzeig et S. Papert avec un objectif éducatif.
- simplicité du langage
- feedback visuel permettant le débuggage
- tortue robotisée (Irving)
NetLogo est un environnement de programmation pour la modélisation/simulation de phénomènes collectifs naturels. Il est bien adapté à la modélisation de systèmes complexes composés de centaines, de milliers d’agents agissant en parallèle. Il offre la possibilité de jouer avec de nombreuses simulations en sociologie, biologie, médecine, physique, chimie, mathématiques, informatique, économie et psychologie sociale, et également de créer ses propres modèles.
Le monde de Netlogo est constitué d’Agents, qui peuvent suivre des instructions. Les activités des différents agents s’exécutent simultanément. Il y a 3 types d’agents : tortues, patches et observateur.
Les Tortues sont les agents qui se déplacent dans le monde. Le Monde est en 2D ou 3D, divisé selon une grille (torique ou non) de patchs. Un Patch est une portion de sol sur laquelle les tortues peuvent se situer et se déplacer. L’Observateur regarde de l’extérieur le monde des tortues et des patchs (n’est pas situé dans le monde).
Interface graphique
Description
Onglet de code
5.2 Lancement du modèle de Ruissellement
Sous linux, ouvrir un terminal,
se placer dans le dossier du programme NetLogo via les commandes ls et cd
taper: sh netlogo.sh
Récupérer le modèle Hydrologie.nlogo dans le dossier datasSous windows, double cliquer sur le fichier Hydrologie.nlogoEnsuite, modifier le chemin d’accès vers les données (boutons HDD et set-directory) pour qu’il pointe sur le dossier de votre ordinateur.
5.3 Le modèle de ruissellement
Le modèle est organisé en deux modules, l’un figurant les flux d’eau et l’autre les agriculteurs.
- Chaque goutte d’eau est figurée par un agent, et dont le déplacement dépend de la gravité.
- L’agriculteur s’occupe uniquement des fossés, avec deux modes de pratiques: enherbement ou herbicides.
Diagramme 
5.3.1 Onglet de code
L’onglet de code est structuré en cinq parties:
Description des différents agents, patches et leurs variables associées
Initialisation du modèle, des agents et chargement des cartes
Ordre temporel des actions
Calcul des procédures et fonctions associées aux agents
Visualisation des cartes et des paramètres des agents
L’interface graphique se compose de deux types de boutons, de gauche à droite de la carte: Côté Gauche Les boutons permettent l’initialisation du modèle et la modification de certains paramètres d’importance. Ils vont conditionner les simulations effectuées. Côté Droit Les boutons servent à illustrer les valeurs de certaines variables du modèle, telles que l’occupation du sol, l’élevation, etc..
Interface graphique: Modèle de ruissellement
5.3.2 Pluie et ruissellement
Au début de la simulation, on simule de la pluie de manière aléatoire dans le temps et l’espace. Les gouttes de pluie se dispersent ensuite dans la zone par ruissellement superficiel.
Simulation de la pluie aléatoire
5.3.3 Agriculteurs et pratiques agricoles
Les agriculteurs sont répartis aléatoirement dans la zone, avec deux types de pratiques différentes.
Distribution des agriculteurs sur la zone
Chaque agriculteur applique sa pratique dans un rayon définissant sa zone d’action, pouvant être modifié. Les fossés en vert et en rouge correspondent à des pratiques d’enherbement et l’application de herbicides, respectivement.
Aperçu des deux types de pratiques agricoles dans les fossés
Au cours de la simulation, les agriculteurs appliquent leur pratique au bout de 30 jours de simulation. A partir de cette instant, les gouttes d’eau vont se charger ou se décharger en herbicides au fur et à mesure qu’elles traversent des zones polluées par les herbicides ou enherbées. Les graphiques sous la carte indiquent le niveau de pollution de l’eau, sur toute la zone (à gauche) ou uniquement à l’exutoire du bassin versant (à droite).
Niveau de pollution de l’eau
D’autres cartes ont été chargées dans le modèle, telles que l’élevation et la distance à une occupation du sol donnée.
Cartes d’élevation
Cartes de distance à une occupation du sol donnée
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?