# climate_indices

Reference implementations of climate and drought indices for scientific monitoring workflows.

Generated by `uv run scripts/generate_llms_txt.py` from committed repository documentation.

## README.md

![Banner Image](https://raw.githubusercontent.com/monocongo/climate_indices/main/assets/multi-index-compare.png)

# climate_indices

[//]: # ([![Coverage Status]&#40;https://coveralls.io/repos/github/monocongo/climate_indices/badge.svg?branch=main&#41;]&#40;https://coveralls.io/github/monocongo/climate_indices?branch=main&#41;)
[//]: # ([![Codacy Status]&#40;https://api.codacy.com/project/badge/Grade/48563cbc37504fc6aa72100370e71f58&#41;]&#40;https://www.codacy.com/app/monocongo/climate_indices?utm_source=github.com&amp;utm_medium=referral&amp;utm_content=monocongo/climate_indices&amp;utm_campaign=Badge_Grade&#41;)
[![Actions Status](https://github.com/monocongo/climate_indices/workflows/tests/badge.svg)](https://github.com/monocongo/climate_indices/actions)
[![License](https://img.shields.io/badge/License-BSD%203--Clause-green.svg)](https://opensource.org/licenses/BSD-3-Clause)
[![Python | 3.10-3.14](https://img.shields.io/badge/Python-3.10--3.14-blue?logo=python)](#supported-python-versions)

#### Python library of climate indices for drought, wildfire, and flood monitoring

`climate_indices` provides reference implementations of climate index algorithms that give a
geographical and temporal picture of the severity and duration of precipitation, temperature,
evaporative-demand, fire-weather, and wet-extreme anomalies for climate monitoring and research.
It began as a drought library; it now covers drought, wildfire, and flood potential, and an
energy family is in development.

### Indices

**Drought and moisture balance**

- [SPI](https://climatedataguide.ucar.edu/climate-data/standardized-precipitation-index-spi),
  Standardized Precipitation Index, utilizing gamma and Pearson Type III distributions, with
  selectable placement of zero accumulations (`zero_handling`)
- [SPEI](https://climatedataguide.ucar.edu/climate-data/standardized-precipitation-evapotranspiration-index-spei),
  Standardized Precipitation Evapotranspiration Index, utilizing gamma, Pearson Type III, and generalized-logistic (log-logistic) distributions
- Standardized index (`indices.standardized_index()`), the SPI fitting pipeline for any
  non-negative series, such as runoff or streamflow for
  [SRI/SSI work](docs/standardized-hydrologic-indices.md)
- [PNP](https://www.droughtmanagement.info/percent-of-normal-precipitation/),
  Percentage of Normal Precipitation
- [PCI](https://doi.org/10.1111/j.0033-0124.1980.00300.x), Precipitation Concentration Index
- [EDDI](https://psl.noaa.gov/eddi/), Evaporative Demand Drought Index
- [Palmer indices](https://www.droughtmanagement.info/literature/USWB_Meteorological_Drought_1965.pdf),
  including PDSI, PHDI, PMDI, Z-Index, and [scPDSI](https://doi.org/10.1175/1520-0442(2004)017%3C2335:ASPDSI%3E2.0.CO;2)

**Evapotranspiration**

- [PET](https://www.ncei.noaa.gov/access/monitoring/dyk/potential-evapotranspiration),
  Potential Evapotranspiration, utilizing the [Thornthwaite](https://doi.org/10.2307/210739),
  [Hargreaves](https://doi.org/10.13031/2013.26773), or
  [FAO-56 Penman-Monteith](https://www.fao.org/4/x0490e/x0490e00.htm) equations

**[Wildfire](https://climate-indices.readthedocs.io/en/latest/wildfire_applications.html)** (`climate_indices.fire`)

- [KBDI](https://www.drought.gov/data-maps-tools/keetch-byram-drought-index), the Keetch-Byram Drought Index
- [CFFWIS](https://natural-resources.canada.ca/forests-forestry/wildland-fires/canada-fire-weather-index-system),
  the Canadian Forest Fire Weather Index System: the FFMC, DMC, and DC moisture
  codes (with Drought Code overwintering) and the ISI, BUI, FWI, and DSR behavior indices
- [Fosberg Fire Weather Index (FFWI)](https://wikifire.wsl.ch/tiki-indexe575.html?page=Fosberg+fire+weather+index)
- [Hot-Dry-Windy Index (HDW)](https://research.fs.usda.gov/treesearch/56562)
- [Haines Index](https://forecast.weather.gov/glossary.php?word=haines+index)

**[Flood potential](https://climate-indices.readthedocs.io/en/latest/flood_applications.html)** (`climate_indices.flood`)

- [Effective Precipitation (PE)](https://digitalcommons.unl.edu/droughtfacpub/32/) and its derived
  Effective Drought Index (EDI, distinct from EDDI), the
  [Flood Index (I_F)](docs/algorithm_refs/flood_index.md), and the
  [Antecedent Precipitation Index (API)](https://forecast.weather.gov/glossary.php?word=Antedecent+Precipitation+Index):
  daily accumulated-wetness measures that indicate flood potential, not flooding
- `flood_events()`, which extracts events with onset, duration, peak, and severity from a
  daily flood index

**Energy (in development, [#1157](https://github.com/monocongo/climate_indices/issues/1157))**

Planned, none shipped yet in the package:

- Heating, cooling, and energy degree days (HDD / CDD / EDD)
- Hub-height wind-speed profile, and wind and temperature-derated PV capacity factors
- Standardized energy production and demand indices: SREPI and a standardized demand index,
  plus SRLI for signed residual-load series, following
  [Allen and Otero (2023)](https://doi.org/10.1016/j.renene.2023.119206)
- Run-theory energy-drought events (onset, duration, severity), and a capacity-weighted
  precipitation proxy for hydropower
- These indices describe the climate-driven component of energy supply and demand, not
  actual generation or load

**Supporting tools**

- Run theory (`climate_indices.runs`): runs above or below a threshold with duration,
  magnitude, intensity, and peak, per series or per grid cell
- Distribution fit diagnostics (`fit_diagnostics()`): the fitted parameters, probability of
  zero, and Kolmogorov-Smirnov statistic per calendar step and cell
- Probability-scale output (`output_scale`): standardized indices as normal scores, fitted
  cumulative probabilities, or a bounded `[-1, 1]` scale

### Highlights

- **Rust backend.** An optional compiled extension (`climate_indices._native`)
  runs Rust ports of the numerical kernels behind SPI, SPEI, the standardized index, EDDI,
  PET, PNP, PCI, the Palmer family, the KBDI and CFFWIS moisture-code recurrences, and the
  flood family. It preserves the public API; a cross-backend parity registry checks every
  ported kernel against the Python reference at `rtol = atol = 1e-10` with matching NaN
  positions, not bit-for-bit equality; see the [parity contract](docs/architecture.md#parity-tolerance).
  Speedups vary by kernel: the flood family
  measured 1.2x to 6.7x faster than its Python path and Thornthwaite PET 1.4x, while
  Hargreaves measured at parity and Penman-Monteith slightly slower (0.89x) at the benchmark
  size; see [the measurements](benchmarks/README.md).
- **Vectorized gridded computation.** SPI, SPEI, EDDI, PNP, PET, and PDSI fit or rank a
  whole `(time, *cells)` block in one call instead of looping over cells. On a 38x87-cell,
  40-year reference grid that is 4.4x faster for SPI, 3.6x for SPEI, 53x for Thornthwaite
  PET, 342x for EDDI, and 111x for gridded PDSI. On fully populated monthly grids, gamma
  SPI/SPEI, EDDI, and PNP are bit-identical to per-cell computation; Thornthwaite PET and
  successful Pearson fits agree within `1e-12`. Pearson SPI's block-wide gamma fallback is
  not per-cell equivalent and can depend on chunk layout; daily and missing-data bounds
  are [documented separately](docs/performance.md#numerical-equivalence).
- **xarray, Dask, and Zarr.** Supported xarray entry points accept `xarray.DataArray` input,
  preserve coordinates, attach CF metadata, and stay lazy on Dask-backed arrays for
  chunk-by-chunk Zarr output. Temporal inference and core-dimension chunk requirements
  vary by function; `indices.standardized_index()` and scPDSI remain NumPy-only.
  See the [compatibility matrix](docs/xarray_compatibility.md),
  [writing NetCDF and Zarr outputs](docs/writing-outputs.md), the
  [Zarr/Dask notebook](notebooks/zarr_dask_spi_spei.ipynb), and
  [chunking and scheduler guidance](docs/performance.md).
- **Validation.** Documented evidence distinguishes external checks against NOAA, NRCan,
  SPEIbase, and literature reference data from regression-only coverage; see
  [Validation](#validation).
- **Command-line interface.** `climate_indices --index
  {spi,spei,pnp,scaled,pet,palmers,kbdi,pe,edi,flood_index,api,all}` computes indices from
  NetCDF inputs.

### Project goals

This Python implementation of the above climate index algorithms is being developed
with the following goals in mind:

- to provide an open source software package to compute a suite of
  climate indices commonly used for climate monitoring, with well
  documented code that is faithful to the relevant literature and
  which produces scientifically verifiable results
- to provide a central, open location for participation and collaboration
  for researchers, developers, and users of climate indices
- to facilitate standardization and consensus on best-of-breed
  climate index algorithms and corresponding compliant implementations in Python
- to provide transparency into the operational code used for climate
  monitoring activities at NCEI/NOAA, and consequent reproducibility
  of published datasets computed from this package
- to incorporate modern software engineering principles and scientific programming
  best practices


This is a developmental/forked version of code that was originally developed by NIDIS/NCEI/NOAA.
See [drought.gov](https://www.drought.gov/drought/python-climate-indices).

- [__Documentation__](https://climate-indices.readthedocs.io/en/latest/)
- [Climate Indices for Wildfire Applications](https://climate-indices.readthedocs.io/en/latest/wildfire_applications.html)
- [Climate Indices for Flood and Wet-Extreme Applications](https://climate-indices.readthedocs.io/en/latest/flood_applications.html)
- [__Validation status__](VALIDATION.md)
- [__Changelog__](CHANGELOG.md)
- [__License__](https://github.com/monocongo/climate_indices/blob/main/LICENSE)
- [__Disclaimer__](https://github.com/monocongo/climate_indices/blob/main/DISCLAIMER)

## Installation

This README describes 3.0.0, which is not yet published on PyPI. The released 2.x
package does not include the `fire`, `flood`, or `runs` modules or the Rust backend.
For 3.0.0 features, install from a [development checkout](docs/development-guide.md#installation).

Install the released 2.x package from PyPI:

```bash
pip install "climate-indices>=2.3"
```

`uv` users can run `uv pip install "climate-indices>=2.3"`. The 2.3.0 floor matches the
[quickstart](https://github.com/monocongo/climate_indices/blob/main/docs/quickstart.md),
which uses the xarray API added in 2.3.0; see
[Supported Python Versions](#supported-python-versions) for interpreter support.

Starting with 3.0.0, wheels for Linux (x86-64 and aarch64, manylinux_2_28, so glibc 2.28 or
newer), macOS (Apple silicon and Intel), and Windows (x86-64) include the optional Rust
acceleration backend; every other platform, including a Linux system with older glibc,
installs the pure-Python wheel. Both carry the same version and the same public
API; results agree within the [cross-backend parity tolerance](docs/architecture.md#parity-tolerance),
not necessarily bit for bit. Calling code need not depend on which one pip selected.
To see which one you have:

```bash
python -c "import climate_indices._native"   # succeeds when the extension is installed
```

A missing extension is not an error: `climate_indices` imports, and every computation
runs the Python implementation. Installing from source needs no Rust toolchain either —
`pip install --no-binary climate_indices climate_indices` installs from the sdist, which
builds the pure-Python wheel. Developers who
want the extension from a checkout build it with `uv run maturin develop --release`; see
[the Rust backend section of the architecture notes](docs/architecture.md#optional-rust-backend).

## Developer Workflow

This project uses trunk-based development. `main` is the trunk and should always
be releasable.

1. Start from current trunk:
   `git switch main && git pull --ff-only origin main`
2. Create a short-lived branch:
   `git switch -c feature/<short-topic>`
3. Make focused changes with tests.
4. Run validation:
   `uv run ruff check src/ tests/`
   `uv run ruff format --check src/ tests/`
   `uv run mypy src/ tests/test_type_checking.py`
   `uv run pytest`
5. Open a PR into `main`.
6. The maintainer merges the PR after review and passing CI; agents never merge.

Use `feature/<topic>`, `fix/<topic>`, `docs/<topic>`, `chore/<topic>`,
`perf/<topic>`, `refactor/<topic>`, `test/<topic>`, `ci/<topic>`, or
`hotfix/<topic>` branch names. Release branches are avoided; use maintenance
branches only for approved older-version support.

The methodology behind this workflow — parallel agent sessions, worktree
isolation, and the evidence each change must carry — is documented in
[AI-assisted development](docs/ai-assisted-development.md).

## Release Recipe

Releases are tag-based. The Git tag, package version, GitHub Release, and PyPI
version must match.

- Git tag: `v1.2.3`
- Package version: `1.2.3`
- GitHub Release: `v1.2.3`
- PyPI release: `1.2.3`

1. Prepare a release PR that updates `pyproject.toml`, `CHANGELOG.md`,
   and release notes/docs; the maintainer merges it after review and passing
   CI.
2. Confirm `main` is green.
3. Create an annotated tag from `main`.
4. Push the tag. The release workflow builds, validates, publishes to PyPI, and
   creates the GitHub Release.

Tag creation and publishing require maintainer approval. See
[`docs/release-process.md`](docs/release-process.md) for the full checklist.

### Maintainer Quick Commands

Read-only preflight:

```bash
git status --short
git branch --show-current
git log --oneline --decorate -5
uv run pytest tests/test_release_integrity.py
```

Safe PR branch setup:

```bash
git switch main
git pull --ff-only origin main
git switch -c chore/issue-667-release-docs
```

Approval-required release tag commands:

```bash
git switch main
git pull --ff-only origin main
git tag -a vX.Y.Z -m "Release vX.Y.Z"
git push origin vX.Y.Z
```

## Supported Python Versions

Supports Python **3.10–3.14**.

All supported versions are tested on Linux (ubuntu-latest) for every pull request and push
to `main`. Python 3.14 is additionally tested on macOS on every event, and 3.10 on macOS on
pushes to `main`, merge groups, the weekly schedule, and manual dispatch. Both latest and
minimum declared dependency versions are tested in CI.

### Version Support Policy

This project provides **12 months notice** before dropping support for a Python version.
When a version approaches end-of-life, removal will be announced via the CHANGELOG and a
GitHub issue, and implemented no sooner than 12 months after announcement with a version bump.

Python 3.9 support was dropped in v2.2.0 (August 2025) due to `scipy>=1.15.3` requiring 3.10+.

### API Stability

| API Surface | Status | Guarantee |
|:------------|:------:|:----------|
| NumPy array functions (`indices.spi`, `indices.spei`, `indices.pet`, and the `flood` indices) | **Stable** | No breaking changes in minor versions |
| xarray DataArray entry points (any index function accepting `xr.DataArray`) | **Beta** | No breaking changes in patch versions |

**Stable API**: The NumPy-based computation functions follow strict semantic versioning.

**Beta API**: The xarray adapter layer provides automatic parameter inference, coordinate
preservation, CF metadata, and Dask support. While beta, computation results follow the
stable NumPy API's documented [numerical-equivalence bounds](docs/performance.md#numerical-equivalence)
— only the interface surface (parameter names, metadata attributes,
coordinate handling) may evolve. Beta features are marked in docstrings, with
``BetaFeatureWarning`` as their public warning category. The adapter stays Beta in 3.0.0
and is promoted no earlier than 3.1.0, once the 3.0.0 calendar alignment has a release
of soak time ([ADR-0012](docs/adr/0012-xarray-api-stays-beta-through-3.0.0.md)).

See `docs/xarray_compatibility.md` for the 3.0.0 compatibility matrix, including
Dask chunking constraints, metadata behavior, and the Palmer xarray adapter.
See `docs/performance.md` for a runnable parallel SPI/SPEI example and
the measured speedups behind the chunk and scheduler guidance.

## Validation

[`VALIDATION.md`](VALIDATION.md) classifies the evidence for the indices in its per-index
inventory, separating independent external validation from regression coverage and recording
their known gaps. Fosberg Fire Weather Index and Haines Index are not yet classified there.
`uv run pytest -m validation` runs the marked reference-fixture and specification checks.
FAO-56 worked-example checks run in the core suite via `uv run pytest tests/test_pm_eto.py`.
The evidence includes:

- **EDDI**: NOAA PSL monthly reference ET/EDDI for 1-, 3-, and 6-month timescales
  (1979–2023); maximum observed error `2.43e-6`.
- **SPI**: characterized against NOAA NCEI's operational climate-divisional SPI for all 344
  divisions and 7 timescales; zero placement compared with the R `SEI` and `SCI` packages.
- **SPEI**: reproduces SPEIbase v2.11 from its own CRU TS inputs (correlation
  0.9996–0.99999) and matches the R `SPEI` package's log-logistic fit within `1e-6`.
- **Standard Palmer** (PDSI, PHDI, PMDI, Z-Index): qualified external validation against
  NOAA NCEI nClimDiv for all 344 divisions, 1895–2022 (median absolute difference 0.0127).
- **scPDSI**: cross-validated against Wells-lineage reference fixtures (`atol=5e-5`).
- **CFFWIS**: all seven components match the NRCan `cffdrs` reference within `1e-9`.
- **KBDI**: reproduces the Keetch & Byram (1968) worked example.
- **HDW**: reproduces event timing and tracks the daily series in the Srock et al. (2018)
  Cedar Fire case; published magnitudes are characterization only, not numerical validation.
- **PET**: selected FAO-56 helper examples and the daily end-to-end Example 18, plus
  Thornthwaite and Hargreaves literature examples; monthly Example 17 is not encoded.
- **Flood family**: specification-level and regression coverage only; no external numeric
  oracle has been adopted yet.

The Rust backend has its own gate: the parity registry in `tests/parity_registry.py` runs
each ported entry point through the extension and through the Python implementation on
property-generated inputs, and fails if a kernel is added without an entry.

## Upgrading to 3.0.0

3.0.0 changes some computed values and exception types without a deprecation period, most
notably corrected daily calendar alignment in the xarray adapter, the PCI February fix, and
the gamma probability of zero now taken from the Calibration Period. It also removes the `spi`
console script in favor of `climate_indices --index spi`. Each change, how to detect it, and
what to change is in [`docs/deprecations/api-changes.md`](docs/deprecations/api-changes.md);
the full list is in the [CHANGELOG](CHANGELOG.md). Upgrading from 2.1.x or earlier also
crosses the 2.2.0 switch from `None` return tuples to exceptions.

## Citation
You can cite `climate_indices` in your projects and research papers via the BibTeX
entry below.
```
@misc {climate_indices,
    author = "James Adams",
    title  = "climate_indices, an open source Python library providing reference implementations of commonly used climate indices",
    url    = "https://github.com/monocongo/climate_indices",
    month  = "may",
    year   = "2017--"
}
```

## VALIDATION.md

# Validation Status

This document records the validation state for the 3.0.0 release artifacts. It
separates external scientific validation from regression coverage so users can
see which results are backed by independent reference data and which are still
guarded only by internal fixtures.

## Test Commands

| Scope | Command | Expected result |
| --- | --- | --- |
| Core suite | `uv run pytest -m "not benchmark and not validation"` | Unit, property, xarray, release guardrail, and regression tests pass. |
| Validation marker suite | `uv run pytest -m validation` | External validation tests, including committed NOAA EDDI fixtures, the SPEIbase v2.11 SPEI plausibility floors and input-matched log-logistic comparison, the R `SPEI` log-logistic cross-implementation fixtures, the standard-Palmer nClimDiv comparison, the Srock et al. (2018) HDW event-discrimination fixture, and the flood-family source-identity checks (specification-level, not external; see "Flood Evidence Classification"), pass; scPDSI oracle cross-validation passes. |
| Lint | `uv run ruff check src/ tests/` | No lint findings. |
| Format | `uv run ruff format --check src/ tests/` | No formatting changes needed. |
| Notebooks | `uv run jupyter nbconvert --execute --to notebook --inplace notebooks/xarray_getting_started.ipynb notebooks/palmer_indices_xarray.ipynb notebooks/eddi_xarray.ipynb` | All 3.0.0 notebooks execute from a clean checkout. |
| Notebook smoke (Zarr/Dask tutorial, local only) | `bash scripts/smoke_e2e_notebook.sh` | `notebooks/zarr_dask_spi_spei.ipynb` executes end to end from a fresh kernel against the prepared sample inputs. |

## Per-Index Evidence

| Index | Status | Tolerance | Evidence | Known gap |
| --- | --- | --- | --- | --- |
| SPI | Validated (internal) + qualified external comparison (Pearson III) | Gamma `atol=1e-8`, Pearson `atol=1e-5` in xarray equivalence tests; NOAA NCEI climdiv characterization ceilings documented in `tests/test_ncei_spi_reference.py` | Legacy NumPy fixtures and xarray equivalence tests compare gamma and Pearson outputs across scales. `tests/test_ncei_spi_reference.py` characterizes `indices.spi()` (Pearson III, full-period-of-record calibration) against NOAA NCEI's operational climdiv SPI product for all 344 climate divisions and all 7 published timescales; marked `validation`. | `tests/fixture/ncei_spi/provenance.json` records that NCEI's climate-divisional SPI product's independence from this codebase's lineage could not be confirmed as airtight from public documentation alone (see `docs/research/spi-dataset-survey.md`), and that NCEI's actual calibration behavior is a full/expanding period-of-record window rather than its documented fixed 1931-1990 baseline (Baldwin & Chen 2020). WMO SPI User Guide worked-example extraction remains outstanding (issue #778). |
| SPI zero placement | SEI independent-implementation comparison; SCI estimator cross-check; closed-form and property regression | Zero-placement constants `atol=1e-12`; parameter-matched full series `atol=1e-6`; end-to-end `atol=0.01`; SCI centre-of-mass `atol=0.2` (Weibull-vs-empirical `p0` estimator difference) | `tests/test_zero_handling.py` pins each ADR-0015 mode's normal-scale closed form at `p0 ∈ {0.1, 0.5, 0.9}`, the ideal sample means, and the zero/positive ordering. `tests/test_zero_handling_reference.py` (marked `validation`) compares the modes to committed fixtures generated by R `SEI` 0.2.0 (`std_index(..., lower = 0, cens = "prob" | "normal")`) and R `SCI` 1.0-3 (`p0.center.mass = TRUE`): the SEI zero constants match exactly, the full series matches when run on SEI's fitted shape/rate, and the end-to-end difference is at most 0.0032 (SEI's MLE gamma versus this library's method-of-moments approximation). SCI's committed `p0` estimator is independently produced and its documented `(np+1)/(2(n+1))` placement is compared to the exercised mode. The fixture also asserts its SEI shape/rate differ from this library's own fit, so a self-derived fixture cannot pass. Fixture provenance and tolerances are in `tests/fixture/zero_handling/provenance.json`; regenerate with `scripts/prepare_zero_handling_fixtures.py`. | An all-zero calibration step diverges by design: SEI's censored modes assign their constant without a fit, while this library follows ADR-0015 decision 5 and scores the step's zeros as classic (SEI classic gives `+Inf`). SCI estimates `p0` with the Weibull plotting position `np/(n+1)` and places zeros at `(np+1)/(2(n+1))`, so its constants differ from this library's `np/(2n)` by at most 0.16 in z for this fixture and more for sparser steps; `sci_center_of_mass_probability` records the estimator it uses. |
| SPEI | Validated (internal) + input-matched cross-implementation comparison (log-logistic, SPEIbase v2.11); gamma/Thornthwaite plausibility check | Gamma `atol=1e-8`, Pearson `atol=1e-5` in xarray equivalence tests; log-logistic like-for-like floors (correlation >= 0.97) in `tests/fixture/speibase_cru_ts/provenance.json`; gamma plausibility floors in `tests/fixture/speibase/provenance.json` | Legacy NumPy fixtures and xarray equivalence tests compare gamma and Pearson outputs across scales. `tests/test_speibase_like_for_like.py` feeds `indices.spei()` the CRU TS 4.09 precipitation and FAO-56 Penman-Monteith PET that SPEIbase v2.11 was computed from, standardizes each 0.5-degree cell with the log-logistic distribution over 1901-2024, averages the per-cell series inside three CONUS climate divisions, and reproduces the CSIC SPEIbase v2.11 grids (correlation 0.9996-0.99999, mean absolute difference 0.003-0.016 across divisions and timescales). `tests/test_speibase_reference.py` retains the earlier gamma/Thornthwaite comparison against grid-cell-mean SPEIbase as a plausibility check. Both are marked `validation`. | The like-for-like comparison removes the PET-method, distribution-family, and combined precipitation/spatial-support confounds of the plausibility check. `climate_indices.lmoments.fit_glo` is a port of the same `lmom` PELGLO routine R `SPEI` uses, so close agreement is expected, not algorithm independence. The small residual is dominated by a PET day-length convention difference: the test uses leap-aware month lengths (matching public SPEIbase `R/functions.R`), while a near-uniform month reproduces the committed v2.11 grids two to four times more closely. The gamma/Thornthwaite comparison remains plausibility only; see "SPEI Plausibility Classification" below. |
| SPEI log-logistic | Validated (external independent implementation) | Series `atol=1e-6`; GLO parameters `rtol=1e-8`, `atol=1e-6` | `tests/test_spei_loglogistic_reference.py` (marked `validation`) compares `indices.spei()` (`Distribution.loglogistic`) against the R `SPEI` package (`distribution = "log-Logistic"`, `fit = "ub-pwm"`) for the same three CONUS climate divisions at timescales 1/3/6/12, plus a short highly skewed synthetic series whose constant calendar step is unfittable at scale 1. Both sides fit the identical `(P − PET) + 1000` scaled series. | R `SPEI` fits each calendar step with unbiased-PWM L-moments (`TLMoments::PWM` + `lmom::pelglo`), the estimator `lmoments.fit_glo` implements, so the standardized series agree at ~1e-12 and the parameters to ≤~1e-10; the recorded tolerance absorbs cross-platform wobble only. `SPEI` does not clip to `[-3.09, 3.09]`, so the tests exclude the clipped tail and assert the clip separately. See "SPEI Log-Logistic Cross-Implementation" below. |
| PET Thornthwaite | Validated | Synthetic-fixture `atol=0.001`; literature worked example `atol=4.0` mm/month | `tests/test_eto.py::test_eto_thornthwaite` covers synthetic regression fixtures; `test_eto_thornthwaite_literature_watson` compares against a Thornthwaite (1948) worked example (Watson & Burnett, 1995) in `tests/fixture/pet_literature/`. See `docs/algorithm_refs/pet.md`. | The source supplies daylight hours, not latitude; coordinate reconstruction is not independently validated. |
| PET Hargreaves | Validated | Synthetic-input sanity checks; literature worked example `atol=0.05` mm/day | `tests/test_eto.py::test_eto_hargreaves_with_fixtures` checks synthetic daily-input output shape and bounds; `test_eto_hargreaves_literature_mehta` compares against a Hargreaves-Samani (1985) worked example (Mehta, 2006) in `tests/fixture/pet_literature/`. See `docs/algorithm_refs/pet.md`. | The source supplies extraterrestrial radiation, not latitude/day of year; coordinate reconstruction is not independently validated. |
| PET Penman-Monteith (FAO-56) | Validated | `atol=0.05` mm/day for the final ETo; `0.05`-`0.1` MJ m-2 day-1 for the printed radiation intermediates | `tests/test_pm_eto.py` checks the Chapter 3 radiation, clear-sky, net-radiation, and wind helpers against FAO-56 worked Examples 10-16, and the full derived-input chain (`penman_monteith_eto`) against the daily grass-reference Example 18 (Uccle, ETo = 3.88 mm/day). `tests/test_pm_eto_xarray.py` checks the xarray path against the NumPy path. See `docs/algorithm_refs/pet.md`. | One daily worked example supplies the end-to-end oracle; the monthly Example 17 is not encoded, and the radiation-helper checks reproduce FAO-56's rounded printed intermediates rather than an independent computation. |
| PNP | Validated | Existing fixture tolerances in `tests/test_indices.py` and xarray PNP tests | Percent-of-normal fixture and xarray wrapper tests cover output shape and metadata. | None blocking 3.0.0. |
| PCI | Validated | Existing fixture tolerances in `tests/test_indices.py` and xarray PCI tests | Daily rainfall fixture and scalar xarray wrapper tests cover PCI. | None blocking 3.0.0. |
| EDDI | Validated | NOAA reference tests use `rtol=1e-5`, `atol=1e-5` | Committed paired NOAA PSL monthly reference ET/EDDI fixtures for 1-, 3-, and 6-month Timescales, 1979–2023, pass in `tests/test_noaa_eddi_reference.py`; maximum observed error is `2.43e-6`. | This is one fixed latitude/longitude subset, not a CONUS-wide assessment. |
| KBDI | Reference-reproduced, not independently validated | Published Figure 1 series reproduced with `atol=4.0` hundredths of an inch (table-vs-equation discretization, measured maximum 2.91); GHCN station regression `rtol=1e-12`, `atol=1e-9` mm | `tests/test_fire_kbdi_reference.py` reproduces the Keetch & Byram (1968) Figure 1 worked example from a committed fixture, audits the published net-rain wet-spell rule, and runs a complete 30-year GHCN-Daily PRCP/TMAX record for Fresno, California against frozen values from an independent plain-Python calculator; the xclim (Finkele variant) and WFAS cross-checks skip when their optional inputs are absent. Recurrence, gap, and unit contracts are covered by `tests/test_fire_kbdi.py`. | Figure 1 is an exact oracle only for the 1968 integer table workflow, not the corrected continuous Equation 18 the library implements; xclim is a later Australian/FFDI variant, and no reproducible WFAS point oracle was obtainable. |
| CFFWIS (FFMC, DMC, DC, ISI, BUI, FWI, DSR) | Validated against the NRCan reference implementation | All seven components `atol=1e-9` against full-precision `cffdrs_py` reference outputs; xclim sanity check at `atol=1e-9` for DC/DMC/BUI and `atol=0.5` for the FFMC lineage | `tests/test_fire_cffwis_reference.py` runs the 48-day Van Wagner & Pickett (1985) calibration sample committed as `tests/fixture/cffwis_vwp1985/` against reference values generated at `cffdrs_py` commit `0f57fcc`; the xclim cross-check skips when xclim is absent. Recurrence, gap, unit, and seasonal-carry contracts remain covered by `tests/test_fire_cffwis_moisture.py`, `tests/test_fire_cffwis_behavior.py`, and `tests/test_fire_dc_overwintering.py`. | The committed record is one calibration sample at latitude 40 N, so the external check exercises the northern DMC/DC day-length band only; the other latitude bands are covered by frozen regression vectors, and no independent (non-NRCan-lineage) CFFWIS oracle exists. |
| HDW | Event-discrimination reproduced; published magnitudes characterization only | Daily-series correlation floor `r >= 0.80` (measured 0.8679; floor stored in `tests/fixture/hdw_srock_cedar/provenance.json`); digitized reference uncertainty +/-5 hPa m s-1 | `tests/test_fire_hdw_reference.py` computes `fire.hot_dry_windy()` on the CFSR reanalysis profiles in `tests/fixture/hdw_srock_cedar/`, extracted at the grid point of Srock et al. (2018) Figure 4a (33.0 N, 116.5 W; 12 October through 9 November 2003), and confirms the daily maximum falls on 26 October 2003 -- the paper's significant fire-behavior day for the Cedar Fire -- 17% above the next-highest day. The digitized Figure 4a series is committed alongside the profiles, with its marker pixel coordinates for audit. | The paper publishes only figures, so the reference series is digitized (+/-5 hPa m s-1). Its adiabatically adjusted, independently maximized formulation (Section 3) is not the library's per-level VPD x wind product, and the digitized magnitudes are not reproducible from public metadata: the library's daily maximum exceeds them on five of 29 days, supplied by the 0000 UTC analysis on all five and by the 1200 UTC series alone on three, while the library's 1800 UTC values alone correlate at 0.945 with no exceedances. Magnitudes are therefore characterization, not validation; unit and example-value contracts remain covered by `tests/test_fire.py`. |
| PDSI, PHDI, PMDI, Z-Index (standard Palmer) | Qualified independent external-product validation | nClimDiv comparison ceilings documented in `tests/test_nclimdiv_reference.py`; library fixtures `atol=5e-5`, `rtol=0` | `tests/test_nclimdiv_reference.py::test_pdsi_vs_noaa_nclimdiv_qualified_external_validation` compares `pdsi()` against the operational NOAA NCEI nClimDiv arrays for all 344 climate divisions, January 1895 through December 2022 (measured median absolute difference 0.0127; 86.2% of months within 0.05). `tests/test_palmer.py` retains the library-generated standard-Palmer fixtures as regression coverage. Both are marked `validation`. | Qualified rather than absolute: NCEI publishes two decimal places, and the committed `precips.npy`/`pet.npy` inputs are not byte-identical to NCEI's operational inputs, so the practical agreement floor is ~0.013 rather than 0.005. See "Palmer Validation Classification" below. |
| scPDSI | Independent implementation cross-validation; nClimDiv comparison is characterization only | Wells-lineage oracle fixtures `atol=5e-5`, `rtol=0` in `tests/test_scpdsi.py`; nClimDiv characterization ceilings in `tests/test_nclimdiv_reference.py` | `tests/test_scpdsi.py` compares the four self-calibrating outputs and fitted duration factors against Wells-lineage reference fixtures for all 344 climate divisions. `tests/test_nclimdiv_reference.py::test_scpdsi_vs_noaa_nclimdiv_characterization` characterizes aggregate agreement against nClimDiv, which applies standard Palmer's fixed national calibration rather than per-division self-calibration. Both are marked `validation`. | No qualified external-product reference yet; assessing DRI/WRCC scPDSI for that role remains open (issue #780). The NCEI Fortran and Wells C++ recursion lineages disagree on roughly 6.3% of months by up to 8 PDSI units, and on those months the NCEI-lineage values sit closer to NOAA operational data 88% of the time (`tests/fixture/palmer/provenance.json`). |
| PE (effective precipitation) | Regression only; specification-level (no external oracle adopted) | Eq. (2) identity `atol=0`; endpoint weights at `pytest.approx` defaults; naive double-sum equivalence `rtol=1e-14`, `atol=1e-13` | `tests/test_flood_reference.py` (marked `validation`) checks the Byun & Wilhite (1999) Eq. (2) two-day identity and the harmonic endpoint weights `w₁ = H_D`, `w_D = 1 / D`, and freezes the public signature; `tests/test_flood_pe.py` checks the harmonic double sum against an independent naive sum, full-window NaN behavior, missing-day recovery, and Spatial Block layouts. | No numeric oracle is committed. The tabulated Byun & Wilhite (1999) Table 5 follows the paper's variable summation duration, not the fixed 365-day contract of ADR-0014, and its Hickman, Nebraska input record is not yet acquired (issue #1135); see "Flood Evidence Classification" and `tests/fixture/flood/README.md`. |
| EDI | Regression only; specification-level (no external oracle adopted) | Fixed-window algebra checked against closed-form expectations at NumPy defaults; contract tolerances in `tests/test_flood_edi.py` | `tests/test_flood_edi.py` checks Eq. (9)'s fixed-window algebra (the harmonic factor cancels, leaving `DEP / SD(PE)` per calendar day), calibration-period validation, per-calendar-day baselines, and missing-data behavior; `tests/test_flood_reference.py` freezes the public signature. | The published Table 5 end-to-end values need the variable-duration convention this implementation does not use and the Hickman, Nebraska record (issue #1135); no other numeric oracle was located (`docs/research/flood-oracle-survey.md`). |
| I_F (Flood Index) | Regression only; specification-level (no external oracle adopted) | Declared normalization and year boundaries checked against closed-form expectations in `tests/test_flood_index.py` | `tests/test_flood_index.py` checks `I_F = (PE − mean(PE_max)) / SD(PE_max)` against complete annual maxima, `year_start_month` boundaries, partial-period exclusion, and missing maxima; `tests/test_flood_reference.py` freezes the public signature. | The Deo et al. (2015) full text is paywalled and its abstract describes an exponential kernel the library does not implement — ADR-0014 selects the harmonic double sum — so no reproducible I_F values are available. |
| API (Antecedent Precipitation Index) | Regression only; specification-level (no external oracle adopted) | Exact fixture `api([4, 0, 2], 0.5) == [4, 2, 3]`; constant-input limit `abs=1e-15` | `tests/test_flood_reference.py` (marked `validation`) checks the analytic `P / (1 − k)` limit, and `tests/test_flood_antecedent.py::test_source_recursion_and_closed_form` checks the Kohler & Linsley (1951) Eq. (3) recursion via the exact fixture; the same file covers the state round-trip, gap policies, and Spatial Block contracts outside the `validation` marker. | No accessible source tabulates index values; the retrieved Kohler & Linsley pages are definitional only and the located MIT `ahrapi` implementation was not adopted as a fixture (`docs/research/flood-oracle-survey.md`). |

## Flood Evidence Classification

The flood family — effective precipitation (PE), the Effective Drought Index
(EDI), the Flood Index (I_F), and the Antecedent Precipitation Index (API) — is
**not externally validated**. No reproducible numeric oracle from a published
source is committed, so every flood row above is classified **regression
only**: the tests encode primary-source algebra and public contracts, not an
external reference dataset. Under the vocabulary that
`docs/research/flood-oracle-survey.md` uses and maps to this document, no flood
index reaches *exact reference*, *independent implementation*, or *digitized
figure* status:

- **Exact reference**: Byun & Wilhite (1999) Table 5 is the only located
  tabulated end-to-end reference. It is not usable as this library's oracle: it
  follows the paper's variable summation duration while
  [ADR-0014](docs/adr/0014-flood-family-scientific-conventions.md) fixes the
  365-day harmonic window, and the Hickman, Nebraska input record it requires
  is not yet acquired (issue #1135).
- **Independent implementation**: `tidyindex` (a harmonic-weight rolling
  z-score) and `rrk4910/EDI` (a weighted double-sum variant that omits the
  variable-duration standardization), plus the MIT `ahrapi` package (fixed-`k`
  API), implement related but non-identical algorithms; none reproduces this
  library's conventions and none is committed as a fixture.
- **Digitized figure**: the Byun & Wilhite Figure 2 index panels are
  supplementary lower-tier evidence the survey recommends only behind the
  tabulated values; no figure has been digitized.
- **Regression only**: what is committed. `tests/test_flood_reference.py`
  (marked `validation`) pins the Byun & Wilhite (1999) Eq. (2) two-day
  identity, the harmonic endpoint weights `w₁ = H_D` and `w_D = 1 / D`, and the
  `P / (1 − k)` constant-input limit, and freezes the public signatures. The
  per-index contract suites, run outside the `validation` marker, cover the
  Eq. (3) API recursion, the fixed-window EDI algebra, annual-maximum
  normalization, missing data, and Spatial Block layouts, and
  `tests/test_flood_xarray.py` checks NumPy/eager/Dask equivalence for the beta
  xarray route.

Even the tabulated Table 5, if it is eventually committed, would show
reproduction of the originating group's own arithmetic rather than independent
correctness of the method — the same distinction the scPDSI entry draws between
independent-implementation cross-validation and external-product validation.
`tests/fixture/flood/README.md` lists the retrieval work each index needs before
a numeric oracle can be committed. Flood indices describe flood potential,
not flooding: they are wetness measures upstream of terrain, soils, and
routing, so even a validated value would not be a flood prediction.

## EDDI Fixture Policy

The paired NOAA PSL EDDI fixtures are committed and run in CI. They contain
monthly reference ET and 1-, 3-, and 6-month EDDI for latitude 39.75–39.875
and longitude -105.0–-104.875, with the 1979–2023 table baseline.

Required fixture layout:

```text
tests/fixture/noaa-eddi-1month/
tests/fixture/noaa-eddi-3month/
tests/fixture/noaa-eddi-6month/
```

Each directory must contain:

- `provenance.json`
- `metadata.json`
- `pet_input.npy`
- `eddi_reference.npy`

The provenance file follows `tests/fixture/provenance_schema.json`, cites the
NOAA PSL EDDI time-series source, records checksums, and states the calibration
period. The validation tests compare non-NaN reference values and separately
verify NaN placement. `scripts/prepare_noaa_eddi_fixtures.py` refuses an
unexpected table baseline so fixture refreshes are intentional.

## SPEI Plausibility Classification

`indices.spei()` (gamma distribution, Thornthwaite PET) is compared against
grid-cell-mean CSIC SPEIbase v2.11 series in `tests/test_speibase_reference.py`,
for three CONUS climate divisions spanning an aridity gradient (humid Alabama
0101, subhumid Oklahoma 3405, arid southwest Arizona 0205) and timescales 1,
3, 6 and 12 months. The stored reference is the mean of the per-cell SPEI
values, not an index computed from averaged inputs.

This comparison is a **plausibility check, not external numerical
validation**, and is deliberately not enforced with `atol`/`rtol`. Three
confounds separate the two products:

- SPEIbase v2.11 uses FAO-56 Penman-Monteith PET from CRU TS 4.09, while the
  compared series uses Thornthwaite PET computed from the committed nClimDiv
  temperatures; the two PET families diverge with climate aridity (van der
  Schrier et al. 2011).
- SPEIbase standardizes with the log-logistic distribution, while this
  gamma/Thornthwaite comparison still fits gamma. The separate input-matched
  log-logistic comparison below removes this confound.
- The reference is a 0.5-degree grid average inside a climate division polygon
  while the compared series is the division's station-derived areal average,
  from different precipitation inputs (CRU TS vs. nClimDiv).

The tests therefore assert per-division, per-timescale floors on correlation,
sign agreement, and drought-category agreement, derived from the measurements
recorded in `tests/fixture/speibase/provenance.json`;
`scripts/prepare_speibase_fixtures.py` re-measures those values from the arrays
it just built, refusing to publish a refresh that drifts beyond `0.001` of the
recorded expectations, and regenerates the floors from the new measurements;
`test_floors_keep_documented_slack` pins every floor to within a documented
band below its measurement so a floor cannot be widened to hide a regression.
Measured correlation is weakest in the arid Arizona division (0.82-0.86 across
timescales) -- the direction the PET-family mismatch predicts -- and strongest
in subhumid Oklahoma (0.95-0.97), with humid Alabama (0.92-0.94) between them.
With only three divisions the ordering cannot isolate the PET-family confound
from the distribution and precipitation-support confounds, which is precisely
why no tight numerical gate is defensible: it would either fail on that
legitimate climate-dependent bias or have to be widened until it no longer
tested anything.

The compared climate_indices series is calibrated against the full 1901-2022
period of record, matching the committed inputs, while SPEIbase v2.11 is
standardized on its own 1901-2024 record; the fixtures retain the 1901-2022
overlap. The dataset survey behind this comparison lives at
`docs/research/spei-dataset-survey.md` on branch `research/spei-dataset-survey`.

## SPEI Log-Logistic Like-for-Like Comparison

`tests/test_speibase_like_for_like.py` removes the three confounds above by
feeding `indices.spei()` the CRU TS 4.09 monthly precipitation and FAO-56
Penman-Monteith PET that SPEIbase v2.11 was computed from, committed as
per-cell fixtures in `tests/fixture/speibase_cru_ts/`. It standardizes each
0.5-degree cell with the log-logistic (generalized logistic) distribution over
1901-2024 and averages the per-cell series inside the same three climate
divisions, matching SPEIbase's `R/functions.R` (`spei.nc`) per-cell-then-mean
pipeline.

This is a cross-implementation check of a port, not of an algorithmically
independent estimator: `climate_indices.lmoments.fit_glo` is a Python
translation of the same `lmom` PELGLO routine R `SPEI` calls, using the same
unbiased-PWM L-moments, so close agreement is expected.

The agreement is near-exact -- correlation 0.9996-0.99999 and mean absolute
difference 0.003-0.016 across divisions and timescales -- and the tests assert
tight per-series floors (correlation >= 0.97) instead of the loose plausibility
floors. The residual is small and dominated by a PET day-length convention
difference, not by the fit: this test multiplies CRU TS PET (mm/day) by the
leap-aware calendar month lengths, matching SPEIbase's public `R/functions.R`
(`etp * ndays` with `Hmisc::monthDays`), while a near-uniform month length
reproduces the committed SPEIbase v2.11 grids two to four times more closely,
so v2.11's effective generation convention could not be confirmed from public
code. Floors are derived from
`tests/fixture/speibase_cru_ts/provenance.json`, and
`scripts/prepare_speibase_cru_ts_inputs.py` re-measures the agreement from the
arrays it builds, refusing to publish a refresh that drifts beyond `0.001` of
the recorded expectations; `test_floors_keep_documented_slack` pins every floor
to a documented band below its measurement.

The earlier gamma/Thornthwaite comparison is retained because it exercises a
different distribution and PET path against the same reference; it remains a
plausibility check, not numerical validation.

## SPEI Log-Logistic Cross-Implementation

`indices.spei()` with `Distribution.loglogistic` is compared against the R
`SPEI` package (`distribution = "log-Logistic"`, `fit = "ub-pwm"`) in
`tests/test_spei_loglogistic_reference.py`, over the same three CONUS climate
divisions and timescales 1/3/6/12 as the SPEIbase comparison and over a short,
highly skewed synthetic series with one constant calendar month. The committed
R output lives under `tests/fixture/spei_loglogistic/`; `provenance.json`
records the pinned R package versions, the fixture checksum, and the
tolerances, and `scripts/prepare_spei_loglogistic_fixtures.py` refuses a refresh
from a different `SPEI`/`TLMoments`/`lmom` version.

This is an **independent-implementation comparison**, not a plausibility check:
`SPEI` is the reference implementation of the log-logistic SPEI and the
distribution the CSIC SPEIbase is built from (Vicente-Serrano et al., 2010;
Beguería et al., 2014). Both sides receive the identical `(P − PET) + 1000`
water balance -- the offset `indices.spei` applies internally -- so the
comparison isolates the distribution fit from the scaling and offsetting steps,
and `SPEI` fits each calendar step with the same unbiased-PWM L-moments
estimator (`TLMoments::PWM` + `lmom::pelglo`) that `lmoments.fit_glo`
implements. The tolerances (`series_atol = 1e-6`; `parameter_rtol = 1e-8`,
`parameter_atol = 1e-6`) therefore absorb only cross-platform numerical wobble,
not algorithmic differences: measured agreement is at ~1e-12 for the
standardized series and ≤~1e-10 for the fitted parameters. Together the
fixtures span the edge cases the ticket names -- a short sample and extreme
L-skewness (GLO shape to about -0.8) from the synthetic series, a negative
underlying location parameter from the real arid divisions, and a constant
calendar step that is unfittable at scale 1 (R reports NA, this library a
zeroed scale and NaN values; at scale 6 the constant month enters non-constant
multi-month sums and is fittable). `SPEI` does not clip its output to
`[-3.09, 3.09]`, so the comparison covers finite values inside that range and
asserts this library's clip on the out-of-range tail separately (the synthetic
scale-1 fit reports one value of `-inf`).

## Palmer Validation Classification

The Palmer evidence splits by product family, because the two consume different
reference artefacts.

**Standard Palmer (`pdsi()`: PDSI, PHDI, PMDI, Z-Index)** is qualified
independent external-product validation against the operational NOAA NCEI
nClimDiv arrays in `tests/fixture/nclimdiv/`, enforced by
`tests/test_nclimdiv_reference.py::test_pdsi_vs_noaa_nclimdiv_qualified_external_validation`.
Measured over 528,384 division-months, `pdsi()` agrees with nClimDiv to a median
absolute difference of 0.0127, with 86.2% of months within 0.05. The
qualification is deliberate and recorded rather than averaged away: NCEI
publishes two decimal places, and the committed `precips.npy`/`pet.npy` inputs
are not byte-identical to NCEI's operational inputs, so the practical agreement
floor is near 0.013 rather than the 0.005 that published precision alone would
permit. The library-generated standard-Palmer fixtures (`tests/fixture/palmer/`)
remain in the suite as regression coverage only -- their provenance identifies
the source as "climate_indices library reference output" -- and they are no
longer the family's only evidence.

**scPDSI (`scpdsi()`)** has independent implementation cross-validation, not
external-product validation: `tests/test_scpdsi.py` compares the four
self-calibrating outputs and fitted duration factors against the Wells-lineage
C++ oracle fixtures at `atol=5e-5`, `rtol=0`. The nClimDiv comparison in
`tests/test_nclimdiv_reference.py::test_scpdsi_vs_noaa_nclimdiv_characterization`
stays explicitly classified as characterization, because nClimDiv uses standard
Palmer's fixed national K-factors while `scpdsi()` self-calibrates per division.
The two recursion lineages disagree on roughly 6.3% of months by up to 8 PDSI
units (the NCEI-lineage values are closer to NOAA operational data on 88% of
those months; see `tests/fixture/palmer/provenance.json`), so neither series can
serve as the other's tolerance oracle.

No further standard-Palmer external reference is required for this validation
contract. scPDSI external-product validation remains open pending the DRI/WRCC
assessment in issue #780.

## scPDSI Calibration-Anchor Measurement

Wells' self-calibration methodology tunes each climate division's duration
and K-prime factors so that, by construction, the calibration-period 2nd and
98th percentiles of the resulting scPDSI series land at exactly -4.0000 and
+4.0000. (`scpdsi()` takes those percentiles of the PDSI output and rescales
the Z-index to match; the Z-index is the quantity adjusted, not the quantity
anchored.) Measuring how closely `scpdsi()` reproduces that anchor across all
344 climate divisions (calibration period 1931-1990) verifies the
self-calibration percentile-scaling step independently of the NOAA nClimDiv
comparison in `tests/test_nclimdiv_reference.py` (see the "Per-Index Evidence"
table above), which instead measures agreement with an external product that
uses different (fixed, national) K-factors.

| Percentile anchor | Median | Mean | Median deviation | Max deviation | Within +/-0.05 | Within +/-0.25 |
| --- | --- | --- | --- | --- | --- | --- |
| 2nd (target -4.0000) | -4.0000 | -4.0142 | 0.0102 | 0.9730 | 74% | 93% |
| 98th (target +4.0000) | +4.0000 | +3.9829 | 0.0130 | 1.5421 | 70% | 91% |

The medians land on the target exactly while the maxima do not, because
`scpdsi()` applies a fixed three rescaling passes rather than iterating to a
fixed point. Divisions whose recursion has not settled within three passes
retain the residual deviations shown above; this is the reference behaviour,
not a defect. These figures are regenerated and bounded by
`test_scpdsi_calibration_anchor_lands_on_target` in
`tests/test_nclimdiv_reference.py`, so they fail a test rather than rotting
silently if the rescaling changes.

This is regression coverage over the committed fixture set, not independent
scientific validation: it confirms the self-calibration percentile-scaling
arithmetic behaves as Wells describes for these 344 divisions, not that the
resulting scPDSI values themselves match an external authoritative source.
See "Palmer Validation Classification" above for the latter.

## References

- Beguería, S., Vicente-Serrano, S. M., Reig, F., and Latorre, B. (2014).
  Standardized precipitation evapotranspiration index (SPEI) revisited:
  parameter fitting, evapotranspiration models, tools, datasets and drought
  monitoring. International Journal of Climatology, 34(10), 3001-3023.
  https://doi.org/10.1002/joc.3887
- Byun, H.-R., and Lee, D.-K. (2002). Defining three rainy seasons and the
  hydrological summer monsoon in Korea using Available Water Resources Index.
  Journal of the Meteorological Society of Japan, 80(1), 33-44.
  https://doi.org/10.2151/jmsj.80.33
- Byun, H.-R., and Wilhite, D. A. (1999). Objective quantification of drought
  severity and duration. Journal of Climate, 12(9), 2747-2756.
  https://doi.org/10.1175/1520-0442(1999)012%3C2747:OQODSA%3E2.0.CO;2
- Deo, R. C., Byun, H.-R., Adamowski, J. F., and Kim, D.-W. (2015). A real-time
  flood monitoring index based on daily effective precipitation and its
  application to Brisbane and Lockyer Valley flood events. Water Resources
  Management, 29(11), 4075-4093. https://doi.org/10.1007/s11269-015-1046-3
- Kohler, M. A., and Linsley, R. K. (1951). Predicting the runoff from storm
  rainfall. U.S. Weather Bureau Research Paper No. 34.
- Hobbins, M. T., Wood, A., McEvoy, D. J., Huntington, J. L., Morton, C.,
  Anderson, M., and Hain, C. (2016). The Evaporative Demand Drought Index.
  Part I. Journal of Hydrometeorology, 17, 1745-1761.
  https://doi.org/10.1175/JHM-D-15-0121.1
- Palmer, W. C. (1965). Meteorological Drought. U.S. Weather Bureau Research
  Paper No. 45. https://www.droughtmanagement.info/literature/USWB_Meteorological_Drought_1965.pdf
- NOAA Physical Sciences Laboratory EDDI archive:
  https://downloads.psl.noaa.gov/Projects/EDDI/
- van der Schrier, G., Jones, P. D., and Briffa, K. R. (2011). The sensitivity
  of the PDSI to the Thornthwaite and Penman-Monteith parameterizations for
  potential evapotranspiration. Journal of Geophysical Research: Atmospheres,
  116, D03106. https://doi.org/10.1029/2010JD015001
- Vicente-Serrano, S. M., Beguería, S., and López-Moreno, J. I. (2010). A
  multiscalar drought index sensitive to global warming: the Standardized
  Precipitation Evapotranspiration Index. Journal of Climate, 23(7),
  1696-1718. https://doi.org/10.1175/2009JCLI2909.1
- Van Wagner, C. E., and Pickett, T. L. (1985). Equations and FORTRAN
  program for the Canadian Forest Fire Weather Index System. Canadian
  Forestry Service, Forest Technical Report 33.
  https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/19973.pdf

## docs/xarray_compatibility.md

# xarray Compatibility Matrix

:::{warning}
The xarray API is beta in 3.0.0 and is promoted no earlier than 3.1.0
([ADR-0012](docs/adr/0012-xarray-api-stays-beta-through-3.0.0.md)). Numerical results
are expected to match the stable NumPy API, while parameter inference, metadata,
and coordinate behavior may change with a minor version, never in a patch release.
:::

| Feature | Supported | Coverage |
| --- | --- | --- |
| `DataArray` inputs for SPI | Yes | `tests/test_xarray_equivalence.py` compares 1-D and gridded xarray outputs to NumPy. |
| `DataArray` inputs for SPEI | Yes | `tests/test_xarray_equivalence.py` compares xarray precipitation/PET inputs to NumPy. |
| `DataArray` inputs for EDDI | Yes | EDDI wrapper and metadata tests cover the typed public API. |
| `DataArray` inputs for flood PE, EDI, I_F, and API | Yes | `flood.effective_precipitation()`, `flood.edi()` (also `climate_indices.edi()`), `flood.flood_index()`, and `flood.antecedent_precipitation_index()` accept daily precipitation/PE with CF units converted to mm. For PE, EDI, and I_F, Gregorian dates are converted to all-leap days before each NumPy computation and restored afterward; API consumes the daily series as given. PE must be computed before EDI/I_F and use the same duration. EDI and I_F infer calendar and Calibration Period years when omitted; I_F still requires `year_start_month`. API keeps its recursive-state contract (`initial_state`/`return_state`, `spin_up`, `nan_policy`/`max_gap_days`), and the returned state stays NumPy. A complete `time` Dask chunk is mandatory; spatial chunks run independently. `tests/test_flood_xarray.py` compares NumPy, eager xarray, and Dask outputs. These indices describe flood potential, not observed flooding. When chaining xarray PE into EDI/I_F, the second adapter reconstructs non-leap February 29 by interpolating **PE**, not the rainfall that produced it; the chained result can differ from applying both NumPy kernels to a single all-leap rainfall array ([#1147](https://github.com/monocongo/climate_indices/issues/1147)). |
| `DataArray` inputs for PET Thornthwaite | Yes | xarray adapter tests cover scalar and spatial latitude handling. |
| `DataArray` inputs for PET Hargreaves | Yes | xarray adapter tests cover aligned daily temperature inputs. |
| `DataArray` inputs for PNP | Yes | PNP wrapper tests cover scale handling and metadata. |
| `DataArray` inputs for PCI | Yes | PCI uses a manual scalar-output wrapper. |
| `DataArray` inputs for KBDI | Yes | `fire.kbdi()` resolves the `kbdi`/`kbdi_imperial` CF entry per call from `units`; supports `return_state`/`initial_state`, `spin_up`, `nan_policy`/`max_gap_days`, and CF `units`-attribute unit inference for the weather inputs and an attributed mean annual climatology. An attached time coordinate must be consecutive daily observations per input; a dimension-only time axis is aligned positionally. Overlapping timesteps are inner-aligned, dropped timesteps emit `InputAlignmentWarning`, and alignment that would drop non-time coordinates is rejected. See `tests/test_fire_kbdi.py`. |
| `DataArray` inputs for HDW | Yes | `fire.hot_dry_windy()` reduces the caller-named `level_dim` (default `"level"`) to the layer maximum via `xr.apply_ufunc`; CF metadata comes from the `hdw` registry entry, and a CF `units` attribute on temperature is converted to Celsius. Not recursive: every other dimension, including `time` if present, is a plain passthrough with no cadence requirement. Dimensions shared by the inputs must carry identical, identically ordered coordinate labels, because the inputs are matched with xarray's exact join rather than aligned or reindexed. `level_dim` must be a single Dask chunk; other dimensions chunk freely. See `tests/test_fire.py`. |
| `DataArray` inputs for CFFWIS | Yes | `fire.cffwis()` returns one `Dataset` variable per selected output (`ffmc`, `dmc`, `dc`, `isi`, `bui`, `fwi`, `dsr`), each carrying its own registry CF metadata. Latitude is inferred from a `lat`/`latitude` coordinate when not supplied and month from the datetime `time` coordinate; CF `units` attributes on temperature and precipitation are converted. Per-cell day-length bands come from broadcasting latitude through `xr.apply_ufunc`, one call per Dask spatial block. A daily coordinate clearly not at noon emits `ClimateIndicesWarning`, while sub-daily coordinates and a multi-chunk `time` dimension are rejected. With `return_state=True` the result is a `CFFWISResult` whose selected fields are DataArrays (`None` for components omitted from `outputs`) and whose state stays NumPy. See `tests/test_fire_cffwis_behavior.py`. |
| `DataArray` inputs for Haines | Yes | `fire.haines_index()` scores caller-selected levels elementwise through `xr.apply_ufunc` with no core dimension, so no dimension has a Dask chunk constraint. CF metadata comes from the per-variant `haines_low`/`haines_mid`/`haines_high` registry entry, including `climate_indices_variant`, and a CF `units` attribute on any temperature input is converted to Celsius. Not recursive: every dimension, including `time`, is a plain passthrough with no cadence requirement. Dimensions shared by the inputs must carry identical, identically ordered coordinate labels, because the inputs are matched with xarray's exact join rather than aligned or reindexed. The elevation-driven variant selection of `fire.haines_index_from_profile()` is NumPy-only. See `tests/test_fire_haines.py`. |
| `DataArray` inputs for PDSI (Palmer) | Yes | `climate_indices.pdsi()` returns an `xr.Dataset` with one variable per index (`pdsi`, `phdi`, `pmdi`, `z_index`), each carrying its own registry CF metadata. AWC is a broadcast input like PET's latitude; a 3-D input reaches `palmer.pdsi()` as one `(time, *cells)` block per Dask spatial block. scPDSI has no xarray entry point ([ADR-0011](docs/adr/0011-palmer-spatial-block-and-per-location-scpdsi.md)). See `tests/test_spatial_kernel.py`; the manual-rewrap pattern remains in `notebooks/palmer_indices_xarray.ipynb`. |
| `DataArray` inputs for fit diagnostics | Yes | `climate_indices.fit_diagnostics()` returns an `xr.Dataset` of the fitted parameters, `prob_zero`, `n_valid`, `ks_statistic`, and `ks_p_value` over a `month`/`dayofyear` dimension plus the input's cell dimensions, with `distribution_used` per cell only (no calendar-step dimension). A Pearson Type III request keeps both parameter families, NaN where a block-level fall back makes one inapplicable; the fall back is decided per fitted block, so a Dask chunk layout can change it, as in the index adapters. See `tests/test_fit_diagnostics_xarray.py`. |
| Coordinate preservation | Yes | Adapter tests verify time and spatial coordinates are preserved. |
| CF-style metadata | Yes | `CF_METADATA` registry and adapter tests verify `long_name`, `units`, `references`, version, and history attributes. |
| Dask-backed arrays | Yes, constrained | The time dimension must be a single chunk, except for the PET adapters, which rechunk a split `time` internally. Adapter tests verify detection and validation. Chunk axes the CLI leaves to Dask (`"auto"`, the KBDI and flood-index grid and division inputs) resolve against a 100 MB array chunk budget rather than Dask's own default; axes pinned with `-1` are unaffected. See [Chunking guidance for gridded indices](#chunking-guidance-for-gridded-indices) for block sizes. |
| Spatial (gridded) kernels | Yes, for SPI, SPEI, PET, EDDI, percentage of normal, PDSI, and flood PE/EDI/I_F/API | `spi`/`spei` receive a time-major `(time, *cells)` block and fit every cell in one pass, so a gridded gamma run costs one kernel call per Dask block instead of one call per cell ([ADR-0009](docs/adr/0009-spatial-block-declaration.md)). The PET entry points (`pet_thornthwaite`, `pet_hargreaves`) hand `indices.pet` and `eto.eto_hargreaves` the same block layout with the per-cell latitude array, so a gridded PET run costs one call per block ([#941](https://github.com/monocongo/climate_indices/issues/941)). `indices.eddi` and `indices.percentage_of_normal` rank, and divide by the normals of, every cell of the block in one pass ([#942](https://github.com/monocongo/climate_indices/issues/942)). `palmer.pdsi` receives the block with AWC broadcast per cell, so a gridded PDSI run costs one recursion per block ([#1016](https://github.com/monocongo/climate_indices/issues/1016)). The NumPy API keeps 2-D input on the legacy `(years, periods)` reading; a 3-D or higher input is a block for `spi`/`spei`/`pet`/`pdsi` unless its first cell axis is the one shape that is ambiguous with that reading (12 or 366, i.e. a `(years, periods, *cells)` array), which must be declared, and for EDDI and percentage of normal every 3-D input must be declared with `spatial_time_major=True`. Inputs with a single non-core dimension keep the per-cell path. See `tests/test_spatial_kernel.py`. |
| Calendar semantics | Yes, constrained | Standard/gregorian/proleptic_gregorian `datetime64` only; monthly input must begin in January and daily input on January 1. Daily values are converted to the 366-day calendar (February 29 synthesized from February 28 and March 1) and restored afterward. A partial final year is supported; `cftime` calendars are rejected. See [ADR-0004](docs/adr/0004-xarray-calendar-semantics.md). |
| Automatic temporal inference | Yes | Monthly and daily time-coordinate inference is covered by adapter tests. |
| Multi-input alignment | Yes | SPEI, Hargreaves, and Penman-Monteith align their inputs with an inner join on the time dimension and emit a warning when timesteps are dropped. Shared cell dimensions must carry matching coordinates: intersecting them would silently drop grid cells, so a mismatch raises `CoordinateValidationError` instead. |

## Stateful fire indices

KBDI's adapter (`fire.kbdi()`) and the CFFWIS adapter (`fire.cffwis()`, #807)
are shipped, following the recursive contract in
[ADR-0006](docs/adr/0006-fire-recursive-state-and-execution.md). The weather-only
Fosberg Fire Weather Index is available today too, via the NumPy layer with
no xarray adapter needed (it is elementwise with no dimension to reduce).
Hot-Dry-Windy's adapter (see the main table above) is not recursive either,
but it does reduce a named vertical dimension, so it does not follow the
stateful contract below; it needs only that dimension in a single Dask chunk.
Haines' adapter is elementwise like Fosberg's, but it still goes through
`xr.apply_ufunc` -- with no core dimension, and therefore no chunk constraint
-- so that it can carry per-variant registry metadata and CF `units`
conversion; its elevation-driven variant selection stays in the NumPy layer.

These adapters are recursive: each daily value needs its predecessor.
Dask-backed inputs must therefore keep the complete `time` dimension in one
chunk for every time-varying weather variable. Spatial chunks remain
supported and parallelize: KBDI's and CFFWIS's NumPy cores already vectorize
over an arbitrary spatial shape, so their adapters dispatch one call per Dask
spatial chunk with the full `time` axis, rather than looping per grid cell.
Multi-chunk time input raises `CoordinateValidationError` with
`data = data.chunk({'time': -1})`; adapters never rechunk implicitly, because
doing so can materialize a large daily history.

Callers use the returned state to append later observations without
recomputing the archive. The state is a NumPy-layer value object rather than
an xarray `Dataset`, so its arrays carry the computational spatial shape but
no coordinates, even when `values` on the same result is a `DataArray`.

Missing observations follow
[ADR-0007](docs/adr/0007-fire-missing-data-policy.md): nothing is interpolated
implicitly. The default `nan_policy="propagate"` poisons the state at the
first interior gap and yields NaN outputs from there on;
`nan_policy="bridge"` with `max_gap_days=N` skips runs of at most `N` missing
days. Adapters forward both keyword arguments to the NumPy core so the policy
applies per cell along the time axis.

## Operational Guidance

- Use NumPy APIs for stable production integrations that cannot absorb beta
  interface changes.
- Use xarray APIs for labeled, gridded workflows where coordinate preservation
  and metadata are more valuable than strict interface stability.
- For a complete runnable gridded SPI and SPEI example with the chunking and
  scheduler guidance applied, see [Gridded Performance](docs/performance.md).
- Keep Dask chunks spatial when possible and leave `time` as one chunk before
  calling index functions. Spatial chunks set the parallelism granularity: SPI and
  SPEI fit a whole `(time, *cells)` block per call, and the gridded path requires
  the full `time` axis in that block. Chunk size is also the memory lever — one block
  is held whole inside the fit, with a few `O(years x periods x cells)` temporaries,
  so a dense continental grid passed in one piece costs several times its own size.
  See [Chunking guidance for gridded indices](#chunking-guidance-for-gridded-indices)
  for measured block sizes.
- Axes the CLI leaves to Dask (`"auto"`, currently the KBDI and flood-index grid and division inputs)
  resolve against a 100 MB array chunk budget
  (`climate_indices._cli.DEFAULT_ARRAY_CHUNK_SIZE`) instead of Dask's own 128 MiB
  default; axes pinned with `-1` are unaffected. Change the budget with explicit
  `chunks` when opening your own dataset, with the `DASK_ARRAY__CHUNK_SIZE`
  environment variable, or with `dask.config.set({"array.chunk-size": ...})` — a
  configured budget is honored rather than overwritten.
- Choose the scheduler at the materialization call: the xarray API never imports
  or configures Dask ([ADR-0002](docs/adr/0002-multiprocessing-cli-dask-xarray.md)). The
  gridded kernels are CPU-bound Python/scipy work, and their Python-level portion
  does not run in parallel under the default threaded scheduler, so materialize a
  lazy result on worker processes — given `precip`, whose `time` dimension is a
  single chunk:

  ```python
  from climate_indices import spi
  from climate_indices.indices import Distribution

  spi_lazy = spi(
      values=precip.chunk({"time": -1, "lat": 20, "lon": 20}),
      scale=3,
      distribution=Distribution.gamma,
  )
  spi_grid = spi_lazy.compute(scheduler="processes")
  ```

  In a script rather than a notebook, that call belongs under
  `if __name__ == "__main__":` while Dask spawns its workers — the `processes`
  scheduler's default start method on every platform — because each worker
  re-imports the entry module. A Zarr write takes the same scheduler through the
  delayed write:
  `spi_lazy.to_zarr(path, compute=False).compute(scheduler="processes")`.
  NetCDF writes do not: their backend lock is built for the scheduler that is
  active when the write graph is built, and the default one cannot be pickled to
  worker processes, so load the result first and write it in memory when the full
  result fits, or let a distributed client stream the write for larger results.
- Every `.compute(scheduler="processes")` call builds and tears down its own
  process pool, so a computation short relative to that start-up spends most of
  its wall clock there: on the 38 x 87 reference grid the SPI pass measured 1.37 s
  with the fresh one-worker pool against 0.76 s with a pre-created, warmed pool,
  and 1.26 s against 0.23 s at eight workers
  ([#928](https://github.com/monocongo/climate_indices/issues/928) harness,
  recorded on [#927](https://github.com/monocongo/climate_indices/issues/927)).
  The serial reference and the fresh-pool rows for SPI, SPEI, PET and EDDI are
  tabulated in
  [the #929 before/after table](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md#beforeafter-on-the-reference-grid-929).
  Keep a `dask.distributed.Client` alive across calls — the [tutorial's client
  cell](https://github.com/monocongo/climate_indices/blob/main/notebooks/zarr_dask_spi_spei.ipynb)
  does, and `distributed` ships with the `dev` extra — or pass a pre-created pool
  when the work is short or repeated.
- The threaded scheduler remains the right choice for I/O-bound or small
  in-memory work — opening a store, a reduction such as `.mean("time")`, or a grid
  small enough that serializing each block to a worker process costs more than
  the parallelism saves.
- The canonical lazy xarray/Dask SPI/SPEI workflow is the teaching notebook
  `notebooks/zarr_dask_spi_spei.ipynb`: the public typed API on Dask-backed DataArrays
  (`xr.apply_ufunc(..., dask="parallelized")`), one full time chunk with spatial
  chunks driving task parallelism, and precipitation/PET exact-aligned at
  preparation time so SPEI never relies on coordinate intersection.
- Persist results to a separate consolidated Zarr v2 store rather than back into
  the prepared inputs. The notebook's `data/e2e/climate_indices_output.zarr` is
  the sample path: `float32` index variables keep the API's `long_name`,
  dimensionless `units`, `references`, `scale`, `distribution`, Calibration
  Period years, `climate_indices_version`, and appended `history`, plus a
  `periodicity` attribute that names the Timescale unit. CF defines no
  `standard_name` for drought indices, so inherited input names such as
  `precipitation_amount` are dropped rather than mislabeling the result.
- Rerun persistence is a local, single-writer replacement: stage beside the
  target, then rename the previous output aside and swap, so a failed
  calculation or write leaves the completed store intact. Close readers before
  rerunning, keep input/output/staging paths non-overlapping, and do not run
  concurrent writers against one output path; the final directory swap is not
  crash-atomic.
- Reopen persisted results with `xr.open_zarr(..., consolidated=True)` in a fresh
  handle: inspect metadata and chunks without loading data, compute only the
  selected diagnostics, then close the handle. The reopened store is independent
  of the prepared inputs, so analysis does not recompute SPI/SPEI.
- Run the notebook CI command before publishing examples:
  `uv run jupyter nbconvert --execute --to notebook --inplace notebooks/xarray_getting_started.ipynb notebooks/palmer_indices_xarray.ipynb notebooks/eddi_xarray.ipynb`.
- Execute the end-to-end tutorial from a fresh kernel before publishing it:
  `bash scripts/smoke_e2e_notebook.sh` prepares the pinned sample inputs and
  runs `notebooks/zarr_dask_spi_spei.ipynb` end to end into a scratch directory,
  so a local run never rewrites the committed notebook. CI executes the
  tutorial only via `tests/test_e2e_with_dask.py`, which runs its cells against
  a synthetic store; that suite deliberately skips the Dask Client cell (#829),
  so this local run is the only full-notebook check.

## Chunking guidance for gridded indices

Three rules fix the shape of a chunked input:

1. **Keep `time` in a single chunk.** Distribution fitting needs each cell's
   complete series, so the SPI/SPEI adapter path rejects a split `time` with
   `CoordinateValidationError`
   ([ADR-0003](docs/adr/0003-dask-time-dimension-single-chunk.md)) and the stateful
   fire adapters require the same for their weather inputs. The PET adapters are
   the exception: they pass `allow_rechunk=True` and rechunk a split `time`
   internally, paying that copy inside every call.
2. **Spatial chunks are the parallelism and memory lever.** A Dask block is
   fitted whole, and the fit materializes the reshaped `(years, periods, *cells)`
   block plus its per-period intermediates, so the working set grows with the
   cells in a block, not with the size of the grid. The PET spatial kernels take
   one block per call as well, so the same block sizing applies to them.
3. **Chunk the shared dimensions of a multi-input index the same way.** SPEI
   aligns precipitation and PET with `xr.align(join="inner")`, but the adapter
   validates shared non-time dimensions first: a mismatched `lat`/`lon`
   coordinate raises `CoordinateValidationError`, and only `time` is left for
   the inner join to intersect; the adapter's `InputAlignmentWarning` reports
   dropped `time` steps. That join does not rechunk, and each input's `time`
   chunking is validated independently, so a chunk-layout mismatch on a shared
   dimension survives to the compute: `xr.apply_ufunc(...,
   dask="parallelized")` makes Dask unify the chunks, inserting a
   `rechunk-merge` stage that copies the data when the graph computes — an extra
   copy inside every index that consumes the mismatched pair. Give PET the
   layout of the dimensions it shares with precipitation (the one
   `notebooks/zarr_dask_spi_spei.ipynb` stores):

   ```python
   pet = pet.chunk({dim: blocks for dim, blocks in precip.chunksizes.items() if dim in pet.dims})
   ```

   PET inputs with fewer dimensions than precipitation broadcast, so a
   `(time,)` series or a `(time, lat)` field needs no chunking on the dimensions
   it lacks.

For the 40-year monthly gamma SPI/SPEI path the fit's working set is about
60 KB per cell — roughly 16x the 3.84 KB the cell's own 480-step float64 series
occupies. On the #893 reference grid (38 x 87 = 3,306 cells, 480 months, Python
3.13.13, macOS arm64, `scheduler="synchronous"`, peak RSS measured above a
dry-run baseline of the same process):

| Spatial chunk | Blocks | Working set | SPI-3 wall time |
| --- | ---: | ---: | ---: |
| `38 x 87` (one block) | 1 | 205 MB | 0.98 s |
| `20 x 20` | 10 | 24 MB | 0.30 s |
| `10 x 10` | 36 | 6 MB | 0.32 s |

Recommendations for monthly grids:

- Keep each block under roughly 100 MB of working set — around 1,600 cells at
  this per-cell cost. That figure is a per-block measurement taken with
  `scheduler="synchronous"`, one block resident at a time; a threaded or
  distributed worker can hold several ready blocks plus their inputs and
  outputs, so budget from the memory a worker can spare per block it runs at
  once. The 38 x 87 reference grid is not large, but as one block it costs
  205 MB and is about 3x slower than the same grid in `10 x 10` blocks.
- Size blocks by cell count rather than shape: the fit's working set grows with
  the cells a block holds, and with a fixed per-block budget the task count
  comes from the grid size, so aspect ratio changes neither. Chunk both spatial
  dimensions (the `10 x 10` to `20 x 20` range here) rather than long rows so
  the block size tracks the budget; chunking a single dimension ties it to the
  grid's row length, which can overshoot the budget on wide grids.
- Leave at least a few blocks per worker so the scheduler has work to balance;
  `benchmarks/parallel_scaling.py` re-chunks spatial dims to one block per
  worker for its strong-scaling runs.
- Daily grids carry up to 366 steps per cell-year instead of 12, so the per-cell
  working set is roughly 30x larger; at the same ~100 MB ceiling that is about a
  `7 x 7` block where a monthly grid uses `10 x 10`.

Rechunk once, at read or prepare time, when the stored layout differs from the
shape the computation wants — a Zarr store with one `time` chunk per year, or a
single chunk spanning the whole grid. `.chunk(...)` only sets the layout of a
lazy graph; the copy happens when the graph is computed, so persist the array
(or write the rechunked layout back to the store) before the index calls to pay
the copy once. Without that, every separate index graph that consumes the array
can reread and redo the rechunk work. The teaching notebook
`notebooks/zarr_dask_spi_spei.ipynb` prepares its store with `time` as one chunk
and `10 x 10` spatial blocks, the layout the table above measures. The
`{"lat": 50, "lon": 50}` shape in ADR-0003 is legal chunking, not a size
recommendation: on a 40-year monthly series it is roughly 150 MB per block,
above the ~100 MB budget recommended here.

## docs/algorithm_refs/eddi.md

# Evaporative Demand Drought Index (EDDI)

## Scope

EDDI is a standardized drought index based on evaporative demand. In this
package, EDDI is computed from PET-like evaporative-demand inputs and follows a
non-parametric ranking workflow:

1. Clip physically invalid negative PET values to zero.
2. Aggregate values over the requested time scale.
3. Rank aggregated values within each calendar period of the calibration
   window.
4. Convert empirical probabilities to standardized normal deviates.
5. Clip output to the package-standard range `[-3.09, 3.09]`.

## Public API

Use `climate_indices.eddi()` for the typed public API. It accepts NumPy arrays
and beta xarray `DataArray` inputs. NumPy callers must provide temporal
parameters explicitly; xarray callers can infer them from the time coordinate.

## Validation

The always-on EDDI tests cover shape handling, missing data behavior, clipping,
empirical ranking properties, and xarray metadata preservation. External NOAA
PSL comparison tests live in `tests/test_noaa_eddi_reference.py` and are marked
`validation`.

The committed NOAA comparison fixtures contain paired monthly reference ET and
EDDI at 1-, 3-, and 6-month Timescales for latitude 39.75–39.875 and longitude
-105.0–-104.875, using the 1979–2023 table baseline:

- `tests/fixture/noaa-eddi-1month/`
- `tests/fixture/noaa-eddi-3month/`
- `tests/fixture/noaa-eddi-6month/`

They are generated from the NOAA PSL EDDI time-series table by
`scripts/prepare_noaa_eddi_fixtures.py` and run in CI.

## Tolerance

NOAA reference comparisons use `rtol=1e-5` and `atol=1e-5`. The committed
fixture comparison's maximum observed error is `2.43e-6`, below both tolerances.

## References

- Hobbins, M. T., Wood, A., McEvoy, D. J., Huntington, J. L., Morton, C.,
  Anderson, M., and Hain, C. (2016). The Evaporative Demand Drought Index.
  Part I: Linking Drought Evolution to Variations in Evaporative Demand.
  Journal of Hydrometeorology, 17, 1745-1761.
  https://doi.org/10.1175/JHM-D-15-0121.1
- NOAA Physical Sciences Laboratory EDDI resources:
  https://psl.noaa.gov/eddi/
- NOAA PSL EDDI archive:
  https://downloads.psl.noaa.gov/Projects/EDDI/

## docs/algorithm_refs/pet.md

# Potential Evapotranspiration (PET)

## Scope

`climate_indices` implements three PET methods:

- **Thornthwaite (1948)** — monthly PET from mean monthly air temperature,
  latitude, and calendar-derived daylight hours.
- **Hargreaves-Samani (1985)** — daily PET (ETo) from minimum, maximum, and
  mean daily air temperature, latitude, and day of year (via
  extraterrestrial radiation), following Equation 52 in Allen et al. (1998)
  FAO Irrigation and Drainage Paper 56.
- **FAO-56 Penman-Monteith (1998)** — daily PET (ETo) from minimum and maximum
  air temperature, humidity, wind speed, solar radiation, and elevation,
  following Equation 6 in Allen et al. (1998) FAO Irrigation and Drainage
  Paper 56.

## Public API

Use `climate_indices.pet_thornthwaite()`, `climate_indices.pet_hargreaves()`,
and `climate_indices.pet_penman_monteith()`. All three accept NumPy arrays and
xarray DataArrays; the NumPy path is stable and the xarray path is beta.

## FAO-56 Penman-Monteith

`climate_indices.pet_penman_monteith()` implements the full FAO-56 reference
crop evapotranspiration equation (Eq 6):

$$
ET_0 = \frac{0.408 \Delta (R_n - G) + \gamma \frac{900}{T+273} u_2 (e_s - e_a)}
{\Delta + \gamma (1 + 0.34 u_2)}
$$

The FAO-56 intermediate variables are derived internally from the supplied
meteorology rather than passed in:

- atmospheric pressure (Eq 7) and the psychrometric constant (Eq 8) from
  elevation;
- saturation vapour pressure (Eq 12) and its slope (Eq 13) from daily minimum
  and maximum temperature;
- actual vapour pressure (Eq 14-19) from the best available humidity pathway;
- wind speed at the 2 m standard height (Eq 47) from the measurement height;
- net radiation (Eq 40) from extraterrestrial radiation (Eq 21), clear-sky
  radiation (Eq 37), net shortwave radiation (Eq 38), and net longwave
  radiation (Eq 39) with the cloudiness term.

**Units and return value.** All temperatures are degrees Celsius, wind is
metres per second, radiation is MJ m-2 day-1, elevation is metres, and the
returned PET is in millimetres per day — the same convention as
`pet_hargreaves`.

**Required versus optional inputs.** The required inputs are daily minimum and
maximum temperature, latitude, elevation, wind speed, and (for NumPy input) the
day of year. Humidity and radiation inputs are optional and grouped into
`climate_indices.HumidityInputs` and `climate_indices.RadiationInputs`.
Humidity pathway precedence is dewpoint, then RHmin/RHmax, then RHmax, then RHmean, and finally
the arid-region `e0(Tmin - 2)` estimate. Radiation precedence is supplied solar
radiation, then sunshine hours, then the temperature-range estimate (Eq 50,
with the interior `kRs = 0.16` or coastal `kRs = 0.19` coefficient); the
temperature-range estimate is limited to the clear-sky radiation. A NumPy call
requires an explicit `day_of_year`; an xarray call infers it from the time
coordinate.

**Soil heat flux.** For daily steps `G` is assumed to be zero, matching the
FAO-56 daily convention; a caller may override it via
`soil_heat_flux_mj_m2_day`.

**Missing-data behaviour.** NaN inputs propagate to the output. An xarray input
must be a daily coordinate beginning on January 1, consistent with the other
daily PET adapter; optional time-series inputs are trimmed to the shared time
range and dropped timesteps emit `InputAlignmentWarning`. A shared cell
dimension with differing coordinates raises `CoordinateValidationError`
instead of being intersected.

**SPEI use.** SPEI accepts a caller-supplied PET array, so the output of
`pet_penman_monteith()` can be passed directly to `climate_indices.spei()` as
`pet_mm`.

## Validation

PET has no single authoritative external gridded dataset: comparing against
a reanalysis-driven reference-ET product (e.g. gridMET) confounds algorithm
correctness with genuine climate-dependent bias between PET formulas, so no
defensible numerical tolerance exists (see the decision on
[Identify a gridded reference-ET dataset for PET stretch validation](https://github.com/monocongo/climate_indices/issues/775)).
Literature worked examples are therefore the primary external output-reference
evidence for PET. The sources provide daylight hours or extraterrestrial
radiation rather than a latitude/date pair, so the tests reconstruct those
coordinates with the library's solar geometry; they do not independently
validate that reconstruction.

`tests/test_pm_eto.py` covers the FAO-56 Chapter 3 radiation and wind helpers
against worked Examples 10-16 (Allen et al., 1998) and the full daily chain
against Example 18 (Uccle, 6 July, ETo = 3.88 mm/day, at ±0.05 mm/day).
`tests/test_pm_eto_xarray.py` checks that the xarray path reproduces the NumPy
path, keeps Dask input lazy, lands the CF metadata, and enforces the daily
calendar contract.

`tests/test_eto.py` covers the temperature-based methods against synthetic
regression fixtures (`tests/fixture/temp_celsius.npy`,
`tests/fixture/pet_thornthwaite.npy`, and the `hargreaves_*` fixtures generated
in `tests/conftest.py`) and, since
issue [#774](https://github.com/monocongo/climate_indices/issues/774),
against literature-derived worked examples in
`tests/fixture/pet_literature/`:

- `test_eto_thornthwaite_literature_watson` — Thornthwaite worked example
  from Watson & Burnett (1995).
- `test_eto_hargreaves_literature_mehta` — Hargreaves-Samani worked example
  from Mehta (2006).

Neither the original Thornthwaite (1948) *Geographical Review* article nor
the original Hargreaves & Samani (1985) *Applied Engineering in Agriculture*
article is freely accessible from this environment (JSTOR paywall for the
former; the latter's hosts were unreachable). Both worked examples are
instead taken from secondary sources that reproduce the original equations
with full inputs and outputs — the same two worked examples used by PyETo
(https://github.com/woodcrafty/PyETo) to validate its own Thornthwaite and
Hargreaves implementations, which `climate_indices.eto` credits in its
module docstring as the code it was derived from. Full citations, and the
derivation of the (latitude, day-of-year) inputs each worked example
implies, are documented in `tests/fixture/pet_literature/metadata.json`.

## Tolerance

- Thornthwaite: `atol=4.0` mm/month, matching the tolerance PyETo's own test
  suite established for this exact source (its rounded coefficients and
  intermediate values, and no day-count adjustment, introduce a few mm of
  imprecision independent of any implementation difference). The worked
  example's monthly daylight hours are matched to an approximate
  back-solved latitude (43.0°N) since the source states daylight hours
  directly rather than a latitude; see metadata.json.
- Hargreaves: `atol=0.05` mm/day, matching the source's 1-decimal-place
  precision. The worked example's extraterrestrial radiation input is
  reproduced exactly via a back-solved (latitude, day-of-year) pair; see
  metadata.json.
- FAO-56 Penman-Monteith: `atol=0.05` mm/day, matching the rounded
  intermediate values printed in FAO-56 Example 18. The radiation and wind
  helper tests use a `0.05`-`0.1` absolute tolerance against the 1-decimal
  values printed in Examples 10-16.

## References

- Thornthwaite, C.W. (1948). An Approach toward a Rational Classification
  of Climate. Geographical Review, 38, 55-94.
  https://doi.org/10.2307/210739
- Hargreaves, G.H. and Samani, Z.A. (1985). Reference Crop Evapotranspiration
  from Temperature. Applied Engineering in Agriculture, 1(2), 96-99.
- Allen, R.G., Pereira, L.S., Raes, D., and Smith, M. (1998). Crop
  Evapotranspiration - Guidelines for Computing Crop Water Requirements.
  FAO Irrigation and Drainage Paper 56. https://www.fao.org/4/x0490e/x0490e00.htm
- Watson, I. and Burnett, A.D. (1995). Hydrology: An Environmental Approach.
  Lewis Publishers/CRC Press. ISBN 978-1-56670-087-0.
- Mehta, V.K. (2006). Estimating Evapotranspiration from Weather Data.
  Arghyam / Cornell University.
- Richards, M. PyETo. https://github.com/woodcrafty/PyETo

## docs/algorithm_refs/palmer.md

# Palmer Drought Indices

## Scope

The Palmer implementation computes five monthly drought products from
precipitation, PET, and available water capacity:

- Palmer Drought Severity Index (PDSI)
- Palmer Hydrological Drought Index (PHDI)
- Palmer Modified Drought Index (PMDI)
- Palmer Z-Index
- Self-calibrated Palmer Drought Severity Index (scPDSI)

The public NumPy APIs are `climate_indices.palmer.pdsi()` and
`climate_indices.palmer.scpdsi()`. `climate_indices.pdsi()` accepts labeled
xarray inputs and returns an `xr.Dataset` of the four standard outputs
(`pdsi`, `phdi`, `pmdi`, `z_index`) with CF metadata; scPDSI stays NumPy-only.
The 3.0.0 xarray notebook demonstrates the manual rewrap pattern for callers
who need it.

## Algorithm Notes

The implementation follows the water-balance structure described by Palmer
(1965):

1. Validate precipitation and PET arrays and available water capacity.
2. Initialize monthly water-balance terms.
3. Compute CAFEC coefficients for the calibration period.
4. Compute K weighting factors.
5. Calculate the Z-Index moisture anomaly.
6. Finish the recursive Palmer index calculations for PDSI, PHDI, and PMDI.

`scpdsi()` shares the water-balance and CAFEC stages, then computes
self-calibrating K-prime factors, fits location-specific duration factors,
uses the Wells recursion, and performs three cumulative Z-index rescaling
passes over the requested calibration period.

### Exact-zero comparisons

Both recursions branch on exact comparisons against zero, and every float
equality comparison site in `src/climate_indices/palmer.py` and the pdi.f
recursion in `src/climate_indices/_palmer_pdi.py` carries a rationale comment
plus a `# NOSONAR` marker. Those zeros are state sentinels, not
near-equality: the statements that clear a spell, a candidate, or a calibration
sum assign `0.0` exactly. A tolerance would change the spell state machine and
the committed reference fixtures, so no tolerance is applied to these
comparisons.

| Compared value | Exact zero means | Assigned exact `0.0` by |
| --- | --- | --- |
| `px3`, read as `x3` inside `_case` | no established spell: `_case` selects the near-normal (larger-magnitude incipient) value, PHDI falls back to the PDSI value, and an incipient index may be promoted | the spell-end assignments, `_statement_190`'s `ppr >= 100` clamp, and the initial `0.0` |
| `px1` / `px2` | no incipient wet/dry index exists to promote | the `px1_computed > 0` / `px2_computed < 0` clamps, the post-promotion and `_statement_210` resets, and the initial `0.0` |
| `sx1` / `sx2` | the backtrack trail has no candidate from that index at this step | `_statement_200`'s deferral copy of the clamped-to-zero `px1`/`px2`; the initial zeros only fill rows the `i < k8` guard masks out |
| `pro` / `ppr`, at 100 and 0 | the spell certainly ends, or no abatement is underway | `ppr` clamped to `100.0` when `>= 100`, reset to `0.0` when a spell ends, reset by `_statement_210` when a spell intensifies or an abatement fizzles, and the initial `0.0` |
| CAFEC `numerator` / `denominator` | the month accumulated no water-balance term, so the ratio is undefined | the water-balance accumulators' initial `0.0`: `both_zero` when both sums vanish (`1.0` for alpha/beta/gamma, `0.0` for delta), `0.0` when only the denominator does |
| `k8` | no months are pending a spell flush (integer counter, not a float) | the recursion's own counter |

Exact-zero tests outside these sentinels are guards on a computed divisor, not
sentinels, and raise rather than dividing: `DurationFactors.from_fitted`'s
cross-denominator check and `_palmer_wells._abatement_transition`'s `q` guard
(both `ConvergenceError`), and `_palmer_pdi._weighting_fraction`'s divisor
check (`ValueError`).

## Validation

`tests/test_palmer.py` compares PDSI, PHDI, PMDI, Z-Index, and CAFEC
coefficients against the committed fixture set for climate divisions. The test
module is marked `validation` and uses `atol=5e-5`, `rtol=0`. That fixture set
is regression coverage only: its provenance states the expected outputs were
generated from this library's own Palmer implementation.

For the standard family (`pdsi()`), qualified independent external-product
validation comes from the NOAA NCEI nClimDiv comparison in
`tests/test_nclimdiv_reference.py`: all 344 climate divisions, January 1895
through December 2022, median absolute difference 0.0127 with 86.2% of months
within 0.05. The comparison is qualified because NCEI publishes two decimal
places and the committed precipitation/PET inputs are not byte-identical to
NCEI's operational inputs, which sets a practical floor near 0.013.

`tests/test_scpdsi.py` compares the four self-calibrating index outputs and
fitted duration factors against committed Wells-lineage reference fixtures for
344 climate divisions at the same tolerance. That comparison is independent
implementation cross-validation. The nClimDiv comparison for `scpdsi()` in
`tests/test_nclimdiv_reference.py` is characterization only, because nClimDiv
applies standard Palmer's fixed national calibration; the NCEI Fortran and
Wells C++ recursion lineages disagree on roughly 6.3% of months by up to 8
PDSI units (see `tests/fixture/palmer/provenance.json`).

## Release Decision

Independent external-product validation for the standard Palmer family is
satisfied by the nClimDiv comparison above. scPDSI has independent
implementation cross-validation (Wells oracle) but no qualified
external-product reference yet; assessing DRI/WRCC scPDSI for that role remains
open (issue #780).

## References

- Palmer, W. C. (1965). Meteorological Drought. U.S. Weather Bureau Research
  Paper No. 45.
  https://www.droughtmanagement.info/literature/USWB_Meteorological_Drought_1965.pdf
- Heddinghaus, T. R., and Sabol, P. (1991). A review of the Palmer Drought
  Severity Index and where do we go from here? Preprints, 7th Conference on
  Applied Climatology, American Meteorological Society.
- Wells, N., Goddard, S., and Hayes, M. J. (2004). A self-calibrating Palmer
  Drought Severity Index. Journal of Climate, 17, 2335-2351.
