Zonal-stats speed benchmarks

Archival chunk: stats compute per country

100 SEAS5 issuances × 7 leadtimes (700 COG slices), admin levels 0–2, full pipeline stats path incl. weighted median and row validation. Seconds of compute per country on Databricks job compute — excludes the ~62 s/country remote data fetch, which is identical in both methods (see notes). Laptop figures are in the table view.

Cost vs number of admin units

Core zonal kernel on a Nigeria-scale grid (25×30 native pixels at 0.4°, 8× upsample equivalent), 100 slices, synthetic contiguous admin regions; best of 3 runs after warm-up, on Databricks job compute. Both axes log — the legacy cost climbs linearly with admin count while the weight-matrix cost stays well below it at every size. This is what makes adding many more (smaller) admin units cheap. Laptop figures are in the table view; the ratio narrows at large admin counts on the cluster node (11× at 3,000 vs 23× on the laptop) because the sparse-matrix path is the more memory-bandwidth-sensitive of the two.

Both axes at once: admin units × dates

The two sweeps above each vary one axis. This walks the full grid, so the interaction is visible. Cell colour and label are the speedup of PR #49 over the legacy method, on Databricks job compute, running the real pipeline runners (weighted median and row validation included) over synthetic non-pixel-aligned polygons on a Nigeria-sized SEAS5 grid.

The two axes behave completely differently. Legacy cost is a clean product of both — 3,000 admins over 200 dates costs it 749 s. The new method is close to flat in dates: at 3,000 admins, going from 1 date to 200 moves it only 0.61 s → 0.97 s, because the sparse matrix product over 200 slices costs barely more than over one. So the multiplier comes from the date axis, while the admin axis mostly sets the floor (5–8× at a single date warm, and see the cold table below for the worst case where the weight build is not yet amortised). The practical reading: daily update runs sit in the cheap corner either way, and adding admin units is close to free — 800 → 3,000 admins costs the new method 0.25 s → 0.97 s at 200 dates, where it costs the legacy method 199 s → 749 s.

One asymmetry in the grid above, and the worst case it hides. The legacy timings rebuild their rasterized mask on every repetition, because that is what the legacy code does. The new-method timings load the weight matrix from its disk cache, because that is what the pipeline does after the first chunk touches a country. That is a fair model of a production run, but it amortises the weight build to zero — so the table below re-measures the same grid with the build inside the timed call (cold), which is the honest worst case for a single isolated run.

Two things follow. First, the weight build is a fixed cost per country and admin level, not per date — cold minus warm is 0.021 / 0.018 / 0.017 s at 10 admins across 1, 10 and 50 dates, and 0.659 / 0.647 / 0.648 s at 3,000. It scales with admin units (it is geometry work) and is then cached across every subsequent chunk, country-worker and run. Second, even paying it in full the new method is never slower: the true worst case in this grid — a single date over only 10 admin units, where there is nothing whatsoever to amortise against — is 1.04×, essentially break-even, rising to 2.9× at 3,000 units and past 100× once there are 50 dates to spread it over. Separate run from the grid above, so its warm figures differ slightly (different cluster instance); the cold/warm comparison within it is self-consistent.

Update run: stats compute per country

One SEAS5 issuance (7 leadtimes), admin levels 0–2, warm weight cache, on Databricks job compute. Even at this small scale the method change is worth 9–13× per country. The rest of a production update run is boundary downloads, DB writes, and serial country processing, which PR #49 addresses separately (per-run boundary cache, configurable workers, country-splitting for single-chunk runs).

Reproducing the 2024 benchmark

The 2024 methodology work benchmarked baseline compute_zonal_statistics (python-rasterstats per date) against fast_compute_zonal_statistics (rasterize-once + numpy — today's legacy method). This reruns the original functions from the 2024 exact-extract branch, same setup as far as it is recoverable — ERA5 monthly at native 0.25°, ETH and PHL ADM2 CODs, sweeps over N dates — with the PR #49 kernel as a third line. Outputs of baseline and fast were verified identical before timing. Log y-axis; seconds per run (best of 10 iterations; 3 for the baseline).

What reproduces: every 2024 shape and ordering — all methods linear in dates, baseline linear in polygons and fast near-flat (ETH), and PHL's baseline far slower than ETH's despite a smaller grid (its coastline multipolygons dominate rasterstats' per-polygon cost). Absolute times differ from the 2024 charts (different machine and library versions): this run's baseline is slower and its fast is faster, so the baseline→fast gap measures 59–428× here vs ~9–175× in 2024. One shape did not reproduce: the 2024 PHL polygons chart was nearly flat (1.42→1.65 s over 1→88 polygons) where this run grows linearly (0.07→12.0 s) — consistent with their PHL polygon sweep having run on the upsampled 0.05° grid, where a large fixed per-date cost swamps the per-polygon term; the sweep configuration was never committed, so this is a hypothesis. At the largest size measured here the 2024-fast → PR #49 gap is 5.7× (ETH) and 9.8× (PHL) — and note that both kernels are running on the native grid, faithful to the 2024 harness, whereas the production legacy path first upsamples to 0.05° (25× the pixels for ERA5). Rerunning the same ETH/PHL data with the legacy method configured as production actually runs it widens the gap to 16–17× (ETH 0.7 s vs 0.043 s; PHL 0.5 s vs 0.030 s at 125 dates), and it grows further with slices and admin count — see the pipeline-scale and 2-D sections above.

Correction (2026-08-27). The first version of this figure timed with the mean of repeated runs, which absorbs contention from anything else on the machine rather than rejecting it. Other benchmarks were running concurrently, and because the largest point in each sweep is measured last it took the worst of it: the PR #49 figure at 125 dates was published as 0.076 s (PHL) and 0.042 s (ETH) against true values of 0.014 s and 0.027 s — inflating PHL roughly 6×, and making the gap look like 2–4× when it is 5.7–9.8×. Re-measuring the kernel in isolation on a quiet machine confirmed there is no cliff: per-slice cost is flat at ~0.10 ms from 75 through 200 slices. All timings on this page now use the minimum of repeated runs, which estimates the uncontended cost. The Databricks figures elsewhere on this page were unaffected — a dedicated cluster with nothing else on it, already min-of-N.