-
-
Save j08lue/436162ac8eb90241a98b2702416356ac to your computer and use it in GitHub Desktop.
Test rio-tiler coverage-weighted stats
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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.preview`. | |
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 fc |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment