How to Use the Statistics API¶

The /statistics endpoint can be used to get summary statistics for a geojson Feature or FeatureCollection.

This notebook demonstrates how to use the statistics API using the TROPESS Chemical Reanalysis O3 Monthly 3-dimensional Product dataset.

This notebook has two halves. The first covers continuous data through the xarray backend. The rest covers categorical data through the rasterio backend, using the OPERA Land Surface Disturbance Annual product.

Setup¶

In [1]:

import json
import os
from datetime import datetime, timezone

import httpx2 as httpx

titiler_endpoint = os.getenv(
    "TITILER_CMR_ENDPOINT", "https://openveda.cloud/api/titiler-cmr"
)

import json import os from datetime import datetime, timezone import httpx2 as httpx titiler_endpoint = os.getenv( "TITILER_CMR_ENDPOINT", "https://openveda.cloud/api/titiler-cmr" )

The xarray backend¶

Specifiy a geojson dictionary to use in the request.

In [2]:

geojson_dict = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {},
            "geometry": {
                "coordinates": [
                    [
                        [-20.79973248834736, 83.55979308678764],
                        [-20.79973248834736, 75.0115425216471],
                        [14.483337068956956, 75.0115425216471],
                        [14.483337068956956, 83.55979308678764],
                        [-20.79973248834736, 83.55979308678764],
                    ]
                ],
                "type": "Polygon",
            },
        }
    ],
}

geojson_dict = { "type": "FeatureCollection", "features": [ { "type": "Feature", "properties": {}, "geometry": { "coordinates": [ [ [-20.79973248834736, 83.55979308678764], [-20.79973248834736, 75.0115425216471], [14.483337068956956, 75.0115425216471], [14.483337068956956, 83.55979308678764], [-20.79973248834736, 83.55979308678764], ] ], "type": "Polygon", }, } ], }

TiTiler-CMR will use one or more sel parameter values to select specific subsets from a multidimensional data granule. The TROPESS dataset consists of monthly granules with daily arrays. To load the values from 2021-10-10 we provide that timestamp to the temporal parameter then provide sel=time=nearest::{datetime} to have TiTiler-CMR interpolate the value you provided to the temporal parameter. For fuzzy matching you can provide the nearest:: prefix to the value.

In [3]:

r = httpx.post(
    f"{titiler_endpoint}/xarray/statistics",
    params=(
        ("collection_concept_id", "C2837626477-GES_DISC"),
        # Timestamp or temporal range for CMR granule query
        ("temporal", datetime(2021, 10, 10, tzinfo=timezone.utc).isoformat()),
        # xarray backend query parameters
        ("variables", "o3"),
        ("sel", "time=nearest::{datetime}"),  # gets interpolated
        ("sel", "lev=1000"),
    ),
    json=geojson_dict,
    timeout=None,
).json()

print(json.dumps(r, indent=2))

r = httpx.post( f"{titiler_endpoint}/xarray/statistics", params=( ("collection_concept_id", "C2837626477-GES_DISC"), # Timestamp or temporal range for CMR granule query ("temporal", datetime(2021, 10, 10, tzinfo=timezone.utc).isoformat()), # xarray backend query parameters ("variables", "o3"), ("sel", "time=nearest::{datetime}"), # gets interpolated ("sel", "lev=1000"), ), json=geojson_dict, timeout=None, ).json() print(json.dumps(r, indent=2))

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -20.79973248834736,
              83.55979308678764
            ],
            [
              -20.79973248834736,
              75.0115425216471
            ],
            [
              14.483337068956956,
              75.0115425216471
            ],
            [
              14.483337068956956,
              83.55979308678764
            ],
            [
              -20.79973248834736,
              83.55979308678764
            ]
          ]
        ]
      },
      "properties": {
        "statistics": {
          "b1": {
            "min": 12.448402404785156,
            "max": 30.638206481933594,
            "mean": 25.00276756286621,
            "count": 236.0,
            "sum": 5900.6533203125,
            "std": 3.4685499803764728,
            "median": 25.134218215942383,
            "majority": 12.448402404785156,
            "minority": 12.448402404785156,
            "unique": 236.0,
            "histogram": [
              [
                2,
                1,
                0,
                6,
                33,
                39,
                37,
                42,
                32,
                44
              ],
              [
                12.448402404785156,
                14.267382621765137,
                16.086362838745117,
                17.905344009399414,
                19.724323272705078,
                21.543304443359375,
                23.362285614013672,
                25.18126678466797,
                27.000246047973633,
                28.81922721862793,
                30.638206481933594
              ]
            ],
            "valid_percent": 95.16,
            "masked_pixels": 12.0,
            "valid_pixels": 236.0,
            "description": "0",
            "percentile_2": 18.29606056213379,
            "percentile_98": 30.242084503173828
          }
        },
        "used_assets": [
          "TRPSCRO3M3D.1:TROPESS_reanalysis_mon_o3_2021.nc"
        ]
      }
    }
  ]
}

You can chose a different time slice from the same granule simply by updating the datetime query parameter.

In [4]:

r = httpx.post(
    f"{titiler_endpoint}/xarray/statistics",
    params=(
        ("collection_concept_id", "C2837626477-GES_DISC"),
        # Datetime for CMR granule query
        ("temporal", datetime(2021, 12, 10, tzinfo=timezone.utc).isoformat()),
        # xarray backend query parameters
        ("variables", "o3"),
        ("sel", "time=nearest::{datetime}"),
        ("sel", "lev=1000"),
    ),
    json=geojson_dict,
    timeout=None,
).json()

print(json.dumps(r, indent=2))

