Validating rio-tiler's zonal statistics against a geographical coordinates method

rio_tiler_zonalstats_reader.py

"""
Adapted from the rio-tiler docs: https://cogeotiff.github.io/rio-tiler/advanced/statistics/#area-weighted-statistics
"""

import attr
from typing import Any, Union, Optional, List, Dict

from rio_tiler import io
from rio_tiler.models import BandStatistics
from rio_tiler.constants import WGS84_CRS

from rasterio.crs import CRS

from geojson_pydantic.features import Feature, FeatureCollection
from geojson_pydantic.geometries import Polygon


class Reader(io.Reader):
    """Custom Reader with zonal_statistics method."""

    def zonal_statistics(
        self,
        geojson: Union[FeatureCollection, Feature],
        categorical: bool = False,
        categories: Optional[List[float]] = None,
        percentiles: Optional[List[int]] = None,
        hist_options: Optional[Dict] = None,
        max_size: int = None,
        align_bounds_with_dataset: bool = True,
        shape_crs: Optional[CRS] = WGS84_CRS,
        cover_scale: int = 10,
        **kwargs: Any,
    ) -> Union[FeatureCollection, Feature]:
        """Return statistics from GeoJSON features.

        Args:
            geojson (Feature or FeatureCollection): a GeoJSON Feature or FeatureCollection.
            categorical (bool): treat input data as categorical data. Defaults to False.
            categories (list of numbers, optional): list of categories to return value for.
            percentiles (list of numbers, optional): list of percentile values to calculate. Defaults to `[2, 98]`.
            hist_options (dict, optional): Options to forward to numpy.histogram function.
            max_size (int, optional): Limit the size of the longest dimension of the dataset read, respecting bounds X/Y aspect ratio. Defaults to None.
            kwargs (optional): Options to forward to `self.feature`.

        Returns:
            Feature or FeatureCollection

        """
        kwargs = {**self.options, **kwargs}

        hist_options = hist_options or {}

        fc = geojson
        # We transform the input Feature to a FeatureCollection
        if isinstance(fc, Feature):
            fc = FeatureCollection(type="FeatureCollection", features=[geojson])

        for feature in fc:
            geom = feature.model_dump(exclude_none=True)

            # Get data overlapping with the feature (using Reader.feature method)
            data = self.feature(
                geom,
                shape_crs=shape_crs,
                align_bounds_with_dataset=align_bounds_with_dataset,
                max_size=max_size,
                **kwargs,
            )
            coverage_array = data.get_coverage_array(
                geom,
                shape_crs=shape_crs,
                cover_scale=cover_scale,
            )

            stats = data.statistics(
                categorical=categorical,
                categories=categories,
                percentiles=percentiles,
                hist_options=hist_options,
                coverage=coverage_array,
            )
            
            yield coverage_array, stats

            # Update input feature properties and add the statistics
            feature.properties = feature.properties or {}
            feature.properties.update({"statistics": stats})

        #return fc.features[0] if isinstance(geojson, Feature) else fcimport attr
from typing import Any, Union, Optional, List, Dict

from rio_tiler import io
from rio_tiler.models import BandStatistics
from rio_tiler.constants import WGS84_CRS

from geojson_pydantic.features import Feature, FeatureCollection
from geojson_pydantic.geometries import Polygon

zonal-statistics-validation-south-america-total-fluxes.ipynb