TD2 2.2.1

In [ ]:

from statistics import mean
from qgis.core import (
    QgsProcessingAlgorithm,
    QgsProcessingParameterFeatureSource,
    QgsProcessingParameterFeatureSink,
    QgsVectorLayer,
    QgsFeatureSink,
    QgsCoordinateReferenceSystem,
    QgsField,
    QgsRectangle,
    QgsFeatureRequest,
    QgsFeature,
    QgsProcessingMultiStepFeedback,
    QgsProcessingParameterNumber,
)
import processing
from PyQt5.QtCore import QVariant

from typing import List, Tuple


class LandscapeBiodiversityMetrics(QgsProcessingAlgorithm):

    RPG = "RPG"
    ARTIFICIAL = "ARTIFICIAL"
    METRICS = "METRICS"
    CELL_SIZE = "CELL_SIZE"

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFeatureSource(self.RPG, self.RPG))
        self.addParameter(
            QgsProcessingParameterFeatureSource(self.ARTIFICIAL, self.ARTIFICIAL)
        )

        self.addParameter(QgsProcessingParameterFeatureSink(self.METRICS, self.METRICS))
        self.addParameter(QgsProcessingParameterNumber(self.CELL_SIZE, self.CELL_SIZE))

    def processAlgorithm(self, parameters, context, feedback):
        rpg: QgsVectorLayer = self.parameterAsVectorLayer(parameters, self.RPG, context)
        artificial = self.parameterAsVectorLayer(parameters, self.ARTIFICIAL, context)
        cell_size = self.parameterAsDouble(parameters, self.CELL_SIZE, context)
        multi_feedback = QgsProcessingMultiStepFeedback(2, feedback=feedback)
        multi_feedback.setCurrentStep(0)
        grid: QgsVectorLayer = processing.run(
            "native:creategrid",
            {
                "TYPE": 2,
                "EXTENT": rpg,
                "HSPACING": cell_size,
                "VSPACING": cell_size,
                "HOVERLAY": 0,
                "VOVERLAY": 0,
                "CRS": QgsCoordinateReferenceSystem("IGNF:LAMB93"),
                "OUTPUT": "TEMPORARY_OUTPUT",
            },
            feedback=multi_feedback,
            context=context,
        )["OUTPUT"]

        multi_feedback.setProgress(1)

        new_fields = [
            QgsField("RICHNESS", QVariant.Double),
            QgsField("DIVERSITY", QVariant.Double),
            QgsField("FIELD_SIZE", QVariant.Double),
            QgsField("SEMI_NATURAL", QVariant.Double),
        ]

        metrics = QgsVectorLayer(
            f"Polygon?crs={grid.crs().authid()}", "Results", "memory"
        )

        metrics.dataProvider().addAttributes(new_fields)
        metrics.updateFields()
        new_features = []
        multi_feedback.setCurrentStep(1)
        total_feats = grid.featureCount()
        for i, cell_feat in enumerate(grid.getFeatures()):
            cell = cell_feat.geometry()
            crops_features = self.spatial_query(rpg, cell.boundingBox())
            artificial_features = self.spatial_query(artificial, cell.boundingBox())

            if len(crops_features) == 0:
                continue

            richness, inv_simpson = self.diversity(crops_features, "CODE_CULTU", cell)

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

            feature["RICHNESS"] = richness
            feature["DIVERSITY"] = inv_simpson
            feature["FIELD_SIZE"] = self.mean_field_size(crops_features, cell)
            feature["SEMI_NATURAL"] = self.semi_natural_cover(
                cell, crops_features, artificial_features
            )
            new_features.append(feature)
            multi_feedback.setProgress(int((i / total_feats) * 100))

        metrics.dataProvider().addFeatures(new_features)

        (sink, dest_id) = self.parameterAsSink(
            parameters,
            self.METRICS,
            context,
            metrics.fields(),
            metrics.wkbType(),
            metrics.sourceCrs(),
        )

        for f in metrics.getFeatures():
            sink.addFeature(f, QgsFeatureSink.FastInsert)

        return {self.METRICS: dest_id}

    def name(self):
        return "landscapebiodiv"

    def displayName(self):
        return "Landscape Biodiversity"

    def group(self):
        return "Custom processing"

    def groupId(self):
        return "custom_processing"

    def createInstance(self):
        return LandscapeBiodiversityMetrics()

    def diversity(
        self, 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.
            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 mean_field_size(self, 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 semi_natural_cover(
        self,
        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 spatial_query(
        self, 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))

from statistics import mean from qgis.core import ( QgsProcessingAlgorithm, QgsProcessingParameterFeatureSource, QgsProcessingParameterFeatureSink, QgsVectorLayer, QgsFeatureSink, QgsCoordinateReferenceSystem, QgsField, QgsRectangle, QgsFeatureRequest, QgsFeature, QgsProcessingMultiStepFeedback, QgsProcessingParameterNumber, ) import processing from PyQt5.QtCore import QVariant from typing import List, Tuple class LandscapeBiodiversityMetrics(QgsProcessingAlgorithm): RPG = "RPG" ARTIFICIAL = "ARTIFICIAL" METRICS = "METRICS" CELL_SIZE = "CELL_SIZE" def initAlgorithm(self, config=None): self.addParameter(QgsProcessingParameterFeatureSource(self.RPG, self.RPG)) self.addParameter( QgsProcessingParameterFeatureSource(self.ARTIFICIAL, self.ARTIFICIAL) ) self.addParameter(QgsProcessingParameterFeatureSink(self.METRICS, self.METRICS)) self.addParameter(QgsProcessingParameterNumber(self.CELL_SIZE, self.CELL_SIZE)) def processAlgorithm(self, parameters, context, feedback): rpg: QgsVectorLayer = self.parameterAsVectorLayer(parameters, self.RPG, context) artificial = self.parameterAsVectorLayer(parameters, self.ARTIFICIAL, context) cell_size = self.parameterAsDouble(parameters, self.CELL_SIZE, context) multi_feedback = QgsProcessingMultiStepFeedback(2, feedback=feedback) multi_feedback.setCurrentStep(0) grid: QgsVectorLayer = processing.run( "native:creategrid", { "TYPE": 2, "EXTENT": rpg, "HSPACING": cell_size, "VSPACING": cell_size, "HOVERLAY": 0, "VOVERLAY": 0, "CRS": QgsCoordinateReferenceSystem("IGNF:LAMB93"), "OUTPUT": "TEMPORARY_OUTPUT", }, feedback=multi_feedback, context=context, )["OUTPUT"] multi_feedback.setProgress(1) new_fields = [ QgsField("RICHNESS", QVariant.Double), QgsField("DIVERSITY", QVariant.Double), QgsField("FIELD_SIZE", QVariant.Double), QgsField("SEMI_NATURAL", QVariant.Double), ] metrics = QgsVectorLayer( f"Polygon?crs={grid.crs().authid()}", "Results", "memory" ) metrics.dataProvider().addAttributes(new_fields) metrics.updateFields() new_features = [] multi_feedback.setCurrentStep(1) total_feats = grid.featureCount() for i, cell_feat in enumerate(grid.getFeatures()): cell = cell_feat.geometry() crops_features = self.spatial_query(rpg, cell.boundingBox()) artificial_features = self.spatial_query(artificial, cell.boundingBox()) if len(crops_features) == 0: continue richness, inv_simpson = self.diversity(crops_features, "CODE_CULTU", cell) feature = QgsFeature(metrics.fields()) feature.setGeometry(cell) feature["RICHNESS"] = richness feature["DIVERSITY"] = inv_simpson feature["FIELD_SIZE"] = self.mean_field_size(crops_features, cell) feature["SEMI_NATURAL"] = self.semi_natural_cover( cell, crops_features, artificial_features ) new_features.append(feature) multi_feedback.setProgress(int((i / total_feats) * 100)) metrics.dataProvider().addFeatures(new_features) (sink, dest_id) = self.parameterAsSink( parameters, self.METRICS, context, metrics.fields(), metrics.wkbType(), metrics.sourceCrs(), ) for f in metrics.getFeatures(): sink.addFeature(f, QgsFeatureSink.FastInsert) return {self.METRICS: dest_id} def name(self): return "landscapebiodiv" def displayName(self): return "Landscape Biodiversity" def group(self): return "Custom processing" def groupId(self): return "custom_processing" def createInstance(self): return LandscapeBiodiversityMetrics() def diversity( self, 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. 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 mean_field_size(self, 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 semi_natural_cover( self, 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 spatial_query( self, 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))