r = httpx.post( f"{titiler_endpoint}/xarray/statistics", params=( ("collection_concept_id", "C2837626477-GES_DISC"), # Datetime for CMR granule query ("temporal", datetime(2021, 12, 10, tzinfo=timezone.utc).isoformat()), # xarray backend query parameters ("variables", "o3"), ("sel", "time=nearest::{datetime}"), ("sel", "lev=1000"), ), json=geojson_dict, timeout=None, ).json() print(json.dumps(r, indent=2))

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -20.79973248834736,
              83.55979308678764
            ],
            [
              -20.79973248834736,
              75.0115425216471
            ],
            [
              14.483337068956956,
              75.0115425216471
            ],
            [
              14.483337068956956,
              83.55979308678764
            ],
            [
              -20.79973248834736,
              83.55979308678764
            ]
          ]
        ]
      },
      "properties": {
        "statistics": {
          "b1": {
            "min": 18.230709075927734,
            "max": 37.22050476074219,
            "mean": 27.433626174926758,
            "count": 192.0,
            "sum": 5267.25634765625,
            "std": 4.652167252044063,
            "median": 27.250431060791016,
            "majority": 18.230709075927734,
            "minority": 18.230709075927734,
            "unique": 192.0,
            "histogram": [
              [
                1,
                23,
                39,
                21,
                16,
                16,
                30,
                26,
                12,
                8
              ],
              [
                18.230709075927734,
                20.129688262939453,
                22.028667449951172,
                23.92764663696289,
                25.826627731323242,
                27.72560691833496,
                29.62458610534668,
                31.52356719970703,
                33.42254638671875,
                35.32152557373047,
                37.22050476074219
              ]
            ],
            "valid_percent": 77.42,
            "masked_pixels": 56.0,
            "valid_pixels": 192.0,
            "description": "0",
            "percentile_2": 20.44560432434082,
            "percentile_98": 36.73855972290039
          }
        },
        "used_assets": [
          "TRPSCRO3M3D.1:TROPESS_reanalysis_mon_o3_2021.nc"
        ]
      }
    }
  ]
}

The rasterio backend and asset selection¶

Everything from here on uses one dataset:

OPERA Land Surface Disturbance Annual (C2519119034-LPCLOUD, "DIST-ANN") maps vegetation loss across global land. For each 30 m square of ground, the algorithm estimates what fraction is covered by green vegetation, compares that against the same square's own history, and records where cover has dropped. A companion product, DIST-ALERT, publishes this every few days as new imagery arrives; DIST-ANN summarizes a whole calendar year. One granule covers one 110 km tile for one year, and the collection spans 2023 through 2025. The LP DAAC product page has the full specification.

It suits this section for two reasons:

1. Its pixel values are class codes rather than measurements, which is what categorical=true exists for.
2. Each granule bundles many files, which is what the assets parameters exist for.

The xarray backend reads variables out of a single multidimensional file. The rasterio backend reads GDAL-readable rasters, and a CMR granule usually bundles many of them.

One DIST-ANN granule is 21 separate Cloud Optimized GeoTIFFs, one per measurement:

..._VEG-DIST-STATUS.tif    <- class codes; the one we want
..._VEG-ANOM-MAX.tif       <- how much vegetation cover was lost
..._VEG-HIST.tif           <- the baseline cover it was measured against
..._VEG-DIST-DATE.tif      <- when the loss was first detected

titiler-cmr calls each of these an asset. Two parameters work together to pick one, and they do genuinely different jobs:

• assets_regex decides what the assets are called. Each file's name is matched against the pattern; files that do not match are dropped, and for the ones that do, the matched substring becomes the asset key.
• assets then selects by key.

The first is doing more than it looks. Without a regex, titiler-cmr has no names to work with, so it falls back to numbering the files by their position in CMR's listing: "0", "1", ... "20". So assets_regex is what gives the assets stable, meaningful names in the first place.

file in the granule             assets_regex       asset key (matched part)
-----------------------------   -------------      ------------------------
..._30_v1_VEG-DIST-STATUS.tif   VEG-[A-Z-]+   ->   VEG-DIST-STATUS
..._30_v1_VEG-ANOM-MAX.tif      VEG-[A-Z-]+   ->   VEG-ANOM-MAX
..._30_v1_VEG-HIST.tif          VEG-[A-Z-]+   ->   VEG-HIST
                                                                 ...

then: assets=["VEG-DIST-STATUS"] reads the one asset (file) with this key

Seen that way, assets_regex=VEG-DIST-STATUS alongside assets=VEG-DIST-STATUS stops looking redundant. The first says keep only the status file, and call it VEG-DIST-STATUS; the second says read the asset called VEG-DIST-STATUS. They carry the same string only because we want a single file and named it after itself. Widen the regex and it stops resembling the assets list at all.

In [5]:

# OPERA Land Surface Disturbance Annual, from Harmonized Landsat Sentinel-2.
# Rasterio backend, 30 m COGs, annual granules for 2023 through 2025.
opera_collection = "C2519119034-LPCLOUD"
opera_year = "2024-01-01T00:00:00Z/2024-12-31T23:59:59Z"

# A small area of interest over the eastern San Francisco Bay.
aoi_bay_area = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [-122.30, 37.80],
                [-122.30, 37.95],
                [-122.10, 37.95],
                [-122.10, 37.80],
                [-122.30, 37.80],
            ]
        ],
    },
}


def opera_statistics(feature, temporal, *, endpoint="/rasterio/statistics", **params):
    """POST a GeoJSON feature to a statistics endpoint and return the parsed JSON."""
    r = httpx.post(
        f"{titiler_endpoint}{endpoint}",
        params={
            "collection_concept_id": opera_collection,
            "temporal": temporal,
            **params,
        },
        json=feature,
        timeout=None,
    )
    r.raise_for_status()

    return r.json()

