TD 1 - Prise en main de PyQGIS

Dans ce TD nous allons nous concentrer sur l'interaction avec le logiciel QGIS en utilisant des commandes Python fournies par la librairie PyQGIS. L'interaction se fera directement depuis la console Python de QGIS qui possède l'environnement nécessaire pour utiliser cette librairie. L'objectif est de faire toutes les manipulations d'objets de QGIS seulement en commande Python sans utiliser l'interface graphique.

0 - Mise en place

Ouvrez QGIS puis créez un nouveau projet TD1.qgz. Enregistrez le projet dans le dossier de votre choix.

1 - Intéraction avec l'interface

1.1 Interface QGIS

Depuis la console Python vous avec accès à une instance de la classe QgisInterface via la variable iface qui pointe sur l'interface active de votre fenêtre QGIS. Cette classe contient des méthodes pour intéragir avec l'interface. Presque toutes les actions que vous pouvez faire en cliquant sur l'interface peuvent être réalisées à partir de la classe QgisInterface. Cette classe permet également la communication avec l'utilisateur via des méthodes pour afficher des informations et des erreurs.

1.1.1 Ajout de couches

1.1.1 Ajout de couches

Ajoutez la couche vectorielle RPG à votre projet en utilisant iface. Ensuite récupérez la couche RPG dans une variable Python.

Solution

iface.addVectorLayer("/path/to/RPG.gpkg", "RPG", "ogr")
# Séléctionnez la couche dans l'interface QGIS
rpg = iface.activeLayer()

1.1.2 Zoom

1.1.2 Zoom

Faites une séléction avec votre souris puis utilisez une action depuis QgisInterface pour zoomer sur la séléction.

Tip

Les méthodes comportant 'action' dans leur nom sont des QAction de la librairie PyQt5. L'action est déclenchée en appellant la méthode trigger().

Solution

iface.actionZoomToSelected().trigger()

Info

En POO (Programation Orientée Objet) il est courant que des méthodes de certaines classes renvoient d'autres classes, on appelle ça de la composition. C'est le cas de QgsMapCanvas depuis QgsProject.

1.1.3 Canvas

1.1.3 Canvas

Récupérez l'objet QgsMapCanvas depuis iface et affichez dans la console les éléments suivants:

• Les coordonnées du centre du canvas.
• L'étendue du canvas.
• La liste des couches visibles dans le canvas.
• Les dernières coordonnées de votre curseur de souris.

Solution

canvas = iface.mapCanvas()
canvas.center()
canvas.extent()
canvas.layers()
canvas.mouseLastXY()

1.1.4 Messages

1.1.4 Messages

Affichez un message puis une erreur avec la classe QgsMessageBar.

Tip

Il faut récupérer l'instance QgsMessageBar depuis l'interface active avec la classe iface.

Solution

message_bar = iface.messageBar()
message_bar.pushMessage("Hello world!")
message_bar.pushCritical("Layer error", "There is an error with layer XXXX.")

1.2 Projet QGIS

La classe QgsProject permet d'intéragir avec le projet courant dans QGIS.

1.2.1 Métadonnées

1.2.1 Métadonnées

Récupérez le nom, le chemin du projet courant et son sytème de coordonnées de réference (CRS). Ensuite donnez un titre au projet.

Solution

project = QgsProject.instance()
project.baseName()
project.fileName()
project.crs()
project.setTitle("TD1")

1.2.2 CRS

1.2.2 CRS

Changez le CRS du projet en le mettant sur 4326 puis rebasculez le sur le CRS d'origine.

Tip

La gestion des CRS en pyQGIS passe par la classe QgsCoordinateReferenceSystem. Pensez à récupérer l'id ou l'instance du CRS d'origine pour pouvoir rebasculez sur celui ci après la modification.

Solution

origin_crs = project.crs()
new_crs = QgsCoordinateReferenceSystem(4326)
project.setCrs(new_crs)
project.crs()
project.setCrs(origin_crs)

Info

Vous pouvez également ajouter une couche au projet depuis une variable Python via Qgsproject.instance() en utilisant addMapLayer()

2 - Données vectorielles

Le coeur de l'API Python pyQGIS réside dans la manipulation des couches de données et particulièrement pour les couches vectorielles. On peut créer et modifier les couches et les features avec plus de contrôle qu'en passant par les algorithmes processing.

2.1 QgsVectorLayer

La classe QgsvectorLayer est une classe qui regroupe un ensemble de méthodes pour agir sur une couche vectorielle et qui, en majorité, pointent vers des fonctions C++.

