TD1 6.1¶

Pour résoudre ce problème il faut bien segmenter les différentes tâches et identifier les actions redondantes pour les réaliser dans des fonctions qu'on utilisera ensuite dans le code final. Il existe plusieurs structures possibles dont la suivante:

1. Une fonction pour extraire les features qui sont dans la cellule.
2. Des fonctions pour calculer chacune des métriques:
   • 1 fonction pour la taille moyenne des champs.
   • 1 fonction pour la richesse et la diversité culturale.
   • 1 fonction pour la surface en milieux semi-naturels.
3. Un code final qui utilisera les fonctions définies précedemment.

0. Imports et initalisation de QgsApplication¶

Si vous travaillez en dehors de l'application QGIS il faut initaliser l'application sans interface et importer les classes nécessaires. Si vous travaillez dans l'application QGIS vous pouvez passer à la section suivante.

In [ ]:

from qgis.core import (
    QgsVectorLayer,
    QgsFeatureRequest,
    QgsRectangle,
    QgsFeature,
    QgsField,
    QgsCoordinateReferenceSystem,
    QgsApplication,
    QgsVectorFileWriter,
    QgsCoordinateTransformContext,
)
from PyQt5.QtCore import QVariant
from typing import List, Tuple
from statistics import mean

import processing

QgsApplication.setPrefixPath("/home/lucas/anaconda3/envs/cours_pyqgis", True)
qgs = QgsApplication([], False)
qgs.initQgis()

processing.core.Processing.Processing.initialize()

from qgis.core import ( QgsVectorLayer, QgsFeatureRequest, QgsRectangle, QgsFeature, QgsField, QgsCoordinateReferenceSystem, QgsApplication, QgsVectorFileWriter, QgsCoordinateTransformContext, ) from PyQt5.QtCore import QVariant from typing import List, Tuple from statistics import mean import processing QgsApplication.setPrefixPath("/home/lucas/anaconda3/envs/cours_pyqgis", True) qgs = QgsApplication([], False) qgs.initQgis() processing.core.Processing.Processing.initialize()

1 - Extraction de features¶

Il faut récupérer toutes les features qui intersectent la cellule cible. On utilise une requête spatial pour cela avec setFilterRect. En retour de fonction on veut la liste des features correspondantes.

In [2]:

def spatial_query(layer: QgsVectorLayer, cell: QgsRectangle) -> List[QgsFeature]:
    """Query features from layer that intersects cell.

    Args:
        layer (QgsVectorLayer): Layer that contains wanted features.
        cell (QgsRectangle): Rectangle to query features.

    Returns:
        List[QgsFeature]: List of features from layer that intersects cell.
    """
    query = QgsFeatureRequest().setFilterRect(cell)
    return list(layer.getFeatures(query))

def spatial_query(layer: QgsVectorLayer, cell: QgsRectangle) -> List[QgsFeature]: """Query features from layer that intersects cell. Args: layer (QgsVectorLayer): Layer that contains wanted features. cell (QgsRectangle): Rectangle to query features. Returns: List[QgsFeature]: List of features from layer that intersects cell. """ query = QgsFeatureRequest().setFilterRect(cell) return list(layer.getFeatures(query))

2 - Taille moyenne des champs¶

Pour avoir la moyenne des surfaces de parcelles il faut récupérer toutes les géométries des features. Récupérer uniquement la surface qui recouvre la cellule et calculer la moyenne des aires.

In [3]:

def mean_field_size(features: List[QgsFeature], cell: QgsRectangle) -> float:
    """Compute mean crop plots area.

    Args:
        features (List[QgsFeature]): List of crop features.
        cell (QgsRectangle):  Area extent to compute field size.

    Returns:
        float: Mean field size
    """
    return mean([f.geometry().intersection(cell).area() for f in features])

def mean_field_size(features: List[QgsFeature], cell: QgsRectangle) -> float: """Compute mean crop plots area. Args: features (List[QgsFeature]): List of crop features. cell (QgsRectangle): Area extent to compute field size. Returns: float: Mean field size """ return mean([f.geometry().intersection(cell).area() for f in features])

3 - Surface des milieux semi-naturels¶

Pour calculer la la surface des milieux semi-naturels on calcule la surface totale de la cellule et on soustrait:

• L'aire totale des parcelles.
• L'aire totale des surfaces artificielles.

In [4]:

def semi_natural_cover(
    cell: QgsRectangle,
    crop_features: List[QgsFeature],
    artificial_features: List[QgsFeature],
)-> float:
    """Compute semi natural cover from crops & artificial features.

    Args:
        cell (QgsRectangle): Area extent semi natural cover.
        crop_features (List[QgsFeature]): Crop features
        artificial_features (List[QgsFeature]): Artificial features.

    Returns:
        float: Semi natural cover.
    """
    cell_area = cell.area()
    crops_area = sum([f.geometry().intersection(cell).area() for f in crop_features])
    artificial_area = sum(
        [f.geometry().intersection(cell).area() for f in artificial_features]
    )

    return (cell_area - crops_area - artificial_area) / cell_area

def semi_natural_cover( cell: QgsRectangle, crop_features: List[QgsFeature], artificial_features: List[QgsFeature], )-> float: """Compute semi natural cover from crops & artificial features. Args: cell (QgsRectangle): Area extent semi natural cover. crop_features (List[QgsFeature]): Crop features artificial_features (List[QgsFeature]): Artificial features. Returns: float: Semi natural cover. """ cell_area = cell.area() crops_area = sum([f.geometry().intersection(cell).area() for f in crop_features]) artificial_area = sum( [f.geometry().intersection(cell).area() for f in artificial_features] ) return (cell_area - crops_area - artificial_area) / cell_area

4 - Diversité culturales¶

On calcule en même temps la richesse spécifique et l'inverse de l'index de simpson. Pour cela il faut la surface totale des parcelles et la surface totale de chaque culture.

In [5]:

def diversity(
    features: List[QgsFeature], field: str, cell: QgsRectangle
) -> Tuple[float, float]:
    """Compute richness and inverse simpson index from crop features.

    Args:
        features (List[QgsFeature]): Crop features.
        field (str): Field to extract crop type ("CODE_CULTU").
        cell (QgsRectangle): Area extent to compute diversity metrics.

    Returns:
        Tuple[float, float]: richness & inverse simpson.
    """
    crops = set([f[field] for f in features])
    richness = len(set([f[field] for f in features]))
    total_crop_area = sum([f.geometry().intersection(cell).area() for f in features])
    crops_area = dict(zip(crops, [0] * len(crops)))
    for feat in features:
        crops_area[feat[field]] += feat.geometry().intersection(cell).area()

    inv_simpson = 1 / sum(
        [(area / total_crop_area) ** 2 for area in list(crops_area.values())]
    )
    return richness, inv_simpson

def diversity( features: List[QgsFeature], field: str, cell: QgsRectangle ) -> Tuple[float, float]: """Compute richness and inverse simpson index from crop features. Args: features (List[QgsFeature]): Crop features. field (str): Field to extract crop type ("CODE_CULTU"). cell (QgsRectangle): Area extent to compute diversity metrics. Returns: Tuple[float, float]: richness & inverse simpson. """ crops = set([f[field] for f in features]) richness = len(set([f[field] for f in features])) total_crop_area = sum([f.geometry().intersection(cell).area() for f in features]) crops_area = dict(zip(crops, [0] * len(crops))) for feat in features: crops_area[feat[field]] += feat.geometry().intersection(cell).area() inv_simpson = 1 / sum( [(area / total_crop_area) ** 2 for area in list(crops_area.values())] ) return richness, inv_simpson

5 - Code final¶

Maintenant on peut tout regrouper pour calculer les métriques sur chaque cellule et regrouper ces résultats dans une nouvelle couche vectorielle. Par simplicité on ignore les cellules qui ne contiennent aucune parcelle agricole.

In [ ]:

#### Paramètres

rpg_file = "/home/lucas/Documents/cours/cours-pyqgis/dev/data/rpg.gpkg"
artificial_file = "/home/lucas/Documents/cours/cours-pyqgis/dev/data/aritificial.gpkg"

#### Algorithme

rpg = QgsVectorLayer(rpg_file)
artificial = QgsVectorLayer(artificial_file)

# Création de la grille avec un algorithme processing
grid: QgsVectorLayer = processing.run(
    "native:creategrid",
    {
        "TYPE": 2,
        "EXTENT": rpg,
        "HSPACING": 1000,
        "VSPACING": 1000,
        "HOVERLAY": 0,
        "VOVERLAY": 0,
        "CRS": QgsCoordinateReferenceSystem("IGNF:LAMB93"),
        "OUTPUT": "TEMPORARY_OUTPUT",
    },
)["OUTPUT"]

# Création de la couche de résultats
results_fields = [
    QgsField("RICHNESS", QVariant.Double),
    QgsField("DIVERSITY", QVariant.Double),
    QgsField("FIELD_SIZE", QVariant.Double),
    QgsField("SEMI_NATURAL", QVariant.Double),
]