# OPERA Land Surface Disturbance Annual, from Harmonized Landsat Sentinel-2. # Rasterio backend, 30 m COGs, annual granules for 2023 through 2025. opera_collection = "C2519119034-LPCLOUD" opera_year = "2024-01-01T00:00:00Z/2024-12-31T23:59:59Z" # A small area of interest over the eastern San Francisco Bay. aoi_bay_area = { "type": "Feature", "properties": {}, "geometry": { "type": "Polygon", "coordinates": [ [ [-122.30, 37.80], [-122.30, 37.95], [-122.10, 37.95], [-122.10, 37.80], [-122.30, 37.80], ] ], }, } def opera_statistics(feature, temporal, *, endpoint="/rasterio/statistics", **params): """POST a GeoJSON feature to a statistics endpoint and return the parsed JSON.""" r = httpx.post( f"{titiler_endpoint}{endpoint}", params={ "collection_concept_id": opera_collection, "temporal": temporal, **params, }, json=feature, timeout=None, ) r.raise_for_status() return r.json()

In [6]:

r = opera_statistics(
    aoi_bay_area, opera_year, assets_regex="VEG-DIST-STATUS", assets="VEG-DIST-STATUS"
)

# `statistics` is keyed by band; this collection's asset is a single band, `b1`.
print(json.dumps(r["properties"]["statistics"], indent=2))

r = opera_statistics( aoi_bay_area, opera_year, assets_regex="VEG-DIST-STATUS", assets="VEG-DIST-STATUS" ) # `statistics` is keyed by band; this collection's asset is a single band, `b1`. print(json.dumps(r["properties"]["statistics"], indent=2))

{
  "b1": {
    "min": 0.0,
    "max": 10.0,
    "mean": 0.05040227249264717,
    "count": 319888.375,
    "sum": 16123.1005859375,
    "std": 0.5926251466974075,
    "median": 0.0,
    "majority": 0.0,
    "minority": 10.0,
    "unique": 7.0,
    "histogram": [
      [
        318114,
        0,
        0,
        120,
        0,
        0,
        35,
        2049,
        67,
        75
      ],
      [
        0.0,
        1.0,
        2.0,
        3.0,
        4.0,
        5.0,
        6.0,
        7.0,
        8.0,
        9.0,
        10.0
      ]
    ],
    "valid_percent": 100.0,
    "masked_pixels": 0.0,
    "valid_pixels": 320460.0,
    "description": "VEG-DIST-STATUS_Vegetation_disturbance_status",
    "percentile_2": 0.0,
    "percentile_98": 0.0
  }
}

Here the two part company: a wider regex puts several assets on the menu, and assets picks two of them.

In [7]:

several = opera_statistics(
    aoi_bay_area,
    opera_year,
    assets_regex="VEG-[A-Z-]+",
    assets=["VEG-DIST-STATUS", "VEG-ANOM-MAX"],
)

for band, stats in several["properties"]["statistics"].items():
    print(f"{band}  {stats['description']}")

several = opera_statistics( aoi_bay_area, opera_year, assets_regex="VEG-[A-Z-]+", assets=["VEG-DIST-STATUS", "VEG-ANOM-MAX"], ) for band, stats in several["properties"]["statistics"].items(): print(f"{band} {stats['description']}")

b1  VEG-DIST-STATUS_Vegetation_disturbance_status
b2  VEG-ANOM-MAX_Maximum_vegetation_loss_anomaly

Two assets in, two bands out, in the order assets listed them. The regex defined the menu; assets chose from it.

Each description is two parts: the asset key, an underscore, then the band description stored inside the GeoTIFF itself. So VEG-ANOM-MAX is the key your regex produced, and Maximum_vegetation_loss_anomaly is text the data producer wrote into the file. The prefix is how you tell which file a band came from, which makes it worth reading before you trust a result.

Leaving the asset parameters off altogether does not fail, which is the part worth seeing:

In [8]:

# No asset parameters at all. The request still succeeds, on the wrong terms.
fallback = opera_statistics(aoi_bay_area, opera_year)
print(fallback["properties"]["statistics"]["b1"]["description"])

# No asset parameters at all. The request still succeeds, on the wrong terms. fallback = opera_statistics(aoi_bay_area, opera_year) print(fallback["properties"]["statistics"]["b1"]["description"])

0_Vegetation_disturbance_status

0_Vegetation_disturbance_status, not VEG-DIST-STATUS_Vegetation_disturbance_status. With no regex the assets are numbered rather than named, so the reader fell back to asset "0", whichever file CMR happened to list first. Here that is coincidentally the file we want, which is exactly what makes it a silent trap on any other multi-asset collection. Always name the asset explicitly.

Two other ways to get this wrong, both of which do fail loudly:

• A full-filename regex matches nothing. assets_regex is applied to the granule filename after the extension has been stripped, so a pattern anchored on \.tif$ can never match. Zero assets survive and the request 400s with the rather misleading No valid asset found. Asset's media types not supported. Match the bare layer name instead. This matters in practice: catalogs sometimes publish an assets_regex in exactly that form.

• assets_regex without assets fails as well. The regex renamed the available assets, so the "0" fallback no longer resolves and the API answers 404 with {"detail": "'0' is not valid, should be one of ['VEG-DIST-STATUS']"}.

  The helper above calls raise_for_status(), so what you actually see is an HTTPStatusError. One caveat worth recognizing: a gateway in front of the API may swallow the 404 and serve its own HTML page with a 200 status, in which case the call fails a moment later while decoding JSON instead. That is a bug in the deployment rather than a behavior of the API, and this notebook deliberately does not work around it.