2.1.1 Chargez une couche

2.1.1 Chargez une couche

Chargez la couche vectorielle RPG dans une variable puis extrayez:

• Le nom la couche.
• Le CRS de la couche.
• Le chemin du fichier sur le disque.
• L'identifiant QGIS de la couche.

Puis changez le nom et le CRS de la couche.

Tip

La classe QgsVectorLayer hérite de la classe QgsMapLayer qui défini des méthodes génériques pour tout les types de couches dans QGIS (vecteur et raster).

Solution

rpg = QgsVectorLayer("/path/to/RPG.gpkg", "RPG", "ogr")
rpg.name()
rpg.crs()
rpg.source()
rpg.id()
origin_crs = rpg.crs()
rpg.setName("Parcelles agricoles")
rpg.setCrs(QgsCoordinateReferenceSystem(2154))
rpg.setCrs(origin_crs)

Info

Les champs des couches vectorielles sont gérées par les objets QgsField (un champ) et QgsFields (conteneur de champs). Les types de données sont gérés par la librairie PyQt5 via les classes QVariant et QMetaType ce qui permet une bonne interaction avec l'interface graphique. Notez que ce sont des classes d'énumerations et donc la valeur des types sont toujours des int. Essayez par exemple de rentrer QVariant.String dans la console.

2.1.2 Champs

2.1.2 Champs

Affichez dans la console tout les champs de la couche RPG avec leur noms et leur types.

Solution

fields = rpg.fields()
for f in fields:
print(f"name: {f.name()}, type: {f.type()}")

Info

Il y a deux manières de modifier les couches vectorielles. La première est de passer par un buffer d'édition (startEditing, commitChanges, rollBack). Cela permet de stocker en mémoire toutes les modifications et de les annuler si besoin ou si on rencontre une erreur lors d'un traitement. C'est la méthode sûre mais elle présente les inconvénients d'être plus lente et plus complexe à écrire que lorsque l'on applique les modifications directement via le fournisseur de données (dataProvider). En passant par le fournisseur de données les éditions sont plus rapides mais le retour en arrière n'est pas possible.

2.1.3 Modification de champs

2.1.3 Modification de champs

Ajoutez un champ AREA aux champs de la couche RPG puis mettez à jour les champs de la couche. Observez ensuite la table d'attributs.

Solution

# With edition buffer
rpg.startEditing()
rpg.addAttribute(QgsField("Area", 6))
rpg.commitChanges()
rpg.updateFields()

# edition buffer with context edit
with edit(rpg):
    rpg.addAttribute(QgsField("Area", 6))
    rpg.updateFields()

# with dataProvider
rpg.dataProvider().addAttributes([QgsField("Area", 6)])
rpg.updateFields()

2.1.4 Index de champs

2.1.4 Index de champs

Retrouvez la position du champ AREA et du champ CODE_CULTU à partir de leur noms.

Solution

area_id = rpg.fields().indexFromName("AREA")
code_cultu_id = rpg.fields().indexFromName("CODE_CULTU")

Info

De nombreuses manipulations des attributs des features de la couches se font avec la connaissance de la position du champ. Attention si la position vaut -1 cela signifie que le nom demandé n'est pas retrouvé dans la liste des champs. Il est fréquent que pyqgis ne remonte pas d'erreurs lorsque des actions impossibles sont demandées par l'utilisateur, c'est le cas ici.

2.1.5 Suppression de champs

2.1.5 Suppression de champs

Supprimez la champ AREA.

Solution

rpg.dataProvider().deleteAttributes([area_id])

2.2 QgsFeature

Les éléments de la couche sont appellés entités et représentés par la classe QgsFeature.

Info

Les requêtes de features depuis une couche vectorielle renvoient toujours des iterateurs. Transformez l'itrérateur en liste pour avoir toutes les features en mémoire. Attention il n'est pas recommandé de faire ça sur de grandes couches car cela peut devenir lent.

2.2.1 Extraire les features

2.2.1 Extraire les features

Extrayez le nombre et la liste des features de la couche RPG.

Solution

count = rpg.featureCount()
features = list(rpg.getFeatures())        

2.2.2 Extraire des features

2.2.2 Extraire des features

Récupérez la feature avec l'identifiant de feature 5. Puis récupérez en une commande les features 10, 200 et 1500.

Solution

feat = rpg.getFeature(5)
feat_10, feat_200, feat_1500 = tuple(rpg.getFeatures([10,20,1500]))

2.2.3 Métadonnées de feature