result_layer = QgsVectorLayer(
    f"Polygon?crs={grid.crs().authid()}", "Results", "memory"
)  # en mémoire

result_layer.dataProvider().addAttributes(results_fields)
result_layer.updateFields()

# Calcul des métriques et enregistrement des features de la couche de résultats

new_features = [] # conteneur de features résultats

## itération sur les cellules de la grille
for cell_feat in grid.getFeatures():
    # requête spatiale des features qui intersectent la cellule
    cell = cell_feat.geometry()
    crops_features = spatial_query(rpg, cell.boundingBox())
    artificial_features = spatial_query(artificial, cell.boundingBox())
    # Cellule ignorée si aucune parcelle agricole
    if len(crops_features) == 0:
        continue

    # calcul des différentes métriques
    richness, inv_simpson = diversity(crops_features, "CODE_CULTU", cell)
    mfs = mean_field_size(crops_features, cell)
    snc = semi_natural_cover(cell, crops_features, artificial_features)

    # Création de la feature résultat

    feature = QgsFeature(result_layer.fields())
    feature.setGeometry(cell)

    feature["RICHNESS"] = richness
    feature["DIVERSITY"] = inv_simpson
    feature["FIELD_SIZE"] = mfs
    feature["SEMI_NATURAL"] = snc
    
    new_features.append(feature)

# Ajout de toutes les features résultats
result_layer.dataProvider().addFeatures(new_features)

#### Paramètres rpg_file = "/home/lucas/Documents/cours/cours-pyqgis/dev/data/rpg.gpkg" artificial_file = "/home/lucas/Documents/cours/cours-pyqgis/dev/data/aritificial.gpkg" #### Algorithme rpg = QgsVectorLayer(rpg_file) artificial = QgsVectorLayer(artificial_file) # Création de la grille avec un algorithme processing grid: QgsVectorLayer = processing.run( "native:creategrid", { "TYPE": 2, "EXTENT": rpg, "HSPACING": 1000, "VSPACING": 1000, "HOVERLAY": 0, "VOVERLAY": 0, "CRS": QgsCoordinateReferenceSystem("IGNF:LAMB93"), "OUTPUT": "TEMPORARY_OUTPUT", }, )["OUTPUT"] # Création de la couche de résultats results_fields = [ QgsField("RICHNESS", QVariant.Double), QgsField("DIVERSITY", QVariant.Double), QgsField("FIELD_SIZE", QVariant.Double), QgsField("SEMI_NATURAL", QVariant.Double), ] result_layer = QgsVectorLayer( f"Polygon?crs={grid.crs().authid()}", "Results", "memory" ) # en mémoire result_layer.dataProvider().addAttributes(results_fields) result_layer.updateFields() # Calcul des métriques et enregistrement des features de la couche de résultats new_features = [] # conteneur de features résultats ## itération sur les cellules de la grille for cell_feat in grid.getFeatures(): # requête spatiale des features qui intersectent la cellule cell = cell_feat.geometry() crops_features = spatial_query(rpg, cell.boundingBox()) artificial_features = spatial_query(artificial, cell.boundingBox()) # Cellule ignorée si aucune parcelle agricole if len(crops_features) == 0: continue # calcul des différentes métriques richness, inv_simpson = diversity(crops_features, "CODE_CULTU", cell) mfs = mean_field_size(crops_features, cell) snc = semi_natural_cover(cell, crops_features, artificial_features) # Création de la feature résultat feature = QgsFeature(result_layer.fields()) feature.setGeometry(cell) feature["RICHNESS"] = richness feature["DIVERSITY"] = inv_simpson feature["FIELD_SIZE"] = mfs feature["SEMI_NATURAL"] = snc new_features.append(feature) # Ajout de toutes les features résultats result_layer.dataProvider().addFeatures(new_features)

6 - Sauvegarde la couche¶

On sauvegarde la couche dans un fichier avec QgsVectorFileWriter si vous travaillez hors QGIS. Sinon vous pouvez juste ajoutez la couche au projet depuis la console.

In [7]:

save_options = QgsVectorFileWriter.SaveVectorOptions()
transform_context = QgsCoordinateTransformContext() 
error = QgsVectorFileWriter.writeAsVectorFormatV3(result_layer,"/home/lucas/Documents/cours/cours-pyqgis/dev/results.gpkg",transform_context,save_options)

save_options = QgsVectorFileWriter.SaveVectorOptions() transform_context = QgsCoordinateTransformContext() error = QgsVectorFileWriter.writeAsVectorFormatV3(result_layer,"/home/lucas/Documents/cours/cours-pyqgis/dev/results.gpkg",transform_context,save_options)