Categorical statistics¶

The pixel values in VEG-DIST-STATUS are class codes, not measurements. 6 is not twice 3; they are labels wearing numbers.

By default, /statistics builds a histogram the continuous way: take the data range, cut it into ten equal-width buckets, and count what falls in each. For class codes that is meaningless.

categorical=true replaces that with one bin per distinct value, counted exactly.

In [9]:

# Settled from here on: name the status layer, and select it.
dist_status_asset = {"assets_regex": "VEG-DIST-STATUS", "assets": "VEG-DIST-STATUS"}

default_stats = opera_statistics(aoi_bay_area, opera_year, **dist_status_asset)
categorical_stats = opera_statistics(
    aoi_bay_area, opera_year, **dist_status_asset, categorical=True
)

for label, response in [
    ("categorical unset", default_stats),
    ("categorical=true", categorical_stats),
]:
    band = response["properties"]["statistics"]["b1"]
    counts, bins = band["histogram"]
    print(label)
    print(f"  counts: {counts}")
    print(f"  bins  : {bins}\n")

# Settled from here on: name the status layer, and select it. dist_status_asset = {"assets_regex": "VEG-DIST-STATUS", "assets": "VEG-DIST-STATUS"} default_stats = opera_statistics(aoi_bay_area, opera_year, **dist_status_asset) categorical_stats = opera_statistics( aoi_bay_area, opera_year, **dist_status_asset, categorical=True ) for label, response in [ ("categorical unset", default_stats), ("categorical=true", categorical_stats), ]: band = response["properties"]["statistics"]["b1"] counts, bins = band["histogram"] print(label) print(f" counts: {counts}") print(f" bins : {bins}\n")

categorical unset
  counts: [318114, 0, 0, 120, 0, 0, 35, 2049, 67, 75]
  bins  : [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]

categorical=true
  counts: [318114, 120, 35, 2049, 67, 70, 5]
  bins  : [0, 3, 6, 7, 8, 9, 10]

categorical=true changes the histogram and nothing else. Every other field is identical between the two responses above.

That matters more than it looks. The response still carries mean, median, std, sum, and percentiles, and those are arithmetic performed on class codes. They always come back, they always look plausible, and they never mean anything for this kind of layer.

In [10]:

band = categorical_stats["properties"]["statistics"]["b1"]

meaningful = [
    "majority",
    "minority",
    "unique",
    "valid_pixels",
    "masked_pixels",
    "valid_percent",
]
meaningless = [
    "min",
    "max",
    "mean",
    "median",
    "std",
    "sum",
    "percentile_2",
    "percentile_98",
]

print("Meaningful for a class-code layer:")

for key in meaningful:
    print(f"  {key:<15}{band[key]}")

print("\nReturned, but MEANINGLESS (arithmetic over class codes):")

for key in meaningless:
    print(f"  {key:<15}{band[key]}")

band = categorical_stats["properties"]["statistics"]["b1"] meaningful = [ "majority", "minority", "unique", "valid_pixels", "masked_pixels", "valid_percent", ] meaningless = [ "min", "max", "mean", "median", "std", "sum", "percentile_2", "percentile_98", ] print("Meaningful for a class-code layer:") for key in meaningful: print(f" {key:<15}{band[key]}") print("\nReturned, but MEANINGLESS (arithmetic over class codes):") for key in meaningless: print(f" {key:<15}{band[key]}")

Meaningful for a class-code layer:
  majority       0.0
  minority       10.0
  unique         7.0
  valid_pixels   320460.0
  masked_pixels  0.0
  valid_percent  100.0

Returned, but MEANINGLESS (arithmetic over class codes):
  min            0.0
  max            10.0
  mean           0.05040227249264717
  median         0.0
  std            0.5926251466974075
  sum            16123.1005859375
  percentile_2   0.0
  percentile_98  0.0

c is read only when categorical=true. Sent on its own it is silently ignored and you get the ten equal-width bins back, with no error. The c parameter is covered below.

To turn counts into percentages, divide by valid_pixels, and not by the sum of the counts. There are two reasons.

First, c can narrow which classes get reported. When it does, the counts sum to less than valid_pixels, and dividing by that sum silently inflates every class.

Second, the response mixes two families of number:

• Unweighted, a plain tally of pixels: valid_pixels, masked_pixels, and the categorical histogram counts.
• Coverage-weighted: count, sum, mean, std, and the percentiles. Your polygon's edge does not follow the pixel grid, so each pixel is weighted by the fraction of its area that falls inside the polygon. A pixel lying half inside contributes 0.5 rather than 1.

That is why count above came back fractional, and a little below valid_pixels. The difference is entirely an edge effect, so it grows as the area of interest shrinks relative to a 30 m pixel: negligible over a city, but several percent over a few hectares. Dividing an unweighted count by a weighted one mixes the two families and bakes that error in.

Reading the OPERA class codes¶

A categorical histogram is only as useful as your understanding of the codes. For VEG-DIST-STATUS, verbatim from the OPERA DIST Product Specification Document v1.0, Table 4-2:

Value	Official label	Plain English
0	No disturbance	Assessed; no confirmed event this year
3	confirmed <50% ongoing	Loss under 50 points, still accumulating at year end
6	confirmed >=50% ongoing	Loss of 50+ points, still accumulating at year end
7	confirmed <50% finished	Loss under 50 points, stopped before year end
8	confirmed >=50% finished	Loss of 50+ points, stopped before year end
9	confirmed previous year <50%	Under 50 points, first detected the previous calendar year
10	confirmed previous year >=50%	50+ points, first detected the previous calendar year
255	No data	Fill. Not a class.