2.2.3 Métadonnées de feature

Depuis une des features récupérez:

• L'identifiant de la feature (fid).
• La géometrie de la feature.
• La liste des champs (idem que pour les couches).
• La liste des attributs (les valeurs pour chaque champs).
• La valeur pour le champ CODE_CULTU.

Solution

feature.id() 
feature.fields() 
feature.attributes() 
feature.geometry()
feature["CODE_CULTU]

2.2.4 Modification de feature

2.2.4 Modification de feature

Modifiez la feature en changeant la valeur pour le champ CODE_CULTU et appliquez un buffer de 5m sur sa géométrie. Ensuite remplacez la feature initiale par la nouvelle feature et zommez sur cette feature avec iface

Solution

feat = rpg.getFeature(5)
feat["CODE_CULTU"] = "YO"
geom = feat.geometry().buffer(5,5)
feat.setGeometry(geom)
feat_id = feat.id()
rpg.selectByIds([feat_id])

with edit(rpg):
    rpg.updateFeature(feat)

2.2.5 Ajout de feature

2.2.5 Ajout de feature

Créez une nouvelle feature avec les mêmes champs que ceux de la couche RPG puis ajoutez la à la couche RPG.

Tip

Le champ fid a une contrainte d'unicité. Il ne peut pas y avoir 2 features avec le même identifiant de feature donc pour ajoutez une feature à une couche pensez à mettre un fid > featureCount(). Le champ fid est un champ protégé.

Solution

feat = QgsFeature(rpg.fields())
feat.setAttributes([rpg.featureCount() + 1, "PPH", "156543514"])
rpg.dataProvider().addFeature(feat)

2.3 QgsGeometry

Les géométries en pyQGIS sont des instance des sous-classes de la classe QgsAbstractGeometry. Pour être associées avec des features de la couche elles sont contenues dans la classe QgsGeometry qui agit comme un wrapper. Dans la majorité des cas nous allons directement intéragir avec la classe QgsGeometry mais il est important de comprendre ce fonctionnement et les différentes sous classes pour avoir une manipulation fine des objets géométriques.

2.3.1 Métadonnées

2.3.1 Métadonnées

Récupérez la géometrie de feature 100 de la couche RPG et inspectez:

• Son type.
• Son type WKB.
• Son type WKT.

Ensuite créez une géométrie depuis la représentaation WKB de la géométrie de la feature 100.

Solution

geometry = feature.geometry()
geometry.type() # Retour le type qui correspond à une des sous classes de QgsAbstractGeometry
geometry.asWkb()
geometry.asWkt()

# Création
new_geom = QgsGeometry()
new_geom.fromWkb(geometry.asWkb())

2.3.2 Manipulation de géométries

2.3.2 Manipulation de géométries

Transformez la géometrie en un polygone avec uniquement un point sur deux.

Tip

Utilisez la méthode asPolygon pour récupérer une liste de points et la classe QgsLineString pour créer un polygone à partir de points. Assurez vous que le polygone n'est pas en plusieurs parties au niveau de la geometry avec isMultipart() ou au niveau des points (liste de plus d'une liste de points → multi-parties). Si le polygone est multi-parties prenez uniquement la première partie.

Solution

geom = feature.geometry()
geom = geom.constGet().geometryN(0)
geom = QgsGeometry(geom)
points = geom.asPolygon()[0]
new_points = points[::2]
linestring = QgsLineString(new_points)
new_poly = QgsPolygon(linestring)
new_geom = QgsGeometry(new_poly)

2.3.2 Géométries dérivées

2.3.2 Géométries dérivées

Dessinez la plus petite boîte englobante orientée autour du nouveau polygone.

Solution

bbox = geom.orientedMinimumBoundingBox()

2.4 Requêtes

On peut faire des requêtes sur les couches vectorielles afin de récupérer les features qui respectent certaines conditions, spatiales ou autres. Ces requêtes se font via la classe QgsFeatureRequest.

Info

La méthode getFeatures des couches vectorielles prend comme paramètre une QgsFeatureRequest qui par défaut est vide mais si vous passez comme paramètre une requête avec une condition la méthode renverra uniquement les features qui respectent la condition. Vous pouvez également créez une couche à partir d'une requête et d'une couche d'origine en utilisant la méthode materialize().

2.4.1 Requête spatiale

2.4.1 Requête spatiale

Récupérez une feature depuis la couche RPG. Appliquez un buffer de 5km à la géométrie. Extrayez toutes les features de la couche RPG qui intersectent la géométrie avec buffer et calculez l'aire totale des géométries qui la superposent. Calculez le temps d'execution de la requête spatiale en écrivant le code dans un script à côté de la console.

Tip

Les requêtes spatiales se font souvent depuis des boîtes englobantes, vous devez ensuite découper les géométries de chaque feature requêtée pour ne garder que la partie qui intersecte la vraie géométrie (et non la boîte englobante) avec la méthode intersection de QgsGeometry.

Solution

import time

feat = rpg.getFeature(537)
geom_buffer = feat.geometry().buffer(5000,5)
t1 = time.time()
query = QgsFeatureRequest().setFilterRect(geom_buffer.boundingBox())
features = rpg.getFeatures(query)
t2 = time.time()
print(t2-t1)
geometries = [f.geometry() for f in features]
geometries = [g.intersection(geom_buffer).area() for g in geometries]
area_tot = sum(geometries)

2.4.2 Index spatial

2.4.2 Index spatial

Faites la même chose qu'en 2.4.1 en utilisant un index spatial QgsSpatialIndex. Comparez le temps de requête spatiale (sans le coût de création de l'index).

Solution

import time

feat = rpg.getFeature(537)
geom_buffer = feat.geometry().buffer(5000,5)
spatial_index = QgsSpatialIndex(rpg.getFeatures())
t1 = time.time()
feats_ids = spatial_index.intersects(geom_buffer.boundingBox())
features = rpg.getFeatures(feats_ids)
t2 = time.time()
print(t2-t1)
geometries = [f.geometry() for f in features]
geometries = [g.intersection(geom_buffer).area() for g in geometries]
area_tot = sum(geometries)

Info

La création d'index spatial présente un coup en temps mais son utilisation accélère les requêtes spatiales lorsque le nombre de features est élevé.

2.4.3 Requête sur les attributs

2.4.3 Requête sur les attributs

Créez une couche qui ne contient que les features dont le champ CODE_CULTU vaut PPH.

Solution

query = QgsFeatureRequest().setFilterExpression("CODE_CULTU = 'PPH'")
rpg_pph = rpg.materialize(query)

2.4.4 Requête par séléction

2.4.4 Requête par séléction

Avec votre curseur sélectionnez plusieurs features depuis l'interface. Ensuite, en ligne de commande, créez une couche qui ne contient que les features séléctionneées.

Tip

Utilisez les id des features séléctionnées pour construire votre requête.

Solution

query = QgsFeatureRequest().setFilterFids(rpg.selectedFeatureIds())
rpg_selected = rpg.materialize(query)

2.5 Expressions

Les expressions sont un bon vecteur de communication entre un utilisateur et un algortihme. Nous avons vu que via les expressions nous pouvons effectuer des requêtes pour filtrer des couches vectorielles. Elles peuvent également servir à faire des modifications sur les couches vectorielles. On utilise les expressions via la classe QgsExpression.

2.5.1 Modifications par expression

2.5.1 Modifications par expression

Ajoutez les champs AREA (Double), FACTOR (Int) et RESULT (Double) à la couche RPG. Uniquement avec des expressions remplissez pour chaque champ respectivement:

• L'aire du polygone.
• Un entier correspondant à un mapping entre une catégorie CODE_CULTU et un entier de 1 à N +1. N étant le nombre de catégorie CODE_CULTU. Le mapping est à créer.
• Le résultat de AREA / FACTOR.

Tip

Utilisez la méthode uniqueValues() de QgsVectorLayer pour avoir le set de catégories de CODE_CULTU. Utilisez QgsExpressionContext pour appliquer vos opérations et QgsExpressionContextScope pour ajouter des variables externes dans vos expressions. Quand vous assignez une valeur à un champ d'une feature il faut réassigner la feature au QgsExpressionContext sinon la valeur utilisée lors de l'évaluation de l'expression sera l'ancienne valeur. Cette question est difficile, n'hésitez pas à demander de l'aide.

Solution

# création du mapping
code_cultu_uniques = rpg.uniqueValues(rpg.fields().indexFromName("CODE_CULTU"))
factor_map = dict(zip(code_cultu_uniques, range(1, len(code_cultu_uniques))))
# Ajout de variable au contexte
expr_context = QgsExpressionContext()
scope = QgsExpressionContextScope()
scope.setVariable("factor_map", factor_map)
expr_context.appendScope(scope)
# Définitions des 3 expressions
expr_area =  QgsExpression("area($geometry)")
expr_factor = QgsExpression("map_get(@factor_map, CODE_CULTU)")
expr_result = QgsExpression("AREA / FACTOR")
# Application des expressions
rpg.startEditing()

for feat in rpg.getFeatures():
    expr_context.setFeature(feat)
    feat["AREA"] =  expr_area.evaluate(expr_context)
    feat["FACTOR"] = expr_factor.evaluate(expr_context)
    expr_context.setFeature(feat)
    feat["RESULT"] = expr_result.evaluate(expr_context)
    rpg.updateFeature(feat)

rpg.commitChanges()

2.6 Création et écriture de couches

2.6.1 Création et sauvegarde de couche

2.2.6 Création et sauvegarde de couche

Créez une nouvelle couche de points en mémoire avec le CRS du projet. Remplissez la couche avec le centroïde des polygones de RPG. Enregistrez la couche avec la classe QgsVectorFileWriter au format GPKG et SHP.

Tip

Utilisez l'URI pour définir le type de géométrie et le CRS.

Solution

project_crs = QgsProject.instance().crs()
points = QgsVectorLayer(f"Point?crs={project_crs.authid()}", "points", "memory")
features = []
for feat in rpg.getFeatures():
    point_feat = QgsFeature()
    point_feat.setGeometry(feat.geometry().centroid())
    features.append(point_feat)

points.dataProvider().addFeatures(features)
QgsProject.instance().addMapLayer(points)

# GPKG
save_options = QgsVectorFileWriter.SaveVectorOptions()
transform_context = QgsProject.instance().transformContext()
error = QgsVectorFileWriter.writeAsVectorFormatV3(layer,"/path/to/points.gpkg",transform_context,save_options)
# SHP
save_options.driverName = "ESRI Shapefile"
error = QgsVectorFileWriter.writeAsVectorFormatV3(layer,"/path/to/points.shp",transform_context,save_options)

3 - Données raster

On peut également manipulez des couches raster avec pyQGIS. Cependant les possibilités sont plus réduites que pour les couches vectorielles et la syntaxe de code plus complexe. Dans la majorité des cas on utilisera des algorithmes processing pour traiter des données raster.

3.1 QgsRasterLayer

3.1.1 Chargement de couche

3.1.1 Chargement de couche

Chargez la couche raster indices.tiff et ajoutez là au projet.

Solution

indices = QgsRasterLayer("/path/to/indices.tiff")

3.1.2 Métadonnées

3.1.2 Métadonnées

Extrayez les informations suivantes:

• Nombre de bandes.
• Nom des bandes.
• Hauteur et largeur du raster.
• Unités par pixel de hauteur/largeur.
• Étendue du raster.

Solution

raster = QgsRasterLayer("/path/to/raster.tiff")
raster.bandCount()
raster.bandName(x)
raster.height()
raster.width()
rasterUnitsPerPixelY()
rasterUnitsPerPixelX()
raster.extent() # Peut aussi être obtenu en multipliant height/width par rasterUnitsPerPixelX/Y.

3.1.3 Valeurs de pixels

3.1.3 Valeurs de pixels

Récupérez la valeur de chaque bande depuis un point QgsPointXY dans le raster.

Tip

Vous pouvez récupérer les coordonnées en QgsPointXY du centre du canvas avec la méthode center() de QgsMapCanvas.

Solution

p = iface.mapCanvas().center()
value, _ = indices.dataProvider().sample(p, 1)

3.2 Numpy et QgsRasterLayer

Pour des analyses plus fine des données du raster on peut transformer la couche raster en array numpy. Cependant une array numpy est chargée en mémoire et peu donc la surcharger si elle trop grande.

3.2.1 Transformation en array

3.2.1 Transformation en array

Extrayez l'array numpy qui correspond à la couche raster et calculez la médiane de chaque bandes.

Solution

import numpy as np
array = indices.as_numpy()
median = np.nanmedian(array, axis=(1,2))

4 - Algorithmes processing

Il est possible d'utiliser des algorithmes de la boîte à outils via pyQGIS grâce au module processing et la méthode run().

Info

Plusieurs options existent pour les sorties des algorithmes (souvent appellées OUTPUT):

