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, in the order you supplied them. Stencil and operator matrix rows follow that same order, including cache hits.

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 global grids, set period=360 on resample_grid, reduce, or the resampler/fused-operator builder to interpolate across the longitude seam. Resolution-based calls then generate one full longitude cycle; explicit target arrays remain unchanged. The default period=None still clamps at the edges. Do not repeat the source endpoint (e.g. include −180° but not also +180°). This wraps sampling, not polygon geometries. See periodic longitude.

For area-weighted coarsening or bounded fields, use resample_grid(da, target_resolution=2, method="conservative"), or prebuild a conservative Resampler for explicit target centres/bounds. It applies spherical cell overlaps with two sparse 1-D factors; no full grid-to-grid matrix is built. The existing method="meanpreserving" remains the default and is intended for smooth refinement. See conservative regridding for footprint alignment, partial coverage, and opt-in NaN handling.

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.

Already have your own reader? Pass its NumPy window arrays, in plan order, to restricted.gather(arrays), then restricted.apply(gathered, how="mean"). No xarray objects or Dask scheduling are needed for this path.

Prebuilt reducers and restricted.apply accept skipna=True to omit missing source cells and renormalize means. NaNs still propagate by default. When resampling is fused, source-cell masking differs from the stencil path's resample-then-mask behavior; see NaN handling.

For float32 grids, opt into smaller coefficients and results independently:

stencil = ghl.Stencil.compute(lats, lons, geoms, dtype=np.float32)
op = ghl.ReduceOperator.compute(stencil, lats, lons, dtype=np.float32)
out = ghl.reduce_with_operator(da.astype(np.float32), op, preserve_dtype=True)

Builders default to float64. All polygon reducers and the operator's NumPy methods accept preserve_dtype=True; floating inputs retain their result dtype, while integer means remain floating-point. Normalizers stay float64, and matrix precision is not lowered at apply time. See dtype controls.

Both resampling methods also accept preserve_dtype=True on resample_grid, resample_grid_with_matrix, and Resampler.apply_grid. Each real floating input retains its result dtype without lowering computation precision; float32 output arrays use half the storage of float64. See resampling result dtypes.

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 · portable NPZ artifacts
• 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.