It is a grid, not a ramp¶

The codes encode two independent axes, lifecycle and intensity:

lifecycle	<50%	>=50%
ongoing at year end	3	6
finished before year end	7	8
first detected previous year	9	10

So the integers are not ordered by severity. Class 8 (50+ points lost, event finished) is more severe than class 9 (under 50 points, carried over from last year), despite the lower number. A continuous color ramp over 0-10 draws the more severe pixel lighter. Use a discrete colormap keyed on the seven valid values, with 255 transparent.

What the percentage thresholds actually measure¶

Absolute percentage points of ground cover, not relative vegetation loss. The product estimates what fraction of each pixel's 900 m2 is covered by green vegetation, and compares it to a baseline built from that same pixel's own history. The anomaly is the subtraction:

anomaly = baseline_cover_percent - current_cover_percent

A pixel dropping from 80% cover to 35% cover lost 45 points, so it is class 3 or 7, even though it lost 56% of its vegetation. Legend text reading "more than half the vegetation lost" is wrong.

A corollary that surprises people: a pixel whose baseline cover is under 50% can never reach class 6 or 8, because losing everything is still under 50 points. Sparse grassland scraped to bare dirt classes as low intensity.

Class 0 does not mean "nothing happened"¶

It means "no event reached confirmed status". A real detection that stalled and was revoked reports 0, indistinguishable from untouched ground.

The same integers mean different things in sibling collections¶

Value	DIST-ANN v1	DIST-ALERT v1	DIST-ALERT v0
1	(undefined)	first <50%	provisional <50%
2	(undefined)	provisional <50%	confirmed <50%
3	confirmed <50% ongoing	confirmed <50%	provisional >=50%
4	(undefined)	first >=50%	confirmed >=50%
6	confirmed >=50% ongoing	confirmed >=50%	(undefined)

Value 3 is low intensity in two of those and high intensity in the third. Key every legend and colormap to the specific collection concept ID; never share one across the DIST family.

Note also that NASA GIBS's published colormap for the annual layer defines only 0, 3, 6, 7, 8, and 255. It has no entry for 9 or 10, so adopting it wholesale leaves some pixels unrenderable.

In [11]:

# The complete VEG-DIST-STATUS encoding, keyed to OPERA_L3_DIST-ANN-HLS_V1.
# Do NOT reuse for DIST-ALERT: the same integers mean different things there.
DIST_ANN_NODATA = 255
DIST_ANN_CLASSES = {
    0: "no confirmed disturbance",
    3: "confirmed <50 pts, ongoing at year end",
    6: "confirmed >=50 pts, ongoing at year end",
    7: "confirmed <50 pts, finished before year end",
    8: "confirmed >=50 pts, finished before year end",
    9: "confirmed <50 pts, first detected previous year",
    10: "confirmed >=50 pts, first detected previous year",
    DIST_ANN_NODATA: "no data (fill; masked out before statistics run)",
}

band = categorical_stats["properties"]["statistics"]["b1"]
counts, values = band["histogram"]
valid_pixels = band["valid_pixels"]

print(f"{'value':>5} {'pixels':>9} {'% of valid':>11}   class")

for value, count in zip(values, counts):
    share = count / valid_pixels
    label = DIST_ANN_CLASSES.get(int(value), "unrecognized class")
    print(f"{value:>5} {count:>9,} {share:>11.2%}   {label}")

disturbed = sum(c for v, c in zip(values, counts) if v != 0)

print(f"\nany confirmed disturbance: {disturbed:,} px ({disturbed / valid_pixels:.2%})")

# The complete VEG-DIST-STATUS encoding, keyed to OPERA_L3_DIST-ANN-HLS_V1. # Do NOT reuse for DIST-ALERT: the same integers mean different things there. DIST_ANN_NODATA = 255 DIST_ANN_CLASSES = { 0: "no confirmed disturbance", 3: "confirmed <50 pts, ongoing at year end", 6: "confirmed >=50 pts, ongoing at year end", 7: "confirmed <50 pts, finished before year end", 8: "confirmed >=50 pts, finished before year end", 9: "confirmed <50 pts, first detected previous year", 10: "confirmed >=50 pts, first detected previous year", DIST_ANN_NODATA: "no data (fill; masked out before statistics run)", } band = categorical_stats["properties"]["statistics"]["b1"] counts, values = band["histogram"] valid_pixels = band["valid_pixels"] print(f"{'value':>5} {'pixels':>9} {'% of valid':>11} class") for value, count in zip(values, counts): share = count / valid_pixels label = DIST_ANN_CLASSES.get(int(value), "unrecognized class") print(f"{value:>5} {count:>9,} {share:>11.2%} {label}") disturbed = sum(c for v, c in zip(values, counts) if v != 0) print(f"\nany confirmed disturbance: {disturbed:,} px ({disturbed / valid_pixels:.2%})")

value    pixels  % of valid   class
    0   318,114      99.27%   no confirmed disturbance
    3       120       0.04%   confirmed <50 pts, ongoing at year end
    6        35       0.01%   confirmed >=50 pts, ongoing at year end
    7     2,049       0.64%   confirmed <50 pts, finished before year end
    8        67       0.02%   confirmed >=50 pts, finished before year end
    9        70       0.02%   confirmed <50 pts, first detected previous year
   10         5       0.00%   confirmed >=50 pts, first detected previous year

any confirmed disturbance: 2,346 px (0.73%)

Aligning timesteps with the c parameter¶

Without c, the bins are "whatever classes happen to appear in this area of interest, at this timestep". That is fine for a single request and a problem for a timeseries, because the bin sets do not line up.