1. memory: pour les sorties en mémoire.
2. TEMPORARY_OUTPUT pour les fichiers temporaires (stockés dans le dossier temp de l'ordinateur et supprimés lorsqu'on l'éteint.)
3. Un chemin vers un nouveau fichier pour écrire la couche sur le disque.

Lorsqu'on travaille dans la console Python de QGIS memory: et TEMPORARY_OUTPUT renvoient des couches en mémoire car QGIS n'a pas besoin d'écrire le fichier pour le retrouver. Quand on travaille sans l'interface alors TEMPORARY_OUTPUT créera un fichier physique sur le disque.

4.1.1 Utilisez processing

4.1.1 Utilisez processing

Calculez la moyenne de NDVI par polygone du RPG en utilisant l'algorithme "Statistiques de zones" via processing. Ajoutez la couche de résultats au projet QGIS.

Solution

zonal_stat_outputs = processing.run("native:zonalstatisticsfb", {'INPUT':rpg,'INPUT_RASTER':indices,'RASTER_BAND':1,'COLUMN_PREFIX':'NDVI_','STATISTICS':[2],'OUTPUT':'TEMPORARY_OUTPUT'})
zonal_stat = zonal_stat_outputs["OUTPUT"]
QgsProject.instance().addMapLayer(zonal_stat)

5 - PyQGIS sans l'interface QGIS

Si on lance des programmes trop lourds depuis l'interface de QGIS (console ou script) l'application QGIS se fige. C'est dû au fait que l'interface ainsi que les scripts lancés sont dans la même thread. Quand QGIS traite les informations de votre programme il ne peut plus gérer les actions que vous demandez de lui faire depuis l'interface. Pour utiliser l'API Python de QGIS sans créer ce problème il faut lancer le programme entièrement en dehors de QGIS. Mais python ne trouvera pas la librairie pyqgis si l'application n'est pas initialisée. Il faut donc initialisée l'application sans interface graphique.

5.1.1 Initialiser QGIS

5.1.1 Initialiser QGIS

Fermez QGIS et ouvrez un script Python avec l'éditeur de votre choix et initialisez l'application QGIS. Chargez les couches indices.tiff et RPG dans votre script réalisez la même action que dans la question 4.1.1.

Tip

• La librairie pyqgis est importée en avec import qgis et la majorité des classes à utilisées sont dans qgis.core.
• La classe QgsApplication vous permet d'initialiser QGIS.
• Pour trouver le chemin de l'installation de QGIS. Ouvrez l'application et dans la console Python entrez QgsApplication.prefixPath().
• Vous pouvez importer processing directement avec import processing mais il faut l'initialiser processing.core.Processing.Processing.initialize() pour qu'il charge les algortihmes.

Solution

from qgis.core import QgsApplication, QgsVectorLayer, QgsRasterLayer
import processing

QgsApplication.setPrefixPath('/installation/path', True)
qgs = QgsApplication([], False)
qgs.initQgis()

processing.core.Processing.Processing.initialize()

rpg = QgsVectorLayer("/home/lucas/Documents/cours/cours-pyqgis/dev/data/rpg.gpkg")
indices = QgsRasterLayer("/home/lucas/Documents/cours/cours-pyqgis/dev/data/indices.tiff")

zonal_stat_outputs = processing.run("native:zonalstatisticsfb", {'INPUT':rpg,'INPUT_RASTER':indices,'RASTER_BAND':1,'COLUMN_PREFIX':'NDVI_','STATISTICS':[2],'OUTPUT':'TEMPORARY_OUTPUT'})
zonal_stat = zonal_stat_outputs["OUTPUT"]

print(zonal_stat)

qgs.exitQgis()

6 - Exercices

6.1 Paysage et biodiversité

6.1 Paysage et biodiversité

Dans le cadre de l'analyse du paysage on essaie de déterminer le potentiel de diversité des paysages. Ce potentiel peut être liée aux composantes du paysage et à sa mosaïque agricole. L'objectif dans cet exercice est d'extraire des métriques du paysages sur des étendues de \(1km^2\). Les métriques attendues sont:

• La richesse spécifique des cultures. Donc le nombre de cultures différentes.
• L'inverse de l'indice de biodiversité simpson calculé avec la formule \(\frac{1}{\sum{{p_i}^2}}\) avec \(p_i\) la proportion de la culture i.
• La moyenne de la taille des parcelles agricoles.
• La surface de milieux semi-naturel \(1 - \frac{A_{artificiel} + A_{parcelles}}{A_{total}}\) avec \(A_{artificiel}\) la surface de structures artificielles et A_{parcelles} la surface de cultures totales.

Les couches à utilisées sont les couches RPG et artificial.

Le résultat attendu est une grille (Couche vectorielle) avec des carrés de \(1km^2\) qui recouvrent la zone d'étude et dans laquelle il y a pour chaque carré les métriques mentionnées.