Skip to content


Form Readers

rio-tiler's Readers provide simple .statistics method to retrieve dataset global statistics

with Reader("my.tif") as src:
    stats = src.statistics()

# Statistics result is in form of Dict[str, rio_tiler.models.BandStatistics]
>>> ["b1", "b2", "b3"]

# rio_tiler.models.BandStatistics is a pydantic model
    # Percentile entries depend on user inputs


You can get statistics from ImageData objects which are returned by all rio-tiler reader methods (e.g. .tile(), .preview(), .part(), ...)

with Reader("cog.tif") as src:
    image = src.preview()
    stats = image.statistics()

Area Weighted Statistics

When getting statistics from a feature, you may want to calculate values from the pixels which intersect with the geometry but also take the pixel intersection percentage into account. Starting with rio-tiler 6.2.0, we've added a coverage option to the statistics utility which enable the user to pass an array representing the coverage percentage such as:

import numpy as np
from rio_tiler.utils import get_array_statistics

# Data Array
# 1, 2
# 3, 4
data =, 2, 3, 4)).reshape((1, 2, 2))

# Coverage Array
# 0.5, 0
# 1, 0.25
coverage = np.array((0.5, 0, 1, 0.25)).reshape((2, 2))

stats = get_array_statistics(data, coverage=coverage)
assert len(stats) == 1
assert stats[0]["min"] == 1
assert stats[0]["max"] == 4
assert (
    round(stats[0]["mean"], 4) == 2.5714
)  # sum of weighted array / sum of weights | 4.5 / 1.75 = 2.57
assert stats[0]["count"] == 1.75

Adjusting geometry align_bounds_with_dataset=True

In rio-tiler 6.3,0 a new option has been introduced to reduce artifacts and produce more precise zonal statistics. This option is available in the low-level reader.part() method used in rio-tiler reader's .feature() and .part() methods.

with Reader("cog.tif") as src:
    data = src.feature(

    coverage_array = data.get_coverage_array(

    stats = data.statistics(coverage=coverage_array)

When passing align_bounds_with_dataset=True to the reader.part() method (forwarded from .feature or .part reader methods), rio-tiler will adjust the input geometry bounds to match the input dataset resolution/transform and avoid unnecessary resampling.

Zonal Statistics method

You can easily extend the rio-tiler's reader to add a .zonal_statistics() method as:

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 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(
        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,
        **kwargs: Any,
    ) -> Union[FeatureCollection, Feature]:
        """Return statistics from GeoJSON features.

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

            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(
            coverage_array = data.get_coverage_array(

            stats = data.statistics(

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

        return fc.features[0] if isinstance(geojson, Feature) else fc