/rasterio/timeseries/statistics expands temporal plus step plus temporal_mode into one /statistics request per interval, and returns them keyed by interval.

In [12]:

# The full span of the collection, sampled one calendar year at a time.
opera_timespan = "2023-01-01T00:00:00Z/2025-12-31T23:59:59Z"
annual_steps = {"step": "P1Y", "temporal_mode": "interval"}

ragged = opera_statistics(
    aoi_bay_area,
    opera_timespan,
    endpoint="/rasterio/timeseries/statistics",
    **dist_status_asset,
    categorical=True,
    **annual_steps,
)

print("without `c`, the bin sets differ per year:")

for interval, stats in ragged["properties"]["statistics"].items():
    counts, values = stats["b1"]["histogram"]
    print(f"  {interval[:10]}  bins={values}")

# The full span of the collection, sampled one calendar year at a time. opera_timespan = "2023-01-01T00:00:00Z/2025-12-31T23:59:59Z" annual_steps = {"step": "P1Y", "temporal_mode": "interval"} ragged = opera_statistics( aoi_bay_area, opera_timespan, endpoint="/rasterio/timeseries/statistics", **dist_status_asset, categorical=True, **annual_steps, ) print("without `c`, the bin sets differ per year:") for interval, stats in ragged["properties"]["statistics"].items(): counts, values = stats["b1"]["histogram"] print(f" {interval[:10]} bins={values}")

without `c`, the bin sets differ per year:
  2023-01-01  bins=[0, 3, 6, 7, 8]
  2024-01-01  bins=[0, 3, 6, 7, 8, 9, 10]
  2025-01-01  bins=[0, 3, 6, 7, 8, 9, 10]

In [13]:

# Derived from the class table above, minus the fill value: the COG
# declares 255 as nodata so those pixels are masked before statistics run, and
# asking for it with `c` would only ever yield an empty bin.
DIST_ANN_CLASS_VALUES = [v for v in DIST_ANN_CLASSES if v != DIST_ANN_NODATA]
print("requesting classes:", DIST_ANN_CLASS_VALUES)

aligned = opera_statistics(
    aoi_bay_area,
    opera_timespan,
    endpoint="/rasterio/timeseries/statistics",
    **dist_status_asset,
    categorical=True,
    c=DIST_ANN_CLASS_VALUES,
    **annual_steps,
)

header = "  ".join(f"{v:>7}" for v in DIST_ANN_CLASS_VALUES)
print(f"{'interval':<12}{header}")

for interval, stats in aligned["properties"]["statistics"].items():
    counts, values = stats["b1"]["histogram"]
    row = "  ".join(f"{c:>7,}" for c in counts)
    print(f"{interval[:10]:<12}{row}")

# Derived from the class table above, minus the fill value: the COG # declares 255 as nodata so those pixels are masked before statistics run, and # asking for it with `c` would only ever yield an empty bin. DIST_ANN_CLASS_VALUES = [v for v in DIST_ANN_CLASSES if v != DIST_ANN_NODATA] print("requesting classes:", DIST_ANN_CLASS_VALUES) aligned = opera_statistics( aoi_bay_area, opera_timespan, endpoint="/rasterio/timeseries/statistics", **dist_status_asset, categorical=True, c=DIST_ANN_CLASS_VALUES, **annual_steps, ) header = " ".join(f"{v:>7}" for v in DIST_ANN_CLASS_VALUES) print(f"{'interval':<12}{header}") for interval, stats in aligned["properties"]["statistics"].items(): counts, values = stats["b1"]["histogram"] row = " ".join(f"{c:>7,}" for c in counts) print(f"{interval[:10]:<12}{row}")

requesting classes: [0, 3, 6, 7, 8, 9, 10]

interval          0        3        6        7        8        9       10
2023-01-01  313,251       74       39    3,055      293        0        0
2024-01-01  318,114      120       35    2,049       67       70        5
2025-01-01  312,933       85       34    2,405       34      972      249

Every year now reports the same seven slots in the same order, with zeros where a class is absent. That is what a chart needs.

Note what is not in DIST_ANN_CLASS_VALUES: the nodata value 255. the source COG declares 255 as its nodata value, so rasterio masks those pixels out before statistics are computed. 255 therefore never appears in a categorical histogram, and asking for it with c produces a bin that is permanently zero. (Passing nodata=255 explicitly changes nothing, which is how you can tell the file already declares it.) This is easy to get wrong, because catalog colormaps for this collection are commonly keyed {0, 3, 6, 7, 8, 9, 10, 255}. If you derive c from colormap keys, strip the fill value first.

In [14]:

# An area of interest over the Pacific, where the granule is all fill.
aoi_ocean = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [-122.80, 37.45],
                [-122.80, 37.60],
                [-122.62, 37.60],
                [-122.62, 37.45],
                [-122.80, 37.45],
            ]
        ],
    },
}

for label, feature in [("land", aoi_bay_area), ("ocean", aoi_ocean)]:
    response = opera_statistics(
        feature, opera_year, **dist_status_asset, categorical=True
    )
    band = response["properties"]["statistics"]["b1"]

    # masked_pixels counts BOTH nodata and everything inside the feature's
    # bounding box but outside the polygon itself.
    print(
        f"{label:<7} valid_pixels={band['valid_pixels']:>10,.0f}  "
        f"masked_pixels={band['masked_pixels']:>10,.0f}  "
        f"bins={band['histogram'][1]}"
    )

