# 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, flood potential, and
climate-driven energy supply and demand.

### 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()` and `flood_event_summary()`, which extract events with onset, duration,
  peak, and severity from a daily flood index and summarize them per period or hydrological year

**[Energy](https://climate-indices.readthedocs.io/en/latest/energy_applications.html)** (`climate_indices.energy`)

The degree-day, wind-profile, and standardized energy index APIs have stable
`xarray.DataArray` support from 3.1.0. The DataArray route for
`energy_drought_events()` remains Beta. Available functions include:

- Heating, cooling, and energy degree days (HDD / CDD / EDD). Validation: exact reference for
  the published cooling example; specification-level for heating and energy degree days
- Hub-height wind-speed profile. Validation: independent-implementation comparison on
  synthetic inputs, not observed-wind validation
- Wind capacity factor (NumPy-only). Validation: regression only; specification-level
- Temperature-derated PV capacity factor (NumPy-only). Validation: independent-implementation
  comparison for the power/temperature term; specification-level for the thermal model and clipping
- 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). Validation:
  independent-implementation comparison against R `SEI` 0.2.0
- Run-theory energy-drought events (onset, duration, severity). Validation: regression only;
  specification-level
- A documented recipe for a capacity-weighted precipitation proxy for hydropower
- These indices describe the climate-driven component of energy supply and demand, not
  actual generation or load

**[Compound hazards](https://climate-indices.readthedocs.io/en/latest/compound-hazards.html)** (`climate_indices.compound`)

- Declared-hazard harmonization onto a shared calendar, joint counts and frequencies, the
  likelihood multiplication factor, conditional probabilities, concurrent and sequential
  event matching, and a monthly block bootstrap
- The NumPy routes are stable; the Dataset routes and the eager event catalogs are Beta
  ([ADR-0022](docs/adr/0022-compound-analysis-contracts.md))
- These statistics describe co-occurrence; they make no causal or significance claim

**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,hdd,cdd,edd,wind_profile,srepi,demand_index,srli,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.1.0, which is not yet published on PyPI. The published 3.0.0
release added the `fire`, `flood`, and `runs` modules and the optional Rust backend.
Relative to 3.0.0, 3.1.0 adds the `compound` module and the `energy` family, promotes the
DataArray index APIs to Stable, and adds optional Icechunk I/O.

Install the current published release from PyPI:

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

`uv` users can run `uv pip install "climate-indices>=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).

Starting with 3.1.0, the optional `icechunk` extra adds Icechunk-backed reads and
writes (`icechunk>=1.1.21,<2`). It needs Python 3.11 or later; the core install is
unchanged. Until 3.1.0 is published, install the extra from a 3.1.0 source checkout:

```bash
pip install ".[icechunk]"
```

After 3.1.0 is published, install it from PyPI:

```bash
pip install "climate_indices[icechunk]>=3.1.0"
```

See [Write NetCDF, Zarr, and Icechunk outputs](docs/writing-outputs.md#read-and-write-icechunk-snapshots)
and the [Icechunk SPI/SPEI notebook](notebooks/icechunk_spi_spei.ipynb).

## Rust crate

The numerical kernels are also available in the standalone `climate-indices`
crate's `0.1.0` release
([crates.io](https://crates.io/crates/climate-indices),
[docs.rs](https://docs.rs/climate-indices)):
`cargo add climate-indices`. The earlier `0.0.0` package reserved the name and
contains no algorithms. The crate is pure Rust: `ndarray` is its only
dependency, and it has no PyO3 or Python bindings. It provides kernels, not complete
index workflows, so the caller prepares calendars, selects the Calibration Period,
validates physical ranges, and applies the final transform.

The crate follows its own `0.x` version track, independent of the Python distribution
version, and needs Rust 1.85 or later. Before 1.0, a minor version can make breaking
changes. See [`crates/climate-indices/README.md`](crates/climate-indices/README.md) for
usage, examples, and benchmarks, and
[`crates/climate-indices/CHANGELOG.md`](crates/climate-indices/CHANGELOG.md) for the
published versions.

## 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`, the `fire`, `flood`, `energy`, and `compound` NumPy routes) | **Stable** | No breaking changes in minor versions |
| xarray index entry points accepting `xr.DataArray` (including Dataset-returning index and diagnostic functions) | **Stable from 3.1.0** | No breaking changes in minor versions |
| Compound xarray Dataset routes and eager event catalogs | **Beta** | No breaking changes in patch versions |

