TD1 6.1

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

#### Launch QGSI APP


#### Fonctions


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 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.
        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(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(
    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


#### Paramètres

rpg_file = "dev/data/rpg.gpkg"
artificial_file = "dev/data/aritificial.gpkg"

#### Algorithme

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

# Création de la grille
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
new_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")

result_layer.dataProvider().addAttributes(new_fields)
result_layer.updateFields()

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

new_features = []

for cell_feat in grid.getFeatures():
    cell = cell_feat.geometry()
    crops_features = spatial_query(rpg, cell.boundingBox())
    artificial_features = spatial_query(artificial, cell.boundingBox())

    if len(crops_features) == 0:
        continue

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

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

    feature["RICHNESS"] = richness
    feature["DIVERSITY"] = inv_simpson
    feature["FIELD_SIZE"] = mean_field_size(crops_features, cell)
    feature["SEMI_NATURAL"] = semi_natural_cover(
        cell, crops_features, artificial_features
    )
    new_features.append(feature)

result_layer.dataProvider().addFeatures(new_features)