# An area of interest over the Pacific, where the granule is all fill. aoi_ocean = { "type": "Feature", "properties": {}, "geometry": { "type": "Polygon", "coordinates": [ [ [-122.80, 37.45], [-122.80, 37.60], [-122.62, 37.60], [-122.62, 37.45], [-122.80, 37.45], ] ], }, } for label, feature in [("land", aoi_bay_area), ("ocean", aoi_ocean)]: response = opera_statistics( feature, opera_year, **dist_status_asset, categorical=True ) band = response["properties"]["statistics"]["b1"] # masked_pixels counts BOTH nodata and everything inside the feature's # bounding box but outside the polygon itself. print( f"{label:<7} valid_pixels={band['valid_pixels']:>10,.0f} " f"masked_pixels={band['masked_pixels']:>10,.0f} " f"bins={band['histogram'][1]}" )

land    valid_pixels=   320,460  masked_pixels=         0  bins=[0, 3, 6, 7, 8, 9, 10]

ocean   valid_pixels=         8  masked_pixels=   289,191  bins=[0]

What changes the numbers: max_size, reproject, and dst_crs¶

/statistics does not read a stored grid and summarize it. It computes an output grid from your request (the bounding box of your feature, in the CRS you asked for), then reads each matching granule onto that grid and summarizes the result.

That has a consequence worth stating plainly: pixel counts are a property of the request, not of the data. Three parameters move them.

The max_size parameter¶

max_size caps the longest dimension of the array that gets read. It defaults to None, meaning no cap: the read happens at the source's native resolution, and the array grows with the area of your feature. It is a ceiling, not a target: if the native read is already smaller, max_size does nothing.

In [15]:

for max_size in [None, 1024, 256]:
    size_param = {} if max_size is None else {"max_size": max_size}
    response = opera_statistics(
        aoi_bay_area, opera_year, **dist_status_asset, categorical=True, **size_param
    )

    band = response["properties"]["statistics"]["b1"]
    counts, values = band["histogram"]
    by_class = dict(zip((int(v) for v in values), counts))

    print(
        f"max_size={str(max_size):<5} valid_pixels={band['valid_pixels']:>9,.0f}  "
        f"class 7={by_class.get(7, 0):>6,}  class 8={by_class.get(8, 0):>5,}  "
        f"class 10={by_class.get(10, 0):>4,}"
    )

for max_size in [None, 1024, 256]: size_param = {} if max_size is None else {"max_size": max_size} response = opera_statistics( aoi_bay_area, opera_year, **dist_status_asset, categorical=True, **size_param ) band = response["properties"]["statistics"]["b1"] counts, values = band["histogram"] by_class = dict(zip((int(v) for v in values), counts)) print( f"max_size={str(max_size):<5} valid_pixels={band['valid_pixels']:>9,.0f} " f"class 7={by_class.get(7, 0):>6,} class 8={by_class.get(8, 0):>5,} " f"class 10={by_class.get(10, 0):>4,}" )

max_size=None  valid_pixels=  320,460  class 7= 2,049  class 8=   67  class 10=   5

max_size=1024  valid_pixels=  320,460  class 7= 2,049  class 8=   67  class 10=   5

max_size=256   valid_pixels=   49,152  class 7=   329  class 8=    9  class 10=   1

max_size=1024 returns numbers identical to the uncapped read, because the native read was already smaller than 1024 on its longest side.

max_size=256 is a different story. Downsampling uses nearest-neighbor, which subsamples rather than aggregates: it keeps every Nth pixel and discards the rest. Common classes survive as percentages. Rare classes become noise and can disappear entirely. And for a question like "how much was disturbed", the disturbed classes are the rare ones.

So: use max_size to bound cost when the area of interest is large, and treat the resulting counts as estimates rather than pixel counts. Report percentages rather than counts when you do.

Why max_size matters in a Lambda¶

titiler-cmr is typically deployed as an AWS Lambda function. Both the memory ceiling and the request timeout are deployment settings, so check yours; the reference deployment allows 10 GB and 30 seconds.

Uncapped, the array scales with area / (30 m)^2, and with the number of intersecting granules, since each one is read onto the same grid before they are merged. max_size=N bounds it at N x N regardless of the area of interest, so peak array memory is roughly N**2 * itemsize * n_bands bytes. VEG-DIST-STATUS is uint8, so 1 byte per pixel; a float32 variable costs four times as much, and rio-tiler carries a boolean mask alongside the data.

The cell below reads a roughly 1-degree area of interest both ways.

In [16]:

import time

# A roughly 1-degree area of interest, about 25x the one used so far.
aoi_large = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [-122.6, 37.2],
                [-122.6, 38.2],
                [-121.6, 38.2],
                [-121.6, 37.2],
                [-122.6, 37.2],
            ]
        ],
    },
}

for max_size in [None, 1024]:
    size_param = {} if max_size is None else {"max_size": max_size}
    started = time.perf_counter()
    response = opera_statistics(
        aoi_large, opera_year, **dist_status_asset, categorical=True, **size_param
    )
    elapsed = time.perf_counter() - started
    band = response["properties"]["statistics"]["b1"]
    total_px = band["valid_pixels"] + band["masked_pixels"]

    print(
        f"max_size={str(max_size):<5} {elapsed:>6.1f} s   "
        f"array={total_px:>12,.0f} px  ~{total_px / 1e6:>5.1f} MB at 1 byte/px"
    )

