from statistics import mean
from qgis.core import (
QgsProcessingAlgorithm,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFeatureSink,
QgsVectorLayer,
QgsFeatureSink,
QgsField,
QgsRectangle,
QgsFeatureRequest,
QgsFeature,
QgsGeometry,
QgsProcessingParameterField,
QgsFields,
QgsProcessingParameterNumber,
Qgis,
)
import processing
from PyQt5.QtCore import QVariant
def spatial_query_clip(layer: QgsVectorLayer, extent: QgsRectangle) -> list[QgsFeature]:
"""Query features that intersect given extent and clip geometries.
:param layer: Layer of features to query.
:type layer: QgsVectorLayer
:param extent: Extent used in spatial query.
:type extent: QgsRectangle
:return: Features with clipped geometries.
:rtype: list[QgsFeature]
"""
request = QgsFeatureRequest().setFilterRect(extent)
features: list[QgsFeature] = list(layer.getFeatures(request))
clipped_features: list[QgsFeature] = []
for feat in features:
geom = feat.geometry()
geom = geom.intersection(QgsGeometry.fromRect(extent))
feat.setGeometry(geom)
clipped_features.append(feat)
return clipped_features
def richness(features: list[QgsFeature], field: str) -> int:
"""Compute richness of specific field for given features.
:param features: Features to compute richness on.
:type features: list[QgsFeature]
:param field: Field to compute richness.
:type field: str
:return: Richness.
:rtype: int
"""
values = [f[field] for f in features]
return len(set(values))
def simpson_diversity(features: list[QgsFeature], field: str) -> float:
"""Compute Simpson diversity index.
:param features: Features to compute diversity on.
:type features: list[QgsFeature]
:param field: Field to compute diversity.
:type field: str
:return: Simpson diversity index.
:rtype: float
"""
total_area = sum([f.geometry().area() for f in features])
uniques = {f[field] for f in features}
components: list[float] = []
for value in uniques:
area = sum([f.geometry().area() for f in features if f[field] == value])
components.append((area / total_area) ** 2)
return 1 / sum(components)
def mean_field_size(features: list[QgsFeature]) -> float:
"""compute
:param features: Fields features.
:type features: list[QgsFeature]
:return: Mean area of fields.
:rtype: float
"""
return mean([f.geometry().area() for f in features]) / 10000
def semi_natural_cover(
crop_features: list[QgsFeature],
artificial_features: list[QgsFeature],
cell_area: float,
) -> float:
"""Compute semi natural cover on given cell area.
:param crop_features: Features of crops.
:type crop_features: list[QgsFeature]
:param artificial_features: Features of artificial layer.
:type artificial_features: list[QgsFeature]
:param cell_area: Area of cell.
:type cell_area: float
:return: Semi natural cover.
:rtype: float
"""
crop_area = sum([f.geometry().area() for f in crop_features])
artificial_area = sum([f.geometry().area() for f in artificial_features])
snc = 1 - (crop_area + artificial_area) / cell_area
return snc
class LandscapeBiodiversityMetrics(QgsProcessingAlgorithm):
RPG = "RPG"
ARTIFICIAL = "ARTIFICIAL"
METRICS = "METRICS"
CROP_FIELD = "CROP_FIELD"
CELL_SIZE = "CELL_SIZE"
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterFeatureSource(self.ARTIFICIAL, "Artificial layer")
)
self.addParameter(QgsProcessingParameterFeatureSource(self.RPG, "Crop layer"))
self.addParameter(
QgsProcessingParameterField(
self.CROP_FIELD, "Crop field", parentLayerParameterName=self.RPG
)
)
self.addParameter(
QgsProcessingParameterNumber(
self.CELL_SIZE,
"Cell size",
Qgis.ProcessingNumberParameterType.Integer,
defaultValue=1000,
minValue=500,
maxValue=2000,
)
)
self.addParameter(
QgsProcessingParameterFeatureSink(self.METRICS, "Metrics layer")
)
def processAlgorithm(self, parameters, context, feedback):
rpg: QgsVectorLayer = self.parameterAsVectorLayer(parameters, self.RPG, context)
artificial = self.parameterAsVectorLayer(parameters, self.ARTIFICIAL, context)
field = self.parameterAsString(parameters, self.CROP_FIELD, context)
cell_size = self.parameterAsInt(parameters, self.CELL_SIZE, context)
metrics_fields = QgsFields(
[
QgsField("RICHNESS", QVariant.Double),
QgsField("DIVERSITY", QVariant.Double),
QgsField("FIELD_SIZE", QVariant.Double),
QgsField("SEMI_NATURAL", QVariant.Double),
]
)
(sink, sink_id) = self.parameterAsSink(
parameters,
self.METRICS,
context,
metrics_fields,
Qgis.WkbType.Polygon,
rpg.sourceCrs(),
)
grid: QgsVectorLayer = processing.run(
"native:creategrid",
{
"TYPE": 2,
"EXTENT": rpg,
"HSPACING": cell_size,
"VSPACING": cell_size,
"HOVERLAY": 0,
"VOVERLAY": 0,
"CRS": rpg.crs(),
"OUTPUT": "TEMPORARY_OUTPUT",
},
feedback=feedback,
context=context,
)["OUTPUT"]
for i, cell_feat in enumerate(grid.getFeatures()):
cell_extent = cell_feat.geometry().boundingBox()
crop_features = spatial_query_clip(rpg, cell_extent)
if len(crop_features) == 0:
continue
artificial_features = spatial_query_clip(artificial, cell_extent)
cell_area = cell_feat.geometry().area()
mfs = mean_field_size(crop_features)
rich = richness(crop_features, field)
simpson = simpson_diversity(crop_features, field)
snc = semi_natural_cover(crop_features, artificial_features, cell_area)
feature = QgsFeature(metrics_fields)
feature.setGeometry(QgsGeometry.fromRect(cell_extent))
feature["RICHNESS"] = rich
feature["DIVERSITY"] = simpson
feature["FIELD_SIZE"] = mfs
feature["SEMI_NATURAL"] = snc
sink.addFeature(feature)
feedback.setProgress(i / grid.featureCount() * 100)
return {self.METRICS: sink_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()