**Stable API**: The NumPy computation functions and xarray index APIs follow strict semantic versioning.
The xarray adapter provides automatic parameter inference, coordinate preservation, CF metadata, and Dask support.
Its interface is stable from 3.1.0, after the 3.0.0 calendar-alignment release
([ADR-0012](docs/adr/0012-xarray-api-stays-beta-through-3.0.0.md)).
Results follow the documented [numerical-equivalence bounds](docs/performance.md#numerical-equivalence).
This promotion changes no numerical algorithm or calendar behavior.

**Beta API**: The compound xarray Dataset routes and eager event catalogs remain Beta under
[ADR-0022](docs/adr/0022-compound-analysis-contracts.md).
Their interfaces may change in minor releases, but not patch releases.
They continue to emit `BetaFeatureWarning`; its public imports and warning filters remain supported.
Stable DataArray index APIs do not emit this warning.

See `docs/xarray_compatibility.md` for the 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 evidence for the library, including the 3.0.0 release artifacts and the 3.1.0 development version.
It separates independent external validation, same-author reference reproduction, synthetic characterization, and regression coverage.
The 3.1.0 additions are the compound-hazard module and the energy family.
Neither has an adopted product-level external oracle.
The compound evidence is a same-author published-method reproduction, a synthetic bootstrap characterization, and interface equivalence.
The energy evidence includes an exact published example, independent-implementation comparisons on synthetic and real-series inputs, and specification-level rows.
Their scope limits are in the sections below.

## 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, the NOAA CPC degree-day worked example and the xclim/windpowerlib wind-profile fixtures (independent implementations on synthetic inputs), the R `SEI` empirical-standardization and energy-wrapper comparison, and the flood-family source-identity checks (specification-level, not external; see "Flood Evidence Classification"), pass; the nine committed compound summer cases reproduce the same-author counts and LMF without independent-validation or significance claims; 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 notebooks/icechunk_spi_spei.ipynb` | All four notebooks execute from a clean checkout; the Icechunk notebook uses the optional extra that the `dev` group installs on Python 3.11 or later. |
| 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 events and per-period summaries (`flood.flood_events`, `flood.flood_event_summary`) | Regression only; specification-level (no external oracle adopted) | Hand-constructed events exact at `atol=0`; NumPy and per-cell xarray agreement at `rtol=0`; the strict `index > threshold` boundary is frozen | `tests/test_flood_events.py` checks onset, end, duration, peak, and raw-sum severity against a plain day-by-day reference over seeded block input, NaN, masked, infinite, empty, and boundary layouts; `tests/test_flood_event_summary.py` checks per-period counts, durations, maxima, and frequency against closed-form expectations, onset-period assignment, hydrological-year labels, and per-cell agreement between the NumPy and xarray paths. | Return period `T` is not implemented: the Deo et al. (2015) estimator is not recoverable, so `frequency` is an event rate and `docs/flood_applications.md` gives a caller-side ranked-severity recipe. The paper's Brisbane values (`I_F_acc = 118`, `I_F_max = 4.4`, `D_F = 104`, `T = 106.2 years`) are **not** claimed as reproduced; no input record or estimator is available (ADR-0014 decision 5). |
| Compound harmonization (`compound.harmonize`, `AlignedHazards`) | Regression only; specification-level (no external oracle adopted) | Exact Boolean and integer equality; the outputs are masks and counts, so no numeric tolerance applies | `tests/test_compound_harmonize.py` checks hand-counted strict-threshold ties and endpoints, float32 widening, month-label normalization, real leap days and interior absent dates, `exact` and `inner` alignment, and daily `any`/`all` reduction (including the complete-month rule) against an independent calendar loop. It also checks rejection of raw z-scores, out-of-range values, unsupported calendars, and conflicting Calibration Periods without clipping or interpolation, plus read-only results and the serializable manifest. These tests do not use the `validation` marker. | The declared-indicator route participates in the scoped summer reproduction below. This does not validate probability-input fitting, general harmonization, or daily-to-monthly reduction. A declared probability is a trust boundary and is not proven to be a PIT value or a good fit. Agreement of the beta xarray route (#1372) with this function is interface equivalence, not validation. |
| Compound joint statistics (`compound.joint_statistics`, `JointResult`) | Scoped reference-reproduced, not independently validated, for the nine public summer cases below; otherwise analytic/specification-level regression | Exact integer equality for counts and exact equality for the balanced Cartesian tables; `atol=1e-12` with `rtol=0` for closed-form probabilities; the Gaussian cases are predeclared 5-sigma Monte Carlo bands, not tolerances | `tests/test_compound_statistics.py` checks pooled and monthly counts against hand-counted contingency tables, one complete-case support for every marginal with different per-source validity masks (two and three sources), balanced two- and three-source independence (`LMF=1` exactly), the closed-form Gaussian joint probability `1/4 + asin(rho)/(2*pi)` at `rho=0` and `rho=0.5` within its declared band, disjoint hazards (`LMF=0`), a zero marginal (`NaN` plus the `undefined_denominator` flag), short support, sparse occurrences, all-missing grid cells, invalid steps never counting as non-exceedances, and pooled-versus-monthly recounting. Hypothesis tests cover source-permutation symmetry, count nesting, and probability bounds on generated validity masks. It also checks option and manual-`AlignedHazards` validation, read-only results, the JSON manifest, and lifecycle logging. These tests do not use the `validation` marker. | The analytic values come from closed form or hand counting. The separate public-study comparison below reproduces the same-author counts and LMF, not an independent oracle. Agreement of the beta xarray route (#1372) with this function is interface equivalence. `LMF=1` in a balanced design is a property of that design, not evidence of independence in real data. This descriptive function makes no inference about significance; see the separate bootstrap evidence below. Known gaps: the plan's nonzero-threshold quadrature and independent-uniform-envelope checks are not implemented, outside the summer reference scope. The nine public cases do not validate other thresholds, three-hazard real-data relationships, upstream SPEI, or the published significance mask. |
| Compound conditional statistics (`compound.conditional_probability`, `ConditionalResult`) | Regression only; specification-level (no external oracle adopted) | Exact integer equality for counts and exact equality for hand-enumerated windows; `atol=1e-12` with `rtol=0` for closed-form frequencies | `tests/test_compound_statistics.py` checks the `(0, 0)` identity against the joint and marginal counts, hand-enumerated forward windows, one counted success per given-positive anchor despite several target occurrences in its window, discarded incomplete trailing horizons (including a window that already holds a target occurrence), a discarded anchor whose window holds an unknown target step, asymmetry under a role swap, the response-window baseline against the one-step marginal and the joint probability, anchor-month grouping across a year boundary, daily civil-day steps across 29 February, a zero conditioning count (`NaN` plus the `undefined_denominator` flag), short support, all-missing grid cells, and a third source that never enters the selected pair's eligibility. Hypothesis tests cover count nesting and probability bounds on generated validity masks. It also checks option, role, and window validation, read-only results, the JSON manifest, and lifecycle logging. These tests do not use the `validation` marker. | No external oracle is used: the analytic values come from closed form or hand enumeration, so this is specification-level evidence, not external validation. A lagged association is not a validated causal cascade and is not an event-onset conversion probability. Agreement of the beta xarray route (#1372) with this function is interface equivalence, and this descriptive function makes no significance or confidence claim; see the separate bootstrap evidence below. The public LMF reproduction below does not validate conditional response windows. |
| Compound events (`compound.compound_events`, `CompoundEventSet`) | Regression only; specification-level (no external oracle adopted) | Exact index, count, lag, and flag equality; the outputs are indices, counts, lags, and flags (NaN where a lag does not apply), so no numeric tolerance applies | `tests/test_compound_events.py` checks the all-pairs concurrent and sequential enumeration, the four boundary-censor flags, and the rejected-bridge count against an independent brute-force interval scan over 24 seeded random scenarios with random unknown steps, plus hand-counted inclusive-end cases, repeated matches, lag-window boundaries, sequential bridges through an unknown step (rejected and counted), a matched and a rejected candidate in one window, empty catalogs with stable dtypes, validation of manually constructed `RunSet` objects, read-only result arrays, lifecycle logging, and the serializable manifest. These tests do not use the `validation` marker. | No external oracle and no published reference values are committed, and this specification-level status stands until an independent matching oracle exists. The seeded brute-force scan is an internal specification oracle, not external validation. Pair counts are not independent observations, so no significance claim is made or supported here. Agreement with the eager xarray catalog (#1373) is interface equivalence, not validation. |
| Compound monthly bootstrap (`compound.bootstrap_dependence`, `BootstrapResult`) | Synthetic characterization and regression only; screening evidence measured group by group, not as one configuration: at the default 1,999 resamples, 200-year records with the seasonal frequencies (0.35 to 0.55) and 1, 6, or 36 months of persistence with independent hazards; at 399 resamples, 200-year records with 1 to 30 months of persistence, monthly frequencies of 0.30 to 0.60 or 0.70 to 0.80, conditional lag windows to `(1, 3)`, 0.5-correlated hazards at 1, 6, or 12 months of persistence, and 100-year records; 36-month persistence at 399 resamples with independent hazards and at the default resamples with correlated hazards, 60- and 30-year records, and frequencies near 0.20 each failed a predeclared gate | Exact counts and seeded repeats; `atol=1e-12`, `rtol=0` for closed-form frequencies and linear quantiles. Characterization uses predeclared exact-binomial acceptance counts, not floating-point tolerances. | `tests/test_compound_bootstrap.py` covers paired versus independent sampling, tiny-block enumeration, month phase, original-chronology lag responses, missing-year segments, earliest ties, separate interval/null guards, zero-joint suppression, inclusive null ties, valid-resample fractions, immutable outputs, lifecycle events, and cell-tile invariance. `scripts/characterize_compound_bootstrap.py` records 128 independent Gaussian records per scenario over 183 cases: the four original 200-year scenarios, record lengths of 30, 60, and 100 years, persistence of 12 to 36 months, conditional lag windows `(0, 1)`, `(0, 3)`, and `(1, 3)`, frequency profiles from 0.20 to 0.80, correlated hazards, and 399 and 1,999 resamples. A 40,000-year simulation supplies the lagged-conditional truth, and a test checks that reference against quadrature at zero lag. The measured report is `docs/design/compound-bootstrap-characterization.json`. | Supported scope (maintainer decision, 2026-10-19), grouped by the resample count that measured each condition because the plan crossed no two of them: at the default 1,999 resamples, 200-year monthly records with the seasonal frequencies (0.35 to 0.55), source overlap of 1, 6, or 36 months, and independent hazards; at 399 resamples, 200-year monthly records with 1 to 30 months of overlap, lag windows `(0, 0)`, `(0, 1)`, `(0, 3)`, and `(1, 3)`, frequencies of 0.30 to 0.60 or 0.70 to 0.80, and hazards correlated 0.5 at the same time with 1, 6, or 12 months of overlap, plus 100-year records with the seasonal frequencies. Unsupported: 36-month persistence at 399 resamples for independent hazards and at the default resamples for correlated hazards, frequencies near 0.20, 60- and 30-year records, and `wide` at 36-month overlap; each failed a predeclared gate when the guards answered, except `wide` at 36 months, which was not run. No support guard was added, because no observable separates the passing from the failing cases without refusing a supported one. No independent-hazard case exceeded its per-month null rejection band, and the 36-month pooled null rejection falls from 7.5% at 399 resamples to 3.5% (3.2% to 3.8%) at the default resamples. No external validation, daily inference, field significance, or arbitrary-persistence guarantee is claimed. See the Compound Bootstrap Characterization section below. |
| Compound beta xarray routes (`compound.harmonize`, `joint_statistics`, `conditional_probability`, `bootstrap_dependence` on `xarray` input) | Regression only; interface equivalence with the NumPy routes (no external oracle adopted) | Exact equality for every count, flag, Boolean, and NaN position; the estimates come from the same NumPy arithmetic on the same counts, so no numeric tolerance applies | `tests/test_compound_xarray.py` compares every result field of each route with its NumPy counterpart, per month and pooled, for in-memory and Dask input over several spatial chunkings; checks bootstrap results are identical across tile boundaries; proves graph construction runs no task under a scheduler that raises; and checks custom and transposed dimensions, NaN versus explicit validity masks, daily-to-monthly reduction, `join="inner"` with its warning, a single location, and the rejected layouts (split time or hazard chunks, mismatched spatial chunks, mixed backends, missing or shifted coordinates, dimension-only or non-Gregorian time, and options that apply only to NumPy input). It also checks registry attributes, bounds, flag masks, no inherited `standard_name` or bounds, the versioned JSON manifest, h5netcdf and Zarr round trips, the `BetaFeatureWarning` and replayed-warning policy, and lifecycle events. `scripts/measure_compound_xarray_memory.py` records per-tile memory in `docs/design/compound-xarray-memory.json`. These tests do not use the `validation` marker. | Equivalence shows the Dataset routes reproduce the NumPy core, not that either is scientifically validated; the NumPy rows above state that evidence. Value-dependent support warnings are not emitted for Dask input, so read `quality_flags`. Memory figures are `tracemalloc` lower bounds from one machine, one platform, and one daily and one monthly configuration, with a manual peak-RSS spot check for `harmonize` only, and bootstrap tile width is set by compute time, which was not characterized. The route is beta (ADR-0012). |
| HDD, CDD, EDD (`energy.heating_degree_days`, `cooling_degree_days`, `energy_degree_days`) | Exact reference for the published cooling example; specification-level for heating and energy degree days | `rtol=0`, `atol=1e-12` °C day after the `(F - 32) * 5/9` conversion | `tests/test_energy_reference.py` (marked `validation`) reproduces NOAA CPC's published seven-day CDD series (`[67, 65, 70, 74, 78, 65, 68]` °F to `[2, 0, 5, 9, 13, 0, 3]` °F day, total 32 °F day) with the exact `(65 - 32) * 5/9` °C base. Heating, EDD, and missing-day cases are hand-calculated specification checks; EDD = HDD + CDD is an algebraic identity, not an independent oracle. `tests/test_energy_degree_days.py` repeats the CPC example outside the marker and covers input validation, dtypes, missing values, and scalar/empty shapes; `tests/test_energy_xarray.py` covers temperature-unit inference. | No public default base ships: C3S and Eurostat pair an activation threshold with a separate reference temperature, which one base per transform cannot reproduce. The functions take a caller-supplied daily mean temperature, so a caller holding Tmin and Tmax computes `(Tmin + Tmax) / 2` first; no other averaging or weighting is implemented, and no C3S numeric comparison is committed. |
| Hub-height wind-speed profile (`energy.hub_height_wind_speed`, `energy.wind_shear_exponent`) | Independent-implementation comparison on synthetic inputs; not observed-wind validation | `rtol=0`, `atol=1e-12` m s⁻¹; two-level shear recovery `atol=1e-15` as consistency evidence only | `tests/test_energy_reference.py` (marked `validation`) compares the power-law and log-law profiles against pinned `xclim` 0.62.0 `wind_profile` and `windpowerlib` 0.2.2 `logarithmic_profile` outputs for nine synthetic cases, including calm, identity, and missing masks. `tests/test_energy_wind.py` and `tests/test_energy_xarray.py` cover the contracts outside the marker. | The inputs are synthetic and heights are generic, so neither comparison establishes skill against observed hub-height wind. Calm handling and the default one-seventh exponent are project conventions. |
| Wind capacity factor (`energy.wind_capacity_factor`) | Regression only; specification-level (no external numeric curve adopted) | Exact rational arithmetic `rtol=0`, `atol=1e-15` for the cubic curve; exact region boundaries, table knots, and hand-calculated linear examples | `tests/test_energy_capacity_factor.py` checks the `generic` and user-supplied curves against `fractions.Fraction` arithmetic, region membership, cut-in and cut-out, air-density scaling, and masking. The `wind_capacity_factor` docstring states the apply-before-averaging rule, but no test exercises it. `docs/research/energy-oracle-survey.md` found no matched numerical C3S curve to adopt. | No turbine-specific curve, observed generation, or model-to-model tolerance is validated. A curve of mean wind is not the mean capacity factor. |
| PV capacity factor (`energy.pv_capacity_factor`) | Independent-implementation comparison for the power/temperature term; specification-level for the thermal model and clipping | `rtol=0`, `atol=1e-12` in dimensionless capacity factor | `tests/test_energy_solar.py` compares nine synthetic cases against `pvlib` 0.15.0 `pvsystem.pvwatts_dc` (pinned commit and SHA-256 in `tests/fixture/energy_solar/`), with the cell temperature supplied through the inverse Jerez thermal equation. The thermal coefficients have hand-calculated tests. | The Jerez thermal model itself, irradiance transposition, optical and inverter losses, and plant output are not validated. Clipping to `[0, 1]` is project policy because the published potential can exceed one. |
| SREPI and standardized demand index (`energy.srepi`, `energy.standardized_demand_index`) | Independent-implementation comparison against R `SEI` 0.2.0; the empirical plotting position, tie rule, and minimum-observation count are fixed to `SEI`'s convention | `scaled_atol=1e-12`, `probability_atol=1e-15`, `normal_atol=1e-12` | `tests/test_energy_sei_reference.py` (marked `validation`) compares the empirical standardization against `SEI` 0.2.0 on the Bloomfield and Brayshaw (2021) German onshore wind capacity-factor series (852 monthly means, 1950-2020, calibration 1981-2010) in two cases and on synthetic capacity, production, demand, and residual arrays in ten cases; it then compares `srepi()` on the capacity, production, and reading cases and `standardized_demand_index()` on the synthetic demand cases. A gamma control is parameter-matched to `SEI`'s MLE values. Fixture provenance, checksums, and session information are in `tests/fixture/energy_sei/`. | The comparison is monthly, so no seasonal-diurnal or hourly calibration is validated against an independent implementation (`#1177` ships with specification tests only); only the two `reading` cases use provider data, so demand, production, and residual coverage is synthetic. The two real-series cases calibrate on 30 values per calendar month, below the 100-value empirical recommendation in `tests/fixture/energy_sei/warnings.txt`, so their probabilities are multiples of 1/32. The gamma control does not validate the library's gamma estimator. No C3S energy-indicator product comparison is committed, the beta distribution has no external reference distribution, and maintainer confirmation of fixture extraction and redistribution remains outstanding. |
| SRLI and residual load (`energy.srli`, `energy.residual_load`) | Independent-implementation comparison against R `SEI` 0.2.0; specification-level for the subtraction | The same `SEI` tolerances as SREPI | `tests/test_energy_sei_reference.py` (marked `validation`) compares `srli()` against `SEI` on the signed series, `tests/test_energy_srli.py` covers the NumPy contracts and negative retention, and `tests/test_energy_standardized_xarray.py` covers the xarray route and the empirical-only distribution restriction. | The committed residual cases are synthetic `demand - production` arrays, so no energy-system demand, capacity, or dispatch model is validated, and no C3S demand comparison is committed. |
| Energy-drought events (`energy.energy_drought_events`) | Regression only; specification-level (no external event catalogue adopted) | Hand-enumerated windows and exact integer counts; NumPy, eager, and Dask summaries compared exactly | `tests/test_energy_events.py` checks run metrics and period summaries on SPI, SREPI, and SRLI series, threshold-relative magnitude, missing steps, period label counts, and gridded Dask summaries that keep the time axis in one chunk. `docs/research/energy-oracle-survey.md` records that `SEI`'s event loop is not a missing-value-policy oracle. | No observed event catalogue validates onset, duration, or severity, and the `SEI` magnitude identity is a qualification, not parity. |

## Compound Published-Method Reproduction

`tests/test_compound_reference_fixture.py::test_compound_api_reproduces_public_summer_reference`
provides nine offline `validation` cases against `compound.harmonize()` and `compound.joint_statistics()`.
The reference is the same-author deterministic calculation in Vieira Passos et al. (2024),
*Identifying regional hotspots of heatwaves, droughts, floods, and their co-occurrences*,
*Stochastic Environmental Research and Risk Assessment* 38, 3875–3893,
<https://doi.org/10.1007/s00477-024-02783-3>.
The public inputs are Vieira Passos (2024), Figshare version 2,
<https://doi.org/10.6084/m9.figshare.25189118.v2>, under CC BY 4.0.

`tests/fixture/compound/` retains source values, missingness, full coordinates, checksums, file IDs, and the adapted code's MIT notice.
The reference pins author commit `6c88f255ee3034279dcebe2a2fc40c954dbe9575`.
The preparation uses month-end forward filling, followed by linear interpolation of remaining daily gaps, outside the production harmonizer.
The API receives declared Boolean indicators with `HWI >= 0` and `SPEI <= -1`, separate validity, and recorded comparators.
The original source Calibration Period is unknown.

The comparison covers June, July, and August separately over 1922–2021 at the three selected Swedish cells.
The measured run gave exact counts and zero maximum absolute error for joint probability, independence product, and LMF across all nine cases.
Counts require exact equality; floating-point comparisons use `rtol=0`, `atol=1e-12`.
These measured errors do not establish cross-platform bitwise identity.

| Cell | Month | Eligible days | Heat days | Drought days | Joint days | LMF, rounded |
| --- | --- | --- | --- | --- | --- | --- |
| Southern | June | 3,000 | 232 | 566 | 43 | 0.982393079079 |
| Southern | July | 3,100 | 530 | 469 | 68 | 0.848050850867 |
| Southern | August | 3,100 | 303 | 590 | 63 | 1.092465178721 |
| Central | June | 3,000 | 178 | 509 | 26 | 0.860908147723 |
| Central | July | 3,100 | 340 | 494 | 45 | 0.830554894022 |
| Central | August | 3,100 | 123 | 435 | 26 | 1.506401270909 |
| Northern | June | 3,000 | 136 | 743 | 36 | 1.068798986620 |
| Northern | July | 3,100 | 264 | 557 | 48 | 1.011914476906 |
| Northern | August | 3,100 | 66 | 523 | 19 | 1.706356104062 |

**Classification: reference-reproduced, not independently validated.**
The author preparation expands monthly drought values to days; it does not create independently observed daily drought resolution.
Inclusive source thresholds intentionally differ from the library's strict probability comparisons.
The author's zero-joint LMF is undefined, while the library returns zero with positive marginals; every empirical case has positive joint counts.
The checks preserve this intentional difference explicitly.
The nine cases do not reproduce the full Figure 8 map, its significance mask, or upstream SPEI correctness.
The author's IID significance procedure and random tie noise are excluded.
Daily inference remains unsupported, so this comparison makes no new bootstrap-significance claim.

| Route | Final evidence scope for this comparison |
| --- | --- |
| Prepared Boolean indicators and joint statistics | Same-author reference reproduction, restricted to the nine cases above. General harmonization remains specification-level. |
| Conditional statistics and event matching | Analytic/specification-level regression; the LMF reference does not validate these routes. |
| Monthly bootstrap | Synthetic characterization and regression only, with the limits below. Further calibration belongs to #1398. |
| Beta xarray routes and catalogs | Interface equivalence with the NumPy core, not an independent scientific oracle. |

The runnable [compound guide](docs/compound-hazards.md) keeps the reproduction and monthly synthetic inference examples separate.
No fixture download, missing-fixture skip, or optional external package substitutes for this reproduction gate.

## Compound Event Catalog Evidence

`compound.compound_events_xarray()` has regression and interface-equivalence evidence
only. `tests/test_compound_catalog.py` compares each event field, timestamp, censor
flag, support count, retained-run count, pair count, and rejected-bridge count exactly
with the existing NumPy matcher. The cases cover daily and monthly inputs,
concurrent and sequential modes, bounded and unbounded lags, duration filtering,
transposed spatial dimensions, one series, and an irrelevant third hazard.
The tests also cover shared missingness, all-invalid cells, no matches, zero-cell
subsets, rejected malformed inputs, and Dask rejection under a scheduler that fails
if any task runs. It also checks the 1,000,000-row cap, a duplicate hazard label, a
non-finite manifest, and dropped boundary attributes on the cell coordinates. Empty and
populated catalogs preserve values, cell labels, flags,
and versioned JSON metadata through h5netcdf and Zarr storage.

These tests do not use the `validation` marker. They establish no new scientific
validation, external oracle, significance calibration, or causal interpretation.
Pair counts are not independent observations. The beta catalog is explicitly eager;
its memory use and runtime have not been characterized. A selected input can be large
enough to exhaust memory. A match that would exceed the 1,000,000-row catalog cap raises
`InvalidArgumentError` before the rows that would overflow it are allocated. The final
public-study reproduction above concerns frequency counts, not event catalogs.

## Compound Bootstrap Characterization

### The #1371 experiment

The synthetic experiment predeclared a familywise error budget of 0.01 across
864 binomial comparisons. Each coverage case required at least 109 covered
records out of 128. Each independent-null case allowed at most 19 rejections
at nominal alpha 0.05. Missing intervals counted as uncovered; missing p-values
failed the null gate. Wilson 95% intervals accompany all measured proportions.

The default six-year blocks passed these broad screening checks. Coverage
counts ranged from 111 to 125 out of 128 across the default cases. This is not
proof of nominal 95% coverage. The twelve-year LMF interval for February in the
independent 36-month moving-average scenario covered only 106 of 128 truths.
Its estimated coverage was 82.8%, with a Wilson 95% interval of 75.3% to 88.4%.
The characterization therefore failed. The case is retained as failed: it was
not discarded and the acceptance band was not widened.

Block-length changes also changed individual test decisions. The report retains
those counts and every failed case.

The reported run follows a fix to p-value tie handling: the first run compared
nulls to the observed statistic exactly, so floating-point rounding dropped some
ties. The rerun used the same seeds and configuration. Interval coverage is
identical, null rejection counts changed by one in 4 of 36 cases, and no gate
outcome changed.

### #1398 measurements (2026-10-19)

[Issue #1398](https://github.com/monocongo/climate_indices/issues/1398) measured 183
cases: the four #1371 scenarios again, plus record lengths of 30, 60, and 100 years,
persistence of 12, 18, 24, and 30 months, conditional lag windows `(0, 1)`, `(0, 3)`,
and `(1, 3)`, monthly frequencies of 0.20 to 0.30, 0.30 to 0.60, and 0.70 to 0.80,
hazards correlated 0.5 at 1, 6, and 12 months of persistence, and the default 1,999
resamples. The gates, bands, and scenario plan were committed before the run, and the
replication group reproduces the #1371 report field by field. The predeclared design,
results, and decision are in `docs/design/compound-bootstrap-calibration.md`.

47 of the 183 cases fail a predeclared gate. The failures fall into three mechanisms.
37 cases miss the interval coverage band, and 45 cases carry a month where the library
suppressed the p-value because fewer than 95% of resamples were valid; 35 cases do both.
**No independent-hazard case exceeded its per-month rejection band.** The pooled
rejection rate across the 24 independent-hazard 200-year scenarios is 3.5% to 7.5%.
At the default 1,999 resamples the 36-month independent scenario passes all nine gate
checks and its pooled rejection rate is 3.5% (3.2% to 3.8%), so most of the 7.5% excess
seen at 399 resamples was a resample-count artifact. Under-coverage is the remaining
limit, and it appears where few effective blocks or few effective occurrences remain.
Coverage passes at 100 years and fails at 30 and 60 years; it passes for 12, 18, 24,
and 30 months of persistence and fails at 36; it passes for monthly probabilities from
0.30 upward and fails near 0.20.

### Supported scope (maintainer decision, 2026-10-19)

This statement supersedes the narrower 2026-10-18 one. The maintainer kept the method
unchanged, added no support guard because no observable separates the passing from the
failing cases without refusing supported ones, and widened the claim only where every
predeclared gate passed:

- **Supported**, listed group by group because the plan crossed none of these
  conditions with each other:
  - At the default 1,999 resamples: 200-year monthly records with the seasonal
    frequencies (0.35 to 0.55) and source overlap of 1, 6, or 36 months, hazards
    independent of each other. The independent 36-month case passes every gate here.
  - At 399 resamples: 200-year monthly records with source overlap of 1 to 30 months;
    conditional lag windows `(0, 1)`, `(0, 3)`, and `(1, 3)` at 1- and 6-month overlap;
    monthly frequencies of 0.30 to 0.60 at 1- and 6-month overlap and 0.70 to 0.80 at 1-,
    6-, and 36-month overlap; hazards correlated 0.5 at the same time with 1, 6, or 12
    months of overlap, where rejection describes power rather than null calibration; and
    100-year records with the seasonal frequencies at 1- and 6-month overlap.
  - Every listed condition passed every gate in its group. Measured coverage is 91.5% to
    94.8%. This is screening evidence, not proof of nominal 95% coverage.
- **Unsupported.** Sources that persist for 36 months: with independent hazards the
  twelve-year LMF coverage gate fails at 399 resamples and passes at the default
  resamples, and with correlated hazards the default-resample case fails one coverage
  gate (3-year blocks, month 4, 107 of 128). Monthly probabilities near
  0.20, where the null p-value is suppressed on too many resamples and coverage also
  fails at 36-month persistence. Records of 60 years, which miss coverage in months 6
  to 8 (94 to 108 of 128) and suppress p-values in months 6 to 9. Records of 30 years,
  where the guards remove 23% of month-records at typical frequencies and 63% near
  0.20, and the answered records fail both gates. Profile `wide` at 36-month overlap,
  which was not run. Do not read an interval or p-value from these cases as nominal.
- **Never claimed.** No external validation, daily inference, field significance,
  nonstationarity, or long memory. The guards report refusal through `quality_flags`
  and effective resample counts; they do not detect the source structure itself.

The beta xarray adapter is #1372, with equivalence-only evidence recorded in the
compound xarray row above. The [worked guide](docs/compound-hazards.md) preserves these
limits. Improving calibration for 36-month persistence or 60-year records needs a
different method or support contract, and a reopen would need its own predeclared
characterization.

## 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
  DataArray index APIs, stable from 3.1.0.

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.

## Energy Evidence Classification

The energy family — heating, cooling, and energy degree days, the hub-height wind
profile, the wind and PV capacity factors, SREPI, the standardized demand
index, SRLI, and the energy-drought event statistics — adopts two classes of
external anchor and no product-level oracle:

- **Exact reference**: NOAA CPC's published seven-day CDD example is transcribed
  in `tests/fixture/energy_cpc/` with its retrieval date and source checksum. It
  is authoritative for that one example, not for the C3S threshold/reference
  scheme or for a default base.
- **Independent implementation**: R `SEI` 0.2.0 pins the empirical
  standardization on a real, attributed ERA5-derived capacity-factor series;
  `xclim` 0.62.0 and `windpowerlib` 0.2.2 pin the wind profile; `pvlib` 0.15.0
  pins the PV power/temperature term. These compare implementations on identical
  inputs, not skill against observations. The empirical seam implements `SEI`'s
  plotting position, tie rule, and minimum count by design, so its close
  agreement is expected rather than an independent check of the method.
- **Not adopted**: the C3S energy-indicator product. Its degree days are
  national monthly CSV, not a matched gridded series, and its capacity factors
  use turbine-specific curves and configurations the library does not reproduce.
  `docs/research/energy-oracle-survey.md` records the missing matched inputs; no
  C3S numeric comparison runs.
- **Not adopted**: beta and the energy-drought event statistics. The `SEI`
  fixture has a finite-bound control on `SEI`'s censored empirical PIT, not a
  beta reference distribution, so the library's beta route has no external
  reference; `SEI`'s event loop uses different missing-value and magnitude
  conventions.
- **Regression only**: the wind capacity factor, the event statistics, and the
  Jerez thermal model and clipping policy. Their tests encode published algebra
  and public contracts, so they are specification evidence.
- **Specification evidence only**: hourly standardization (`#1177`) ranks each
  hour of day x day of year stratum against itself, and its tests check the plotting
  position, the calendar plan, and agreement between the NumPy and xarray paths.
  There is no independent hourly oracle, because `SEI` accepts hourly timestamps
  but does not stratify by hour of day or day of year. The ERA5/C3S gridded
  demonstration remains pending (`#1179`).

Energy indices describe the climate-driven component of energy supply and
demand, not actual generation, load, or residual load. Even an exactly
reproduced standardized index would not be a generation or demand forecast.
`tests/fixture/energy_cpc/`, `tests/fixture/energy_sei/`, and
`tests/fixture/energy_wind/` record that maintainer confirmation of fixture
transcription or generation and redistribution is still outstanding.

## 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
- Allen, S., and Otero, N. (2023). Standardised indices to monitor energy
  droughts. *Renewable Energy*, 217, 119206.
  https://doi.org/10.1016/j.renene.2023.119206
- Bloomfield, H., and Brayshaw, D. (2021). ERA5 derived time series of European
  aggregated surface weather variables, wind power, and solar power capacity
  factors: hourly data from 1950-2020. University of Reading.
  https://doi.org/10.17864/1947.000321
- Copernicus Climate Change Service. Global climate and energy indicators from
  1950 to present derived from reanalysis.
  https://doi.org/10.24381/3bb607bd
- Jerez, S., Tobin, I., Vautard, R., Montávez, J. P., López-Romero, J. M.,
  Thais, F., Bartok, B., Christensen, O. B., Colette, A., Déqué, M., Nikulin,
  G., Kotlarski, S., van Meijgaard, E., Teichmann, C., and Wild, M. (2015).
  The impact of climate change on photovoltaic power generation in Europe.
  *Nature Communications*, 6, 10014. https://doi.org/10.1038/ncomms10014
- NOAA Climate Prediction Center. Degree days explained.
  https://www.cpc.ncep.noaa.gov/products/analysis_monitoring/cdus/degree_days/ddayexp.shtml
- 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

:::{note}
The DataArray index APIs and public adapter interface are **Stable** from 3.1.0
([ADR-0012](docs/adr/0012-xarray-api-stays-beta-through-3.0.0.md)).
They introduce no breaking interface changes in minor versions.
This includes index and diagnostic functions that return Datasets.
The promotion changes no numerical algorithm, calendar behavior, or Dask chunk requirement.
The compound Dataset routes and eager event catalogs remain **Beta** under
[ADR-0022](docs/adr/0022-compound-analysis-contracts.md) and continue to emit `BetaFeatureWarning`.
Its public imports and existing warning filters remain supported.
See the [API change record](docs/deprecations/api-changes.md#dataarray-index-api-promotion-310)
for the promotion and warning compatibility details.
:::

| 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 the generic standardized index | Yes | `climate_indices.standardized_index()` (also `climate_indices.indices.standardized_index()` for NumPy) standardizes a caller-supplied series under the selected distribution's support constraints, so it stamps the neutral `CF_METADATA["standardized_index"]` entry: `long_name` "Standardized index", `units` "1", no `standard_name`. Temporal parameters are inferred from the time coordinate, `aggregation` (`"sum"` or `"mean"`) reaches the fit and is recorded in the output's calculation metadata, and all four distribution families (gamma, Pearson Type III, empirical, beta) stay available. Beta requires source and windowed values in `[0, 1]`. Gamma and Pearson Type III clip negative values to zero by default and reject them with `clip_negatives=False`. Empirical preserves signed values with `clip_negatives=False`. Gridded input is standardized one cell at a time, so every cell matches its own NumPy series exactly, and `fitting_params` reaches every cell unchanged, so cell-specific parameters are not supported. Root re-export, aggregation reporting, and the neutral metadata follow [ADR-0020](docs/adr/0020-standardized-index-public-surface.md). See `tests/test_standardized_index_xarray.py`. |
| `DataArray` inputs for standardized energy indices | Yes | `energy.srepi()`, `energy.standardized_demand_index()`, and `energy.srli()` take the same adapter path as the generic seam, with their own registry entries. `srepi()` defaults to `aggregation="mean"`, which suits a capacity factor, so accumulated production must pass `aggregation="sum"` explicitly; the demand index and SRLI default to `"sum"`. A DataArray input keeps the NumPy wrapper's distribution restrictions: empirical or beta for SREPI and the demand index, empirical only for the signed SRLI series. See `tests/test_energy_standardized_xarray.py`. |
| `DataArray` inputs for energy-drought event statistics | Yes (Beta) | `energy.energy_drought_events()` returns the NumPy `RunStatistics` for a 1-D DataArray and an `xr.Dataset` for a gridded one: `count`, `total_duration`, `max_magnitude`, `valid_count`, and `frequency` over a `period` dimension plus the input's cell dimensions, and one `events` variable holding each cell's `RunSet`. A one-dimensional DataArray is evaluated eagerly, so a Dask-backed series is computed. Dask-backed gridded input must keep the time dimension in a single chunk, because events run the whole time axis; an empty time axis is rejected on the gridded path rather than reported as an empty summary. See `tests/test_energy_events.py`. |
| `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 energy degree days and wind profiles | Yes | `energy.heating_degree_days()`, `energy.cooling_degree_days()`, `energy.energy_degree_days()`, `energy.hub_height_wind_speed()`, and the auxiliary `energy.wind_shear_exponent()` preserve labels and match the NumPy functions exactly. Temperature units are inferred from CF attributes (K, °C, °F); bases remain explicit °C values. Wind units must be m s⁻¹ or absent. HDD/CDD attach `air_temperature_threshold`, EDD attaches both bases, and hub-height wind attaches `height`. Shared wind/alpha coordinates must match exactly; unlabeled non-scalar companions are rejected. Every dimension may have multiple Dask chunks. These stable DataArray calls do not emit `BetaFeatureWarning`. See `tests/test_energy_xarray.py`. These transforms describe climate-driven demand/resources, not actual load or generation. |
| `DataArray` inputs for compound hazards | Yes (Beta) | `compound.harmonize()` takes a mapping of DataArrays and returns a Dataset; `compound.joint_statistics()`, `compound.conditional_probability()`, and `compound.bootstrap_dependence()` take that Dataset. Results equal the NumPy routes exactly; `time` (and `hazard`) must be one Dask chunk, spatial chunks must match across inputs, and nothing is rechunked, loaded, or computed. Dask input leaves the value-dependent support warnings to `quality_flags`. Each call emits `BetaFeatureWarning`. See `tests/test_compound_xarray.py` and the {ref}`compound reference <compound-xarray-route>`. `compound.compound_events_xarray()` adds an eager numeric event catalog: it takes the prepared Dataset, rejects Dask-backed variables and coordinates without computing them, and returns one `event` row per matched pair beside a `cell` lookup table. See `tests/test_compound_catalog.py` and the {ref}`catalog reference <compound-event-catalog>`. |
| `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. `aggregation` reaches the fit and is recorded in each variable's calculation metadata. 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 | Fitting and stateful indices require one time chunk. PET adapters rechunk a split `time` internally. Elementwise Haines and energy transforms accept multiple chunks on every dimension; HDW requires only a single vertical core chunk. 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, 366, or 8784, 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, daily input on January 1, and hourly input at 00:00 on January 1. Daily values are converted to the 366-day calendar (February 29 synthesized from February 28 and March 1) and restored afterward; hourly values use the same calendar at 24 values per day ([#1177](https://github.com/monocongo/climate_indices/issues/1177)). A partial final year is supported; `cftime` calendars are rejected. See [ADR-0004](docs/adr/0004-xarray-calendar-semantics.md). |
| Automatic temporal inference | Yes, constrained | Adapter tests cover monthly and daily time-coordinate inference. Hourly inference applies only to APIs that accept hourly empirical standardization: `standardized_index()`, `energy.srepi()`, `energy.standardized_demand_index()`, and `energy.srli()`. SPI and SPEI reject hourly input. Sub-hourly and twelve-hourly cadences are rejected. |
| 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. |

## Energy transforms

The energy adapters use elementwise `xr.apply_ufunc` calls with no core dimension.
They do not fit a distribution or carry state, so ADR-0003's single-time-chunk
requirement does not apply. They neither import Dask nor choose a scheduler.
The `wind_shear_exponent()` output is an auxiliary diagnostic, not a registry
index; it has dimensionless units and output provenance.

Degree-day inputs represent daily mean temperature. An attached `time`
coordinate must contain consecutive Gregorian `datetime64` daily samples.
A Dask-backed `time` coordinate is rejected without computation; load the
coordinate before calling. The temperature data can remain lazy.
There is no January-start requirement or all-leap conversion, and the adapters
never synthesize February 29. Inputs without a time coordinate are accepted as
already-daily observations. Wind profiles have no calendar requirement.

Wind inputs must agree on shared geographic coordinates, including auxiliary
latitude/longitude coordinates, not only dimension indexes. Each shared
dimension requires coordinate labels on both inputs. Independent dimensions
still broadcast. To keep validation lazy, shared Dask-backed coordinates must
share one graph (the same array or an unmodified copy). Otherwise, load and
reconcile the coordinates before calling.
A labeled power-law `alpha` must have dimensionless units (`1`,
`dimensionless`, or no units attribute). Hub-height outputs replace scalar CF
height aliases rather than retaining the reference height as a contradictory
output coordinate. A varying coordinate named `height` or `hub_height`, or
tagged `standard_name="height"`, cannot describe a single target hub height
and is rejected.

The outputs are DataArrays, not Datasets. They preserve existing coordinates
and their attributes except for physical output coordinates. HDD/CDD replace
`air_temperature_threshold`; EDD replaces `heating_base_temperature` and
`cooling_base_temperature`. Each replacement sets the output value and attributes.
Hub-height wind replaces scalar height coordinates as described above.
The adapters do not invent time bounds from a timestamp whose interval
position is unknown. When writing CF datasets with degree-day time integrals,
retain the source daily time bounds in the enclosing Dataset. The scalar
threshold and height coordinates carry their own physical units.

## 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 plain-array workflows.
- Use the stable DataArray index APIs for labeled, gridded workflows with
  coordinate preservation and metadata.
- If you require stable interfaces, use the NumPy compound routes instead of
  the Beta compound Dataset routes and eager event catalogs.
- 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 xarray `DataArray` inputs. Both APIs are stable from 3.1.0. 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 stable from 3.1.0.

## 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.
