geohalo

Exact-fractional-area zonal statistics over regular lat/lon grids.

Given a regular lat/lon mesh of gridded values — temperature, precipitation, population density, a land-cover fraction, a satellite band, … (loaded with xarray from GRIB, NetCDF, Zarr, …) — and an arbitrary set of polygons, geohalo reduces the spatial dimensions of the mesh to one value per polygon with sub-cell precision and millisecond-scale aggregation in the hot path.

The expensive geometric work happens once; every subsequent grid collapses to a single sparse · dense matmul.

📖 Full documentation: https://campiohe.github.io/geohalo/

How it works

Aggregation is a linear operator:

aggregates = W @ flat_grid_values

where W ∈ ℝ^(N_polygons × N_cells) is a sparse matrix whose entries are the exact fractional area of cell ∩ polygon weighted by each cell's true surface area on a sphere. W (the Stencil) depends only on the grid topology and the polygon set — not on the grid values — so it is built once (and cacheable) and reused for every slice. See Aggregation as a linear operator.

Install

geohalo targets Python ≥ 3.12.

uv add geohalo            # or: pip install geohalo

Optional extras: redis (the RedisCache backend) and matplotlib (the helpers in geohalo.plot).

Quickstart

import numpy as np
import geopandas as gpd
import xarray as xr
from shapely.geometry import box
import geohalo as ghl

# any regular lat/lon DataArray works; a synthetic field so this runs as-is
lats = np.arange(-25.0, -19.0, 0.25)
lons = np.arange(-50.0, -42.0, 0.25)
lon2d, lat2d = np.meshgrid(lons, lats)
field = 290.0 + 5.0 * np.cos(np.deg2rad(4 * lat2d)) + 0.1 * lon2d

da = xr.DataArray(
    field, dims=("latitude", "longitude"),
    coords={"latitude": lats, "longitude": lons}, name="value",
)

geoms = gpd.GeoSeries(
    [box(-49, -24, -47, -22), box(-47, -24, -45, -22), box(-46, -22, -44, -20)],
    index=["SP", "RJ", "MG"],            # the index holds the keys
)

out = ghl.reduce(da, geoms)                  # hot path; ms-scale
out_fine = ghl.reduce(da, geoms, target_resolution=0.05)   # refine the grid first
# out: xr.DataArray over (..., geom)

The output preserves every non-spatial dim of da (time, ensemble member, band, vertical level, …) and replaces (latitude, longitude) with a single geom dim indexed by the GeoSeries keys.

reduce also accepts an xr.Dataset (every spatial data var is reduced), how={"mean", "sum"}, a weight_key naming a per-cell weight variable, and spherical_correction=False to disable the latitude-area correction.

For clean lazy grids, opt into chunk-aware reads with a precomputed fused operator:

restricted = ghl.RestrictedOperator.from_grid(op, da)
out = ghl.reduce_with_restricted_operator(da, restricted)

This loads only the contributing spatial chunks and processes batch chunks separately. See chunk-aware reduction for explicit chunk sizes, caching, and backend limitations.

Documentation

Everything is covered in depth at https://campiohe.github.io/geohalo/:

• Quickstart — runnable examples: batching, datasets, weights, NaN handling, refining, rollups.
• Concepts — linear operator · the stencil · why exact fractional coverage · latitude correction · mean-preserving downscaling · the fused reduce operator · NaN-aware & weighted reduction · hierarchical rollups
• Guides — caching the precompute · resampling grids
• API reference

Performance

The hot path is a single sparse · dense matmul: a 50-member batch over the ~5,570 GADM Brazil L2 municipalities reduces in single-digit milliseconds, and the one-time Stencil precompute is seconds and cacheable. Methodology, full tables, and the fused-operator size win are on the Performance page; re-run the suite with uv run python -m benchmarks.run.

Non-goals

• No reprojection — EPSG:4326 throughout (grids and polygons).
• No per-variable cache — the Stencil depends on grid + polygons only.
• No WGS84-ellipsoidal cell areas — spherical is within ~0.3 % (spherical_correction=False gives planar/equal-area weights).
• No DAG hierarchies — each child has exactly one parent (tree only).
• No how={"min", "max"} — mean and sum only.

Development

uv sync                                          # install deps
uv run pytest                                    # tests
uv run ruff check .                              # lint
uv run --group docs mkdocs serve                 # preview the docs locally
uv run --group docs python docs/gen_figures.py   # regenerate the doc figures

Docs are built with Material for MkDocs and deployed to https://campiohe.github.io/geohalo/ on every push to main.