import time # A roughly 1-degree area of interest, about 25x the one used so far. aoi_large = { "type": "Feature", "properties": {}, "geometry": { "type": "Polygon", "coordinates": [ [ [-122.6, 37.2], [-122.6, 38.2], [-121.6, 38.2], [-121.6, 37.2], [-122.6, 37.2], ] ], }, } for max_size in [None, 1024]: size_param = {} if max_size is None else {"max_size": max_size} started = time.perf_counter() response = opera_statistics( aoi_large, opera_year, **dist_status_asset, categorical=True, **size_param ) elapsed = time.perf_counter() - started band = response["properties"]["statistics"]["b1"] total_px = band["valid_pixels"] + band["masked_pixels"] print( f"max_size={str(max_size):<5} {elapsed:>6.1f} s " f"array={total_px:>12,.0f} px ~{total_px / 1e6:>5.1f} MB at 1 byte/px" )

max_size=None    11.2 s   array=  10,578,756 px  ~ 10.6 MB at 1 byte/px

max_size=1024     1.8 s   array=   1,048,576 px  ~  1.0 MB at 1 byte/px

Note which limit actually binds. At 10 GB of Lambda memory, an array of this size is nowhere near the memory ceiling, but the uncapped read spends a large fraction of the 30 second timeout. In practice you will hit the timeout long before you hit the memory limit, and max_size bounds both.

/timeseries/statistics applies a size guard of its own, before it dispatches any sub-requests. It is judged on the area of interest and the number of timesteps, and bounded by the configurable time_series_max_image_size and time_series_statistics_max_total_size settings.

Asking for a smaller output does not get you past it: neither max_size nor an explicit width and height changes the verdict. The rejection says so:

The AOI for this request is too large for the /statistics endpoint for this
dataset. Try again with either a smaller AOI

The remedy is a smaller area of interest, or fewer timesteps.

The resampling and reproject parameters¶

There are two separate resampling knobs, and they are easy to confuse:

• resampling controls the decimated read, i.e., what happens when max_size shrinks the array.
• reproject controls the warp kernel, i.e., what happens when the output CRS differs from the granule's own CRS.

Both default to nearest, which preserves class codes exactly. Leave both alone for any class-code layer.

A granule is warped whenever the output CRS differs from its native UTM zone, which is the default case, since the default output CRS is geographic. Asking for the granule's own UTM zone skips the warp entirely.

Here is what happens if you change reproject.

In [17]:

for method in ["nearest", "bilinear", "average"]:
    response = opera_statistics(
        aoi_bay_area,
        opera_year,
        **dist_status_asset,
        categorical=True,
        reproject=method,
    )

    bins = response["properties"]["statistics"]["b1"]["histogram"][1]
    invalid = [b for b in bins if int(b) not in DIST_ANN_CLASSES]

    print(f"reproject={method:<9} bins={bins}  INVALID: {invalid}" if invalid else "")

for method in ["nearest", "bilinear", "average"]: response = opera_statistics( aoi_bay_area, opera_year, **dist_status_asset, categorical=True, reproject=method, ) bins = response["properties"]["statistics"]["b1"]["histogram"][1] invalid = [b for b in bins if int(b) not in DIST_ANN_CLASSES] print(f"reproject={method:<9} bins={bins} INVALID: {invalid}" if invalid else "")

reproject=bilinear  bins=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]  INVALID: [1, 2, 4, 5]

reproject=average   bins=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]  INVALID: [1, 2, 4, 5]

bilinear and average interpolate between class codes, manufacturing values that do not exist in this product. Worse, 1, 2, 4, and 5 are undefined in DIST-ANN but are meaningful in DIST-ALERT, so the corrupted output looks like valid data from a sibling collection.

The dst_crs parameter¶

dst_crs chooses the projection of the output grid, and different projections mean differently-shaped pixels covering different amounts of ground.

In [18]:

for label, dst_crs in [
    ("default (EPSG:4326)", None),
    ("equal-area (+proj=cea)", "+proj=cea"),
    ("native UTM (EPSG:32610)", "EPSG:32610"),
]:
    crs_param = {} if dst_crs is None else {"dst_crs": dst_crs}
    response = opera_statistics(
        aoi_bay_area, opera_year, **dist_status_asset, categorical=True, **crs_param
    )

    band = response["properties"]["statistics"]["b1"]
    counts, values = band["histogram"]
    by_class = dict(zip((int(v) for v in values), counts))
    valid = band["valid_pixels"]

    print(
        f"{label:<24} valid_pixels={valid:>9,.0f}  "
        f"class 7={by_class.get(7, 0):>6,} ({by_class.get(7, 0) / valid:.3%})"
    )

for label, dst_crs in [ ("default (EPSG:4326)", None), ("equal-area (+proj=cea)", "+proj=cea"), ("native UTM (EPSG:32610)", "EPSG:32610"), ]: crs_param = {} if dst_crs is None else {"dst_crs": dst_crs} response = opera_statistics( aoi_bay_area, opera_year, **dist_status_asset, categorical=True, **crs_param ) band = response["properties"]["statistics"]["b1"] counts, values = band["histogram"] by_class = dict(zip((int(v) for v in values), counts)) valid = band["valid_pixels"] print( f"{label:<24} valid_pixels={valid:>9,.0f} " f"class 7={by_class.get(7, 0):>6,} ({by_class.get(7, 0) / valid:.3%})" )

default (EPSG:4326)      valid_pixels=  320,460  class 7= 2,049 (0.639%)

equal-area (+proj=cea)   valid_pixels=  296,652  class 7= 1,878 (0.633%)

native UTM (EPSG:32610)  valid_pixels=  326,448  class 7= 2,115 (0.648%)

The counts move by roughly 9% across those three projections; the percentages hold to within about 2%.

Two rules follow:

• Report percentages, not raw pixel counts, or state the projection alongside the count.
• If you need area in km2, ask for an equal-area projection such as +proj=cea. In geographic coordinates a "square" pixel covers less ground east-to-west than north-to-south, and the discrepancy grows with latitude.