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

## docs/project-overview.md

# Project Overview

## Executive Summary

**climate_indices** is a production-grade Python library providing reference implementations of climate and hydrologic indices used for drought monitoring, agricultural planning, flood-potential assessment, and fire-weather assessment. The library computes the Standardized Precipitation Index (SPI), Standardized Precipitation Evapotranspiration Index (SPEI), Evaporative Demand Drought Index (EDDI), Percentage of Normal Precipitation (PNP), Precipitation Concentration Index (PCI), Potential Evapotranspiration (PET), the Palmer Drought Indices, a namespaced family of fire-weather indices, and a namespaced family of flood-potential indices (Effective Precipitation, Effective Drought Index (EDI), Flood Index, Antecedent Precipitation Index).

### Project Status
- **Maturity**: Production/Stable (Development Status 5)
- **Development Version**: 3.1.0 (not yet published)
- **Published Version**: 3.0.0
- **Python Support**: 3.10, 3.11, 3.12, 3.13, 3.14
- **License**: BSD 3-Clause
- **Documentation**: [ReadTheDocs](https://climate-indices.readthedocs.io/)
- **Repository**: [GitHub](https://github.com/monocongo/climate_indices)

### Key Capabilities
- **Drought / Moisture Indices**: SPI, SPEI, PNP, EDDI, PCI
- **Palmer Family**: PDSI, PHDI, PMDI, Palmer Z-Index, and self-calibrated scPDSI
- **Flood-Potential Indices**: `climate_indices.flood` namespace — Effective Precipitation, fixed-window EDI, Flood Index, and the recursive Antecedent Precipitation Index
- **PET Methods**: Thornthwaite, Hargreaves, and FAO-56 Penman-Monteith (`climate_indices.pet_penman_monteith`) are public APIs; lower-level primitives in `pm_eto.py` are internal and back the fire-weather indices
- **Fire-Weather Indices**: `climate_indices.fire` namespace — Fosberg Fire Weather Index (FFWI), Hot-Dry-Windy Index (HDW), Keetch-Byram Drought Index (KBDI), and the Canadian Forest Fire Weather Index System (FFMC, DMC, DC, ISI, BUI, FWI, DSR)
- **Run Theory**: `runs` module for identifying contiguous above/below-threshold events (duration, magnitude, intensity, peak, interarrival)
- **Flexible Input Formats**: NumPy arrays, xarray DataArrays, and Dask arrays
- **Multiple Temporal Scales**: Monthly data at 1–72 months; daily data at 1–2,196 days (where applicable)
- **Distribution Options**: Gamma and Pearson Type III distributions for SPI/SPEI
- **CLI Tools**: A single `climate_indices` CLI for batch processing NetCDF data across the index families
- **Scientific Rigor**: Peer-reviewed methodologies with comprehensive validation against reference datasets

## Project Classification

### Repository Structure
- **Type**: Monolith (single cohesive codebase)
- **Project Type**: Library
- **Primary Language**: Python 3.10+
- **Build System**: Hatchling (PEP 517) + uv for dependency management
- **Package Management**: uv (modern Python package installer/resolver)

### Architecture Pattern
**Layered Library Architecture**:
```
┌─────────────────────────────────────┐
│   CLI Layer                         │  ← __main__.py, _cli.py
├─────────────────────────────────────┤
│   Public API Layer                  │  ← typed_public_api.py, xarray_adapter.py,
│                                     │     cf_metadata_registry.py, fire/, flood/
├─────────────────────────────────────┤
│   Computation Layer                 │  ← indices.py, compute.py, palmer.py,
│                                     │     self_calibration.py, eto.py, pm_eto.py,
│                                     │     runs.py, validation.py
├─────────────────────────────────────┤
│   Math/Statistics Layer             │  ← lmoments.py, _palmer_duration.py,
│                                     │     _palmer_pdi.py, _palmer_wells.py
├─────────────────────────────────────┤
│   Infrastructure Layer              │  ← utils.py, logging_config.py,
│                                     │     exceptions.py, performance.py
└─────────────────────────────────────┘
```

### Entry Points
The installed script names map to the single full-featured CLI:
- **`climate_indices`** / **`process_climate_indices`**: `__main__.main`, the only CLI entry point. `--index` choices are `spi`, `spei`, `pnp`, `scaled`, `pet`, `palmers`, `kbdi`, `pe`, `edi`, `flood_index`, `api`, and `all`. The former specialized `spi` CLI (parameter save/load) was retired in 3.0.0; fit once with `compute.gamma_parameters()`/`compute.pearson_parameters()` and pass the dict to `fitting_params`.

### Technology Stack Summary
| Category | Technology |
|----------|-----------|
| **Core Language** | Python 3.10+ |
| **Scientific Computing** | NumPy, SciPy, xarray, Dask, pandas, cftime, h5netcdf |
| **Logging** | structlog (structured logging) |
| **Testing** | pytest, hypothesis (property-based), pytest-benchmark |
| **Type Checking** | mypy --strict |
| **Linting/Formatting** | ruff |
| **Build** | Hatchling (PEP 517) |
| **CI/CD** | GitHub Actions (`benchmarks.yml`, `release.yml`, `unit-tests-workflow.yml`) |
| **Containerization** | Docker (Python 3.14-slim base) |
| **Documentation** | Sphinx with ReadTheDocs hosting |

## Quick Reference

### Installation
```bash
# From PyPI
pip install climate_indices

# From source (development mode with uv)
git clone https://github.com/monocongo/climate_indices.git
cd climate_indices
uv sync --group dev
```

### Core API Usage

#### Modern xarray API (Recommended)
```python
import climate_indices as ci
from climate_indices import indices
import xarray as xr

# Load precipitation data
precip = xr.open_dataarray("precip.nc")

# Compute 6-month SPI (time parameters inferred from coordinates when omitted)
spi_6 = ci.spi(
    precip,
    scale=6,
    distribution=indices.Distribution.gamma,
    calibration_year_initial=1981,
    calibration_year_final=2010
)
```

#### Legacy NumPy API (Backward Compatibility)

> Note: the calibration-window boundary parameters are `calibration_year_initial`/
> `calibration_year_final` across the Python API. `percentage_of_normal` still
> accepts the older `calibration_start_year`/`calibration_end_year` keywords as
> deprecated aliases (removed in 4.0.0), and the CLI accepts
> `--calibration_year_initial`/`--calibration_year_final` as aliases for its
> `--calibration_start_year`/`--calibration_end_year` flags.
```python
from climate_indices import indices, compute
import numpy as np

# Load data as numpy array
precip_mm = np.load("precip_monthly.npy")

# Compute 6-month SPI
spi_6 = indices.spi(
    precip_mm,
    scale=6,
    distribution=indices.Distribution.gamma,
    data_start_year=1980,
    calibration_year_initial=1981,
    calibration_year_final=2010,
    periodicity=compute.Periodicity.monthly
)
```

### CLI Usage
```bash
# Process gridded NetCDF data for multiple indices
climate_indices \
    --index spi \
    --periodicity monthly \
    --scales 1 3 6 12 \
    --netcdf_precip precip.nc \
    --var_name_precip prcp \
    --calibration_start_year 1981 \
    --calibration_end_year 2010 \
    --output_file_base results/spi
```

### Development Commands
```bash
# Run tests
uv run pytest

# Run tests with coverage
uv run pytest --cov=src --cov-report=term

# Type checking (Python 3.12, the stack CI's lint job uses)
uv run --python 3.12 mypy src/ tests/test_type_checking.py

# Linting and formatting
uv run ruff check src/ tests/
uv run ruff format --check src/ tests/

# Run benchmarks (deselected by default)
uv run pytest -m benchmark
```

## Project Goals and Use Cases

### Primary Use Cases
1. **Drought Monitoring**: Real-time and historical drought assessment using SPI/SPEI/PDSI
2. **Agricultural Planning**: Crop yield forecasting and irrigation scheduling
3. **Water Resource Management**: Reservoir operations and water allocation decisions
4. **Flood Assessment**: Daily flood-potential screening with Effective Precipitation, EDI, Flood Index, and API
5. **Fire-Weather Assessment**: Fuel-dryness and fire-danger indices for wildland fire management
6. **Event Analysis**: Run-theory extraction of drought, wet, and fire-spell events
7. **Climate Research**: Long-term precipitation and drought trend analysis
8. **Operational Meteorology**: Integration into weather services and early warning systems

### Design Philosophy
- **Scientific Correctness**: Implementations follow peer-reviewed methodologies (McKee et al. 1993, Vicente-Serrano et al. 2010, Thornthwaite 1948, Hargreaves & Samani 1985, Palmer 1965, Wells et al. 2004, Allen et al. 1998 FAO-56)
- **Performance**: Optimized for large gridded datasets using Dask parallelization
- **Usability**: Dual API (NumPy/xarray) for different user workflows
- **Reliability**: Comprehensive test coverage (>90%), property-based testing, backward compatibility guarantees
- **Extensibility**: Modular design allows addition of new indices, distributions, or PET methods

## AI-Assisted Development Guidance

### Recommended Starting Points for AI Agents
1. **For understanding computation**: Start with `docs/architecture.md`, then `src/climate_indices/compute.py` and `src/climate_indices/indices.py`
2. **For understanding the public API**: Read `src/climate_indices/typed_public_api.py` (strict mypy typing) and `src/climate_indices/xarray_adapter.py`
3. **For understanding the CLI**: Examine `src/climate_indices/__main__.py` (the index pipeline registry and `--index` handling) and `src/climate_indices/_cli.py` (shared argument registration)
4. **For understanding Palmer calculations**: Read `src/climate_indices/palmer.py`, `src/climate_indices/_palmer_pdi.py`, `src/climate_indices/_palmer_wells.py`, `src/climate_indices/_palmer_duration.py`, and `src/climate_indices/self_calibration.py`
5. **For understanding fire-weather indices**: Read `src/climate_indices/fire/__init__.py` and the subsystem design in `docs/design/fire-subsystem.md`
6. **For understanding flood-potential indices**: Read `src/climate_indices/flood/__init__.py` (public dispatch) and the `flood/_pe.py`, `flood/_edi.py`, `flood/_if.py`, `flood/_antecedent.py` kernels
7. **For understanding run theory**: Read `src/climate_indices/runs.py`
8. **For testing patterns**: Review `tests/conftest.py` (fixtures), `tests/test_xarray_adapter.py` (modern API), `tests/test_property_based.py` (invariants)
9. **For error handling**: Study `src/climate_indices/exceptions.py` (complete hierarchy with attributes)

### Critical Architectural Invariants
- **Time dimension chunking**: SPI/SPEI and stateful fire adapters reject split Dask time chunks; PET adapters rechunk split time internally. See [chunking guidance](docs/xarray_compatibility.md#chunking-guidance-for-gridded-indices).
- **Calibration period**: Default minimum 30 years; violations trigger `ShortCalibrationWarning`
- **Pearson Type III fitting**: Requires at least 4 non-zero values per calendar time step; fewer raise `InsufficientDataError`
- **Coordinate validation**: Time coordinates must be monotonically increasing; xarray inputs undergo automatic validation
- **Backward compatibility**: Legacy NumPy API (`indices.py`) must remain stable; new features go to the xarray API
- **Fire namespace**: Fire-weather indices are exposed only under `climate_indices.fire`; they are not top-level package functions (see ADR-0005)
- **Flood namespace**: Flood-potential indices are exposed only under `climate_indices.flood` (see ADR-0013); `edi` (flood) and `eddi` (drought) are distinct indices

### Code Patterns to Follow
1. **Type annotations**: All public functions require full type hints (package-wide strict mypy checking configured in `pyproject.toml`)
2. **Docstrings**: Google-style docstrings with parameter descriptions, return types, raises, examples
3. **Error handling**: Use specific exception types from `exceptions.py` with context attributes
4. **Logging**: Use structlog with structured key-value pairs; avoid string interpolation in log messages
5. **Testing**: Property-based tests for mathematical invariants, regression tests against reference datasets, benchmark tests for performance

### Dependencies and Constraints
- **Core dependencies**: cftime>=1.6.4.post1, dask>=2025.7.0, h5netcdf>=1.6.3, numpy>=1.24.0, pandas>=2.0.0, scipy>=1.15.3, structlog>=24.1.0, xarray>=2025.6.1
- **Python version**: Must support Python 3.10-3.14 (`requires-python = ">=3.10,<3.15"`)
- **Line length**: 120 characters (ruff configuration)
- **Import order**: stdlib → third-party → local (enforced by ruff)
- **Test markers**: Use `@pytest.mark.benchmark` for performance tests, `@pytest.mark.slow` for long-running tests

### Common Pitfalls
1. **Do NOT** use wildcard imports (`from module import *`) - explicitly forbidden by ruff
2. **Dask time chunks**: SPI/SPEI and stateful fire adapters reject split time chunks; PET adapters accept them and rechunk internally
3. **Do NOT** modify the `indices.py` API in a breaking way - backward compatibility requirement
4. **Do NOT** commit without running tests - pre-commit hooks will fail
5. **Do NOT** use string paths - always use `pathlib.Path` objects

### Key Files for Modification Scenarios
| Task | Primary Files | Test Files |
|------|--------------|------------|
| Add new index | `compute.py`, `indices.py`, `typed_public_api.py`, `xarray_adapter.py` | `test_compute.py`, `test_indices.py`, `test_typed_public_api.py`, `test_xarray_adapter.py` |
| Add new distribution | `compute.py`, `indices.py`, `lmoments.py` | `test_compute.py`, `test_property_based.py` |
| Modify Palmer indices | `palmer.py`, `_palmer_pdi.py`, `_palmer_wells.py`, `_palmer_duration.py`, `self_calibration.py` | `test_palmer.py`, `test_scpdsi.py`, `test_self_calibration.py` |
| Add a fire-weather index | `fire/` subpackage | `test_fire*.py`, reference tests |
| Add a flood-potential index | `flood/` subpackage, `typed_public_api.py` | `test_flood*.py`, `test_main_flood.py` |
| Modify run extraction | `runs.py` | `test_runs.py` |
| Fix CLI bug | `__main__.py`, `_cli.py` | `test_cli_common.py`, `test_main_*.py` |
| Add validation | `validation.py`, `xarray_adapter.py`, `exceptions.py` | `test_input_validation.py`, `test_exceptions.py` |
| Update CF metadata | `cf_metadata_registry.py` | `test_cf_metadata.py`, `test_metadata_validation.py` |
| Performance optimization | `compute.py`, `performance.py`, chunking strategies | `test_benchmark_*.py` |

---

**Next Steps**: See [architecture.md](docs/architecture.md) for detailed technical architecture, [development-guide.md](docs/development-guide.md) for setup instructions, and [reference.md](docs/reference.md) for API reference details.

## docs/architecture.md

# Technical Architecture

## Executive Summary

The **climate_indices** library implements a **layered library architecture** optimized for scientific computing on climate data. The architecture separates concerns across five distinct layers: CLI interfaces for batch processing, dual public APIs (legacy numpy + modern xarray), core computation logic, mathematical/statistical algorithms, and infrastructure services.

### Architectural Highlights
- **Dual API Design**: Maintains backward-compatible numpy API while providing modern xarray/Dask integration
- **Layered Separation**: Clean boundaries between user interfaces, computation, and infrastructure
- **Parallelization Strategy**: Multiprocessing for CLI, Dask for xarray workflows
- **Exception Hierarchy**: Structured error handling with context-rich exceptions
- **Test Architecture**: Comprehensive fixture-based testing with property-based validation

### Design Principles
1. **Scientific Correctness**: Implementations strictly follow peer-reviewed methodologies
2. **Backward Compatibility**: Legacy numpy API remains stable across minor versions
3. **Performance**: Optimized for gridded datasets with 10⁶+ cells
4. **Type Safety**: Package-wide strict mypy checking configured in `pyproject.toml`
5. **Observability**: Structured logging with performance metrics

## Technology Stack

### Core Dependencies
| Dependency | Version | Purpose |
|------------|---------|---------|
| **Python** | 3.10-3.14 | Language runtime |
| **scipy** | >=1.15.3 | Statistical distributions, numerical optimization |
| **xarray** | >=2025.6.1 | Labeled multi-dimensional arrays, CF metadata |
| **dask** | >=2025.7.0 | Parallel computation, lazy evaluation |
| **h5netcdf** | >=1.6.3 | NetCDF file I/O backend |
| **cftime** | >=1.6.4 | Calendar-aware datetime handling |
| **structlog** | >=24.1.0 | Structured logging with context |

### Development Dependencies
| Tool | Purpose |
|------|---------|
| **pytest** | Test runner with fixtures and parametrization |
| **hypothesis** | Property-based testing for invariants |
| **pytest-benchmark** | Performance regression testing |
| **pytest-cov** | Code coverage reporting |
| **mypy** | Static type checking |
| **ruff** | Linting and code formatting |
| **sphinx** | Documentation generation |

### Build and Deployment
- **Build Backend**: Hatchling (PEP 517 compliant); the optional Rust extension is built with maturin (see [Optional Rust Backend](#optional-rust-backend))
- **Package Manager**: uv (modern resolver, lockfile support)
- **CI/CD**: GitHub Actions (3 workflows: unit tests, releases, benchmarks)
- **Container**: Docker with Python 3.14-slim base image
- **Documentation**: Sphinx with ReadTheDocs hosting

## Layered Architecture Pattern

```
┌────────────────────────────────────────────────────────────────────┐
│                         CLI Layer                                  │
├────────────────────────────────────────────────────────────────────┤
│  __main__.py          │  Full-featured CLI (all indices)          │
├────────────────────────────────────────────────────────────────────┤
│                      Public API Layer                              │
├────────────────────────────────────────────────────────────────────┤
│  typed_public_api.py  │  Strict mypy, xarray wrappers            │
│  xarray_adapter.py    │  CF-compliant xarray interface           │
│  validation.py        │  Shared input validation facade          │
│  indices.py           │  Legacy numpy API (backward compat)       │
├────────────────────────────────────────────────────────────────────┤
│                    Computation Layer                               │
├────────────────────────────────────────────────────────────────────┤
│  compute.py           │  Core algorithms (scaling, fitting, PDF) │
│  palmer.py            │  Palmer drought indices implementation    │
├────────────────────────────────────────────────────────────────────┤
│                 Math/Statistics Layer                              │
├────────────────────────────────────────────────────────────────────┤
│  eto.py               │  PET: Thornthwaite, Hargreaves methods   │
│  lmoments.py          │  L-moments for distribution fitting       │
├────────────────────────────────────────────────────────────────────┤
│                  Infrastructure Layer                              │
├────────────────────────────────────────────────────────────────────┤
│  utils.py             │  Utilities, calendar conversions         │
│  logging_config.py    │  Structured logging setup                │
│  exceptions.py        │  Exception hierarchy with context        │
│  performance.py       │  Performance metrics tracking            │
└────────────────────────────────────────────────────────────────────┘
```

### Layer Responsibilities

#### 1. CLI Layer
**Purpose**: Command-line interfaces for batch processing NetCDF datasets.

**Modules**:
- **`__main__.py`**: Full-featured CLI supporting SPI, SPEI, PET, Palmer, and PNP indices
  - Multiprocessing pool for gridded data parallelization
  - NetCDF dimension validation through the validation facade, and coordinate
    conversion
  - Shared memory arrays for worker processes
  - Dataset layout detection (grid, divisions, timeseries)

**Entry Points** (`pyproject.toml`):
```toml
[project.scripts]
climate_indices = "climate_indices.__main__:main"
process_climate_indices = "climate_indices.__main__:main"
```

#### 2. Public API Layer
**Purpose**: User-facing interfaces for programmatic access.

**Modules**:
- **`typed_public_api.py`**: Modern xarray API with strict mypy compliance
  - Type-safe wrappers for SPI and SPEI
  - Enforces keyword-only arguments
  - Full mypy --strict compliance

- **`xarray_adapter.py`**: CF-compliant xarray interface
  - Coordinate validation and alignment
  - CF metadata preservation and generation
  - Dask array support with chunking validation
  - PET computation (Thornthwaite, Hargreaves)

- **`validation.py`**: Public validation facade shared by the CLI, the xarray
  and fire adapters
  - Input type detection (`DataArray`, `ndarray`; `Dataset` is rejected with a
    select-a-variable hint)
  - Dataset-layout classification (grid, divisions, timeseries) and the
    dimension orders each layout accepts
  - Time-dimension, monotonicity, and Dask single-chunk validation

- **`indices.py`**: Legacy numpy API
  - Backward-compatible function signatures
  - Direct numpy array inputs/outputs
  - Distribution enum (`Distribution.gamma`, `Distribution.pearson`)
  - SPI, SPEI, PNP, PET computation

**API Design Decision**: The dual API approach allows:
- **Legacy users**: Continue using numpy arrays without migration
- **Modern users**: Leverage xarray's labeled dimensions, CF metadata, and Dask parallelization
- **Migration path**: xarray API uses numpy implementation internally

#### 3. Computation Layer
**Purpose**: Core mathematical algorithms for climate index calculation.

**Modules**:
- **`compute.py`**: Core computation functions
  - `prepare_scaled()`: Shared flatten/clip/roll-sum/reshape preparation for the fitting-based indices
  - `is_all_missing()`, `reshape_time_major()`, `unfold_time_major()`: single owner of all-missing detection and the time-major block round trip
  - `fit_and_standardize()`: Shared parameter normalization, gamma/Pearson dispatch, and Pearson→gamma fall-back seam for SPI and SPEI
  - `gamma_parameters()`, `pearson_parameters()`: Distribution fitting
  - `transform_fitted_gamma()`, `transform_fitted_pearson()`: CDF transformation
  - `sum_to_scale()`: Optimized sliding window summation
  - `Periodicity` enum: `monthly` (12 steps/year), `daily` (366 steps/year), with a `period_length` property and a `unit()` method

- **`palmer.py`**: Palmer Drought Index family
  - PDSI (Palmer Drought Severity Index)
  - PHDI (Palmer Hydrological Drought Index)
  - PMDI (Palmer Modified Drought Index)
  - ZINDEX (Palmer Z-Index)
  - scPDSI (Self-calibrated Palmer Drought Severity Index), available through the NumPy `palmer.scpdsi()` API

**Key Algorithms**:
1. **SPI/SPEI Computation**:
   ```
   Input: precip (or P-PET)
   ↓
   Scale to an N-step window (sum_to_scale)
   ↓
   Fit distribution (gamma or Pearson Type III) per calendar month/day
   ↓
   Transform to standard normal via CDF (scipy.stats)
   ↓
   Output: SPI/SPEI values
   ```

2. **Distribution Fitting Strategy**:
   - **Gamma**: Method of moments (α, β parameters)
   - **Pearson Type III**: L-moments (location, scale, skew parameters)
   - **Calibration period**: Default 30+ years, user-configurable
   - **Handling zeros**: Probability of zero tracked separately
   - **Fall-back policy**: SPI falls back from a failed Pearson Type III fit, or one that
     loses most of the input's valid values, to gamma; SPEI propagates the failure.
     `compute.fit_and_standardize()` takes this as `fallback_to_gamma`, so the divergence
     is a parameter of one seam rather than a copy of the fit/dispatch branch in each
     index function.

#### 4. Math/Statistics Layer
**Purpose**: Low-level mathematical and statistical functions.

**Modules**:
- **`eto.py`**: Potential Evapotranspiration methods
  - **Thornthwaite (1948)**: Monthly PET from temperature and latitude
    - Heat index computation
    - Day length adjustment based on latitude
  - **Hargreaves (1985)**: Daily PET from temperature range and solar radiation
    - Requires tmin, tmax, and latitude
    - Extraterrestrial radiation calculation

- **`lmoments.py`**: L-moments for robust distribution fitting
  - Implements Hosking (1990) L-moments algorithm
  - Used for Pearson Type III parameter estimation
  - More robust than method of moments for skewed distributions

**Design Note**: This layer has no dependencies on upper layers and could be extracted as standalone utilities.

#### 5. Infrastructure Layer
**Purpose**: Cross-cutting concerns (utilities, logging, error handling).

**Modules**:
- **`exceptions.py`**: Exception hierarchy
  - Base: `ClimateIndicesError` (catch-all for library errors)
  - Computation: `DistributionFittingError`, `InsufficientDataError`, `PearsonFittingError`
  - Validation: `DimensionMismatchError`, `CoordinateValidationError`, `InputTypeError`, `InvalidArgumentError`
  - Warnings: `MissingDataWarning`, `ShortCalibrationWarning`, `GoodnessOfFitWarning`, `InputAlignmentWarning`
  - All exceptions carry context attributes (e.g., `distribution_name`, `input_shape`, `parameters`)

- **`logging_config.py`**: Structured logging configuration
  - `configure_logging()`: Sets up structlog with JSON serialization
  - Console: Human-readable colored output
  - File: JSON-formatted for log aggregators
  - Context binding for tracing

- **`utils.py`**: Utility functions
  - Calendar conversions: `transform_to_366day()`, `transform_to_gregorian()`
  - Calendar planning: `DailyCalendarPlan`
  - Data validation: `is_data_valid()`
  - Array reshaping: `reshape_to_2d()`, `reshape_to_divs_years_months()`
  - Periodicity utilities: `gregorian_length_as_366day()`

- **`performance.py`**: Memory metrics
  - `get_process_memory_mb()`: current process memory usage
  - `check_large_array_memory()`: returns memory metrics when arrays exceed the 1 GB threshold

## Optional Rust Backend

Numerical kernels are being migrated to Rust as an internal acceleration
backend. The public Python API is the compatibility and citation contract:
users import and call `climate_indices` exactly as before, and nothing in the
public signatures or results depends on whether the backend is installed.

```
Cargo.toml                       # Cargo workspace
crates/climate-indices/          # pure Rust numerical kernels, published as `climate-indices`
crates/climate-py/               # PyO3/NumPy bindings -> climate_indices._native
src/climate_indices/_native.pyi  # type stub for the extension
```

- **`climate-indices`** is pure numerical Rust. It has no PyO3, NumPy bindings,
  Python exceptions, or CPython assumptions, receives already-validated
  numbers, and returns numbers. It publishes to crates.io on its own `0.x`
  track, independent of the Python distribution's version
  ([ADR-0021](docs/adr/0021-climate-indices-crate-identity.md)). The reserved
  `0.0.0` release contains no algorithms; the first real release supersedes it.
- **`climate-py`** is the only crate that knows about Python. It converts
  arrays and errors at the boundary and holds no climate algorithm. It stays
  `publish = false`, inherited from the workspace, and reaches a user only as
  the wheel `climate_indices._native`.
- **Dispatch** lives in the Python module that owns each computation. It
  imports `climate_indices._native` inside `try`/`except ImportError`; without
  the extension every computation runs in pure Python. A runtime error raised
  by the extension propagates and is never silently retried in Python.
- **Reference implementation**: the Python implementations stay alive and
  directly testable. They are the oracle for Rust parity tests until a separate,
  explicit decision retires them.

**Ported kernels.** The gamma fit and transform behind SPI were the first port,
followed by EDDI's empirical ranking and inverse normal, the distribution fits
SPEI adds, Pearson Type III (also used by SPI and the standardized index), the
log-logistic (generalized logistic, GLO), the PNP and PCI numerical blocks, the
fire-weather recurrences, the flood family, and the Palmer family's water
balance, Z-index, and spell recursions. The gamma and distribution-fit kernels
replace the
numerical blocks inside `compute.py` functions; the EDDI, PNP, and PCI kernels
replace the blocks inside the `indices.py` functions that own them, so the
validation, calibration-period resolution, data-quality and goodness-of-fit
warnings, the Pearson-to-gamma fallback, logging, zero placement, support-limit
masks, and output scaling around them still run in Python, and every caller of
those functions (SPI, SPEI and the standardized index, `fit_diagnostics`, the
xarray adapter) uses them:

| Python seam | Rust kernel (`climate-indices`) |
|---|---|
| method-of-moments block of `compute.gamma_parameters` | `gamma::gamma_parameters` |
| `scipy.stats.gamma.cdf` and zero-mass mixing in `compute.transform_fitted_gamma` | `gamma::gamma_probabilities` |
| `scipy.stats.norm.ppf` in `compute.transform_fitted_gamma` | `special::norm_ppf` |
| the year loop of `eto.eto_thornthwaite` | `eto::thornthwaite` |
| the day loop of `eto.eto_hargreaves` | `eto::hargreaves` |
| FAO-56 Eq 6 in `pm_eto.pm_eto` | `pm_eto::pm_eto` |
| the intermediate chain of `pm_eto.penman_monteith_eto` | `pm_eto::penman_monteith_eto` |
| sample L-moments and the Pearson Type III fit of `compute.pearson_parameters` (`lmoments.fit`, `fit_spatial`) | `lmoments::sample_lmoments`, `pearson::pearson_parameters` |
| `scipy.stats.pearson3.cdf` in `compute._pearson_fit` | `pearson::pearson_cdf_block` |
| GLO fit of `compute.loglogistic_parameters` (`lmoments.fit_glo`, `fit_glo_spatial`) | `loglogistic::loglogistic_parameters` |
| GLO probability step of `compute._loglogistic_fit` | `loglogistic::loglogistic_cdf_block` |
| rank count and Tukey plotting position in the per-period loop of `indices.eddi` | `eddi::tukey_probabilities` |
| `indices._hastings_inverse_normal` | `eddi::hastings_inverse_normal` |
| calibration normals and ratios in `indices.percentage_of_normal` | `pnp::pnp_normals`, `pnp::pnp_percentages` |
| monthly reduction and ratio in `indices.pci` | `pci::pci` |

The PET entry points follow the same dispatch rule as the gamma kernels. The
public FAO-56 helper functions stay Python callables; `pm_eto.penman_monteith_eto`
resolves the humidity and radiation pathways, and checks the wind measurement
height and the `rh_min`-without-`rh_max` case, before the kernel is reached, so
the Python and Rust paths raise the same error in the same order. `eto` also
keeps a Thornthwaite block with an all-NaN month column in Python, since only the
Python path reports `np.nanmean`'s empty-slice warning. Each kernel copies every
operand it reads before it releases the GIL, so the native route holds the caller's
arrays, flattened real-array inputs, and those copies: a bounded multiple of the
request. Penman-Monteith passes constants as scalars, caches station/astronomy
terms on exact input bits, and evaluates Eq 6 without intermediate arrays.
Thornthwaite/Hargreaves still pass constant latitudes as zero-stride views, and
Hargreaves reports the bytes it copies beside the arrays it is handed, so the
logged memory model covers the route the dispatch selected. Their measured effect on
three representative inputs is in `benchmarks/README.md`; RUST-011 owns whether
each kernel is worth its conversion overhead.

EDDI is non-parametric, so no SciPy special function is ported for it: its
probabilities are a count of the period's climatology values strictly below each
value, converted by the Tukey plotting position. The chunking Python applies
across cells (`_EDDI_RANK_COMPARISON_ELEMENT_BUDGET`) only bounds an intermediate
under the Python path; the Rust kernel walks every column in one pass, which
cannot change a count. The calibration-period resolution, the leading-scale-pad
mask, and the unfolding to the caller's layout stay in Python, as they did.

The special functions are line-by-line ports of the Cephes `igam`, `ndtri`,
`ndtr`, and `lgam` that SciPy 1.17 evaluates, since a generic implementation
would not hold the parity contract in the transformed tails. The L-moment fits
route to Rust only for a block whose values are all NaN or at most `1e100` in
magnitude: an infinity (or an overflowing weighted sum) makes the L-moments NaN,
which the single-series Python fit returns as NaN parameters but the cell-axis
fit marks invalid, so those blocks keep the Python fit. The `climate_indices.lmoments`
log records that the single-series fits write per failed step (an `ERROR` for
invalid L-moments, and a `WARNING` for a step with fewer than four non-NaN
values) are not written by the Rust fit, as they are not by the cell-axis fit;
the failed-fit count and the high-failure-rate warning are unchanged, and a
failure rate at or below that warning's threshold logs one
`distribution_fitting_failures` `WARNING` with the counts on every backend. SciPy's
last-bit results depend on whether its build fuses multiply-adds (aarch64 builds
do, x86-64 wheels do not), and the Pearson fit's `exp(gammaln(a) - gammaln(a + 0.5))`
amplifies one ulp by up to `1e11` for a near-symmetric sample, so the ported
`lgam` and polynomial helpers fuse on aarch64 and only there (`special::mul_add`).
The fire-weather recurrences are the one port that does not replace a block
inside a Python function: each replaces a recurrence's whole day loop, so the
Rust side owns the time axis while Python keeps the calendar, the validation,
and the error surface around it:

| Python seam | Rust kernel (`climate-indices`) |
|---|---|
| `fire.ffmc` recurrence | `fire::ffmc` over `recurrence::run` |
| `fire.duff_moisture_code` recurrence | `fire::dmc` over `recurrence::run` |
| `fire.drought_code` recurrence | `fire::dc` over `recurrence::run` |
| `fire.kbdi` recurrence | `fire::kbdi` over `recurrence::run` |

`recurrence::run` ports the shared day loop in `climate_indices._recurrence`:
the ADR-0007 missing-day policy, the ADR-0010 seasonal carry mask, the spin-up
offset, and the recorded-history NaNs are identical in both.
`climate_indices.fire._native` builds each code's kernel from the same
`_CodeInputs` its Python step reads, so the two paths cannot disagree about
which arrays a code consumes, and the KBDI, FFMC, DMC, and DC daily updates are
line-by-line ports of the Python expressions, including the
`np.maximum`/`np.minimum` NaN propagation and NumPy's operation order. The
combined `cffwis()` orchestrator runs each of its three codes through the same
kernels; the components are independent, so running them outside the shared day
loop cannot change a result. Elementwise fire indices (ISI, BUI, FWI, DSR,
Fosberg, HDW, Haines) stay in Python: they are single NumPy expressions with no
recurrence, and a port would not pay for itself.

The Palmer family ports its loops and state machines. Dispatch for every stage
lives in `palmer.py`, which owns `pdsi()` and `scpdsi()`; the private
`_palmer_pdi`, `_palmer_wells`, and `self_calibration` functions stay pure
Python and remain the oracle when called directly:

| Python seam | Rust kernel (`climate-indices`) |
|---|---|
| `palmer._calc_water_balances` (with `_calc_potential_loss`, `_calc_recharge`) | `palmer::water_balance` |
| `palmer._calc_k_prime_and_dbar` (PDSI K factors and scPDSI K-prime) | `palmer_zindex::k_prime_and_dbar` |
| `palmer._calc_raw_zindex` (with `_calc_cafec_zindex`) | `palmer_zindex::raw_zindex` |
| `_palmer_pdi.calculate`, from `_calculate_pdsi_prepared` | `palmer_pdi::calculate` |
| `self_calibration.duration_factors`, from `_calculate_scpdsi_prepared` | `self_calibration::duration_factors` |
| `_palmer_wells.calculate`, from `_calculate_scpdsi_prepared` | `palmer_wells::calculate` |

The CAFEC ratios (`_calc_cafec_coefficients`), the T ratio
(`_calc_zindex_factors`), the K-factor normalization in `_calc_kfactors`, and the
scPDSI percentile rescaling (`nan_safe_percentile`, `_rescale_scpdsi_zindex`) stay
in Python for the fire indices' reason: each is one NumPy expression. So do the
`DurationFactors`/`PdiDurationFactors` validation, the scPDSI K-prime finiteness
check, masks, fully-missing-cell NaNs, and logging; the Wells kernel takes the
recurrence coefficients Python derived. The PDI state machine Python vectorizes
across cells partitions every cell into exactly one branch each month, so the
kernel runs it one cell at a time with the same branch order and exact-zero
comparisons. A scalar AWC goes to Rust only when it is a Python `int`/`float` or
an `np.float64` (other scalars keep their NumPy type promotion in Python), and an
array AWC only when it is one plain float64 value per cell. An infinite Z value
and a calibration Z series too short for the longest rolling window keep the
Python path, which raises its own error; the kernels' abatement and least-squares
failures raise `_native.NoConvergenceError`, which `palmer.py` re-raises as the
`ConvergenceError` the Python path raises. The native water balance releases its
unused Python output placeholders before allocating Rust outputs; its inputs are
still copied before releasing the GIL. The K-prime and raw Z-index stages borrow
the CAFEC arrays with the GIL held, avoiding five full-record copies per stage.

The flood family ports the computation behind each NumPy entry point in
`climate_indices.flood`, with the Antecedent Precipitation Index running on the
same `recurrence::run` driver as the fire codes rather than a second recurrence
interface:

| Python seam | Rust kernel (`climate-indices`) |
|---|---|
| `scipy.ndimage.correlate1d` window of `flood.effective_precipitation` | `flood::effective_precipitation` |
| calendar-day standardization of `flood.edi` | `flood::edi` |
| annual maxima and standardization of `flood.flood_index` | `flood::flood_index` |
| `flood.antecedent_precipitation_index` recurrence | `flood::antecedent_precipitation_index` over `recurrence::run` |

The effective-precipitation kernel reproduces `correlate1d`'s general loop
(the newest day's term first, then the window from its oldest day) and its
NaN propagation, so a window holding a missing day is NaN. The EDI and Flood
Index calibration sums follow NumPy's axis-0 reduction, which adds the years
sequentially for a block of several columns and pairwise for a single one (a
1-D Flood Index sample); both share the population SD and the `8 * eps * |mean|`
rounding guard. `climate_indices.flood._native` hands each kernel the blocks and
the Calibration Period rows the Python modules resolved; validation, the
all-leap layout, and the xarray and CLI paths stay in Python, and the xarray
adapters reach the kernels through the same NumPy entry points. An API decay
constant whose type would promote the Python step beyond float64 (an extended
`np.longdouble`) keeps the Python recurrence.

Dispatch takes the Rust path only for a
plain, aligned float64 `ndarray` whose fit parameters are aligned and one per
calendar step (and cell). Unaligned arrays, other dtypes, and caller-supplied
parameters that vary by year run the Python implementation. A mask is a
missing-value marker, so the seams that read it as one replace it with NaN
before this check, and the prepared plain float64 array may then use Rust: the
gamma transform (`transform_fitted_gamma` and the parameter resolver it calls),
the GLO transform (`transform_fitted_loglogistic`, whose fit and CDF may then run
natively), the sliding sum that prepares a scaled series, and PNP preparation. A
partially masked SPI/SPEI gamma transform therefore matches the same input
passed through `np.ma.filled(values, np.nan)`. Fully masked SPI/SPEI inputs
short-circuit before fitting or native calls, returning a `MaskedArray`, not a
filled plain NaN array. Existing preparation still applies: SPI flattens a 2-D
series, and SPEI forms its water-balance array, so input object identity is not
a general guarantee. A direct `transform_fitted_gamma` call returns the fully
masked input unchanged. The gamma parameter resolver fills even a
full mask with NaN, so a fully masked fit returns plain NaN parameters without
calling a kernel; supplied `alpha`/`beta` come back unchanged, and only the
computed `prob_zero` is NaN. A direct `gamma_parameters` call also returns plain
NaN parameters for a fully masked input, but retains Python's mask semantics for
a partial mask. Seams that do
not replace a mask keep their Python implementation when handed one: a direct
`gamma_parameters` or `loglogistic_parameters` call, the Pearson Type III fit
and transform, PCI, whose dispatch checks the original 1-D input and requires a
plain 1-D array, and the fire recurrences, which apply the same guard to every
weather array and to the seed they resume from, as the flood kernels do to their
prepared series and the API seed. A layout the kernel cannot take unchanged (an
empty axis, a non-float64 or unaligned array, or, for the fire recurrences, a
time-first array whose spatial axes cannot be viewed as one cell axis without a
copy) stays in Python as well. The kernels take views of the prepared arrays, so
a broadcast month series or season mask is copied once, by the binding, rather
than first materialized in Python; that boundary copy is the native path's own
buffer and is not part of the `array_memory_mb` a recurrence reports, which
counts the one recorded history the kernel returns. The flood kernels differ in
one respect: a time-last (transposed) array reshapes to a view, and a layout
that cannot merge without a copy, such as a regional slice of a larger grid or a
Fortran-ordered array, is copied once in Python and still runs natively, because
the Python path is roughly
three times slower than the kernel plus that copy (365-day PE over 400 cells).
The public flood entry points already hand the dispatch a C-contiguous array,
since their validation copies, so the policy matters to direct callers of
`flood._native`. A loaded extension that predates a kernel,
or a recurrence option wider than the binding's integer parameters, also keeps
the recurrence on its Python steps rather than failing at the boundary.
Native dispatch also requires NumPy floating-point errors to be ignored
(`np.errstate(all="ignore")`); policies that warn, raise, call, log, or print keep
the Python path. A warning filter that promotes `RuntimeWarning` (or a superclass)
to an exception also keeps Python, even under `all="ignore"`. Python 3.14
context-aware warnings conservatively keep the Python path. Default NumPy error
policies therefore use Python even when the extension is installed.

For gamma calibration blocks, the contract is:

- For a plain (or NaN-normalized) calibration block, if any column has no
  positive value after zero replacement (only NaN, zero, or negative values),
  the **whole fit** stays in Python. Its `np.nanmean` emits
  `RuntimeWarning: Mean of empty slice` independently of `np.errstate`, so that
  warning and its promotion to an exception are preserved. A block passed as a
  `MaskedArray` is not normalized here: its masked reductions warn nothing and
  return masked parameters. Entirely missing inputs still return NaN parameters
  before fitting, without this warning.
- Negatives are not removed or newly rejected. A block containing negatives may
  use the native fit only when every column also has a positive value and all
  floating-point errors are ignored. Under `invalid="warn"` or `invalid="raise"`,
  Python preserves `RuntimeWarning: invalid value encountered in log` or
  `FloatingPointError`, respectively.
- Constant and single-positive-value columns may use Rust under `all="ignore"`;
  their degenerate parameters and their transform results match Python.
  A Python fit does not prohibit later eligible native CDF/inverse-normal calls.
  Return values and `climate_indices` warnings match on both backends.

These rules retain existing behavior; [ADR-0017](docs/adr/0017-rust-core-acceleration-backend.md)
records the decision. Direct extension
calls reject unaligned inputs and copy empty arrays without creating Rust views
of caller-owned storage. `tests/test_native_parity.py` (gamma),
`tests/test_native_parity_distributions.py` (Pearson Type III and GLO),
`tests/test_native_parity_fire.py` (the fire recurrences, which compare the
returned state as well), `tests/test_native_parity_palmer.py` (the Palmer
family, which also requires the backtracking's sign pattern to match exactly),
and `tests/test_native_parity_flood.py` (the flood family, including the
returned API state and a resumed run that is bitwise a single pass) explicitly
ignore floating-point errors and compare the
two paths at `rtol = atol = 1e-10` with matching NaN positions, and the
`python_backend` fixture in
`tests/conftest.py` pins any test to the Python reference. The
consolidated suite in `tests/test_native_parity_registry.py` drives the same
comparison from one registry
(`tests/parity_registry.py`: Python entry point, dispatch module, expected
kernels, input family, tolerance), adds Hypothesis draws over lengths, NaN
patterns, zero runs, extreme magnitudes, and spatial shapes, after each entry's
input derivation. Palmer entries replace gaps and clip precipitation to [5, 300],
so their property cases cover lengths and bounded precipitation variation, not
gaps, zero runs, or extreme magnitudes. The suite also asserts the documented
dispatch decision for the input kinds below.
`tests/test_native_e2e_parity.py` extends the same comparison to the surfaces
that orchestrate the kernels: the xarray adapter, threaded and distributed Dask,
the CLI, and `fit_diagnostics`.

### Dispatch routing

Which path an input takes, and why. Every row is asserted by
`tests/test_native_parity_registry.py` against the documented behavior:

| Case | Path | Why |
|---|---|---|
| `plain_daily_series` | Rust | a plain, contiguous float64 ndarray is the kernel's own input|
| `strided` | Rust | a non-contiguous float64 rainfall series reaches the PCI dispatch guard unchanged, then is copied before the GIL is released |
| `float32` | Python | NumPy fits a float32 series in float32, so the fit stays where it was |
| `masked` | Python | a masked array is not a plain float64 ndarray |
| `year_varying_parameters` | Python CDF | only a caller can pass per-step parameters; the kernel takes one per calendar step |
| `overflowing_lmoment_block` | Python fit | an infinity makes the L-moments NaN, which the two fits report differently |
| `oversized_lmoment_block` | Python fit | a finite value above `1e100` exceeds the native L-moment safety bound |

Pearson's single-series goodness-of-fit check also batches sorted calibration
CDFs and KS D statistics in `pearson::pearson_ks_statistics`. Python keeps the
validation, critical values, candidate exact p-values, warnings, and logging.
Only a D statistic below the critical value by the full parity envelope
(`1e-10 * (1 + abs(critical))`, plus dtype epsilon) skips the SciPy check;
candidates and near-boundary statistics use the original oracle decision.
Spatial goodness-of-fit remains Python. Extension errors propagate before the
reference check's exception handler and are never retried.

### Parity tolerance

The contract is `rtol = atol = 1e-10` with matching NaN positions, per kernel,
and `scripts/native_parity_maxima.py` measures what each entry actually deviates
by on its family's fixed sample; the Hypothesis draws in the suite are held to
the same tolerance but are not part of the table. The `test-native` job records
the table in its summary once the consolidated suite passes, so the Linux x86-64
leg and a developer's machine both report their own numbers; the values below
are measured on macOS arm64
(`uv run python scripts/native_parity_maxima.py`, extension built):

| Entry | Kernels | Max absolute | Max relative |
| --- | --- | --- | --- |
| `eddi` | `hastings_inverse_normal`, `tukey_probabilities` | 0.000e+00 | 0.000e+00 |
| `eddi_spatial_block` | `hastings_inverse_normal`, `tukey_probabilities` | 0.000e+00 | 0.000e+00 |
| `fire_drought_code` | `drought_code` | 0.000e+00 | 0.000e+00 |
| `fire_duff_moisture_code` | `duff_moisture_code` | 0.000e+00 | 0.000e+00 |
| `fire_ffmc` | `ffmc` | 0.000e+00 | 0.000e+00 |
| `fire_kbdi` | `kbdi` | 0.000e+00 | 0.000e+00 |
| `fit_diagnostics` | `gamma_parameters` | 0.000e+00 | 0.000e+00 |
| `flood_api` | `antecedent_precipitation_index` | 0.000e+00 | 0.000e+00 |
| `flood_edi` | `edi`, `effective_precipitation` | 8.527e-14 | 7.883e-15 |
| `flood_flood_index` | `effective_precipitation`, `flood_index` | 2.220e-16 | 1.096e-15 |
| `flood_pe` | `effective_precipitation` | 4.547e-13 | 2.781e-16 |
| `hargreaves` | `hargreaves` | 4.441e-16 | 1.975e-16 |
| `palmer_pdsi` | `palmer_k_prime`, `palmer_pdi`, `palmer_raw_zindex`, `palmer_water_balance` | 0.000e+00 | 0.000e+00 |
| `palmer_scpdsi` | `palmer_k_prime`, `palmer_raw_zindex`, `palmer_water_balance`, `palmer_wells`, `scpdsi_duration_factors` | 0.000e+00 | 0.000e+00 |
| `pci` | `pci` | 0.000e+00 | 0.000e+00 |
| `pearson_ks_statistics` | `pearson_ks_statistics` | 3.331e-16 | 1.830e-15 |
| `penman_monteith` | `fao56_eto` | 8.882e-16 | 4.337e-16 |
| `percentage_of_normal` | `pnp_normals`, `pnp_percentages` | 0.000e+00 | 0.000e+00 |
| `pm_eto_intermediates` | `pm_eto` | 0.000e+00 | 0.000e+00 |
| `spei_gamma` | `gamma_parameters`, `gamma_probabilities`, `norm_ppf` | 2.220e-16 | 3.457e-16 |
| `spei_loglogistic` | `loglogistic_cdf`, `loglogistic_parameters` | 0.000e+00 | 0.000e+00 |
| `spi_gamma` | `gamma_parameters`, `gamma_probabilities`, `norm_ppf` | 1.998e-15 | 4.765e-14 |
| `spi_gamma_mean_zero` | `gamma_parameters`, `gamma_probabilities`, `norm_ppf` | 1.998e-15 | 4.765e-14 |
| `spi_gamma_spatial_block` | `gamma_parameters`, `gamma_probabilities`, `norm_ppf` | 7.105e-15 | 1.868e-13 |
| `spi_pearson` | `pearson_cdf`, `pearson_ks_statistics`, `pearson_parameters` | 8.604e-15 | 6.822e-14 |
| `thornthwaite` | `thornthwaite` | 2.842e-14 | 2.208e-16 |

The largest measured absolute deviation is 4.547e-13 (`flood_pe`), more than two
orders of magnitude inside `atol`. The arm64 FMA behavior described above is why the
same entry can report a smaller deviation on an x86-64 leg, and why a deviation
recorded in CI is the evidence for the platform it ran on rather than a new
tolerance.

**Migration policy.** Port expensive numerical kernels, hot loops, and
algorithms that benefit materially from native execution. Keep in Python:
xarray orchestration, CF metadata, validation at the public API boundary,
provenance, the CLI, logging and warnings, file and network I/O, and all other
user-facing behavior. A port reproduces the Python numerics, including their
NaN, zero, and edge-case semantics; it does not improve or change them.

**Building.** Hatchling remains the PEP 517 backend, so the published sdist and
`py3-none-any` wheel are pure Python, install without a Rust toolchain, and run the
Python implementations. Developers with a Rust toolchain build the extension in place:

```bash
uv run maturin develop --release   # builds src/climate_indices/_native.*.so
uv run python -c "import climate_indices._native"
cargo test --workspace
```

`maturin develop` installs the package in editable mode and copies the compiled
extension into `src/climate_indices/`, where it stays importable after `uv sync`
reinstalls the project; delete the `_native.*` file to return to pure Python.

**Packaging.** maturin also builds the published binary wheels, and they ship in
the same release as the sdist and the pure wheel. They are abi3 (`abi3-py310`), so
one wheel per platform validates on every supported Python, and there are five:
manylinux_2_28 `x86_64` and `aarch64`, macOS arm64 and `x86_64`, and Windows
`x86_64`. pip prefers a matching platform wheel over `py3-none-any`; a platform
with neither installs the pure wheel and runs the Python implementations, which is
the supported install path without Rust — installing from source needs no
toolchain either, because hatchling never invokes `cargo`. musllinux is not
published. Which wheel is installed is visible from
`python -c "import climate_indices._native"`, and any `ImportError` from that
import falls back to Python rather than failing. The full comparison of the
options and the failure modes are in
[ADR-0018](docs/adr/0018-optional-rust-packaging.md).

**Rust in CI.** `.github/workflows/unit-tests-workflow.yml` runs three jobs on
every event the workflow handles (pull request, push to `main`, merge group,
schedule, and manual dispatch), next to the pure-Python legs. Those legs install
no Rust toolchain, so they keep proving the fallback.

- `rust`: `cargo fmt --all -- --check`, `cargo clippy --workspace --all-targets
  -- -D warnings`, `RUSTDOCFLAGS="-D warnings" cargo doc -p climate-indices
  --no-deps`, and `cargo test --workspace` (`PYO3_PYTHON` points the
  `climate-py` test binary at the project interpreter).
- `test-native`: `maturin develop --release`, then the same core pytest command
  as `test`, on the boundary legs (oldest and newest Python on Linux, newest on
  macOS). It sets `CLIMATE_INDICES_REQUIRE_NATIVE=1`, which makes the native test
  modules raise on a missing extension instead of skipping, so a broken build
  cannot silently drop the parity suite. `tests/test_native_parity.py`,
  `tests/test_native_parity_distributions.py`,
  `tests/test_native_parity_fire.py`, and `tests/test_native_parity_palmer.py`
  have no other skip. The Python 3.14-only context-aware-warnings routing check lives in
  `tests/test_native_backend.py` and skips on the 3.10 leg. It then runs the consolidated
  parity suite as a named step, so a leg that collected nothing from it fails visibly
  rather than passing silently, and writes the measured parity maxima
  (`scripts/native_parity_maxima.py`, the table above) into the job summary.
- `native-wheel`: `maturin build --release` on Linux and macOS at both boundary
  Pythons, plus a Windows smoke build on the newest. Each wheel is installed
  into a fresh venv and imported from outside the checkout, and SPI must reach
  the bundled `gamma_parameters`, `gamma_probabilities`, and `norm_ppf` kernels.
  This validates the developer build; the published wheel matrix, its abi3 tag,
  and the installation checks for each artifact run in `.github/workflows/release.yml`
  ([ADR-0018](docs/adr/0018-optional-rust-packaging.md)).

`rust-toolchain.toml` pins the compiler so a new stable clippy lint cannot fail
`-D warnings` on an unrelated change. The workspace uses edition 2024, so the
minimum supported Rust is 1.85. To bump the pin, change `channel`, run the clippy
command above, and fix any new lints in the same change.

## Source Code Organization

```
climate_indices/
├── src/climate_indices/          # Main package directory
│   ├── __init__.py               # Public API exports
│   ├── __main__.py               # Full-featured CLI entry point
│   ├── typed_public_api.py       # Strict mypy-compliant API
│   ├── xarray_adapter.py         # Modern xarray interface
│   ├── validation.py             # Shared input validation facade
│   ├── indices.py                # Legacy numpy API (STABLE)
│   ├── compute.py                # Core computation algorithms
│   ├── palmer.py                 # Palmer drought indices
│   ├── eto.py                    # PET: Thornthwaite & Hargreaves
│   ├── lmoments.py               # L-moments for Pearson fitting
│   ├── utils.py                  # Utility functions
│   ├── logging_config.py         # Structured logging setup
│   ├── exceptions.py             # Exception hierarchy
│   └── performance.py            # Performance metrics
│
├── tests/                        # Test suite
│   ├── conftest.py               # Shared fixtures (session-scoped)
│   ├── test_indices.py           # Legacy API tests
│   ├── test_xarray_adapter.py    # Modern API tests
│   ├── test_validation.py        # Validation facade tests
│   ├── test_compute.py           # Computation tests
│   ├── test_property_based.py    # Property-based invariant tests
│   ├── test_backward_compat.py   # Backward compatibility suite
│   ├── test_exceptions.py        # Exception handling tests
│   ├── test_observability.py     # Logging behavior and lifecycle tests
│   ├── test_metadata_validation.py # CF metadata validation tests
│   ├── test_benchmark_*.py       # Performance regression tests
│   └── fixture/                  # Test data (numpy arrays, JSON)
│
├── docs/                         # Documentation
│   ├── conf.py                   # Sphinx configuration
│   ├── index.md                  # Main Sphinx doc (ReadTheDocs)
│   ├── reference.md              # API reference (autodoc)
│   ├── release-process.md        # Maintainer release runbook
│   └── *.md                      # AI-readable project docs
│
├── .github/workflows/            # CI/CD pipelines
│   ├── unit-tests-workflow.yml   # Test matrix (Python 3.10-3.14)
│   ├── release.yml               # Automated PyPI releases
│   └── benchmarks.yml            # Performance tracking
│
├── pyproject.toml                # PEP 517 build config + tool settings
├── uv.lock                       # Reproducible dependency lock
├── Dockerfile                    # Container image definition
├── README.md                     # GitHub landing page
└── CONTRIBUTING.md               # Development guidelines
```

### Critical Directories
- **`src/climate_indices/`**: Production code for the core indices and the fire subsystem
- **`tests/`**: Test suite and fixture data
- **`docs/`**: Sphinx + MyST Markdown project documentation
- **`.github/workflows/`**: CI/CD automation (3 workflows)

## Data Flow and Computation Patterns

### SPI/SPEI Computation Flow

```
┌─────────────────────────────────────────────────────────────────┐
│  User Input: xarray.DataArray or numpy.ndarray                 │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  Input Validation & Type Detection                              │
│  - Check time dimension exists and is monotonic                 │
│  - Validate calibration period length (≥30 years recommended)   │
│  - Check for excessive NaNs                                     │
│  - Detect input type (numpy vs xarray vs Dask)                 │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  Coordinate Alignment (xarray only)                             │
│  - Align precipitation and PET on time coordinate (SPEI)        │
│  - Warn if alignment drops time steps                           │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  Temporal Scaling (compute.prepare_scaled)                      │
│  - Rolling sum over N months/days                               │
│  - Handles NaN propagation                                      │
│  - Output: scaled_values (same shape as input)                  │
│  - Except an all-missing input: returned unreshaped, so a       │
│    2-D all-missing input comes back 1-D                         │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  Distribution Fitting (per calendar month/day)                  │
│  - Gamma: alpha, beta parameters (method of moments)            │
│  - Pearson: loc, scale, skew parameters (L-moments)             │
│  - Fit on calibration period only                               │
│  - Track probability of zero separately                         │
│  - One seam (compute.fit_and_standardize) normalizes fitting    │
│    parameters and dispatches on the distribution                │
│  - SPI falls back from a failed Pearson fit to gamma, SPEI      │
│    propagates the failure (fallback_to_gamma)                   │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  CDF Transformation                                             │
│  - Apply fitted CDF to scaled values                            │
│  - Transform to uniform [0,1] distribution                      │
│  - Uses scipy.stats.gamma or scipy.stats.pearson3              │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  Inverse Normal Transformation                                  │
│  - Apply scipy.stats.norm.ppf() to CDF values                   │
│  - Output: SPI/SPEI values (standard normal distribution)       │
│  - Handle edge cases (CDF=0 → -3.09, CDF=1 → 3.09)             │
└───────────────────────┬─────────────────────────────────────────┘
                        │
                        ▼
┌─────────────────────────────────────────────────────────────────┐
│  Output Formatting                                              │
│  - xarray: Preserve input structure + add CF metadata           │
│  - numpy: Return array with same shape as input                │
│  - Dask: Return lazy Dask array (compute on demand)            │
└─────────────────────────────────────────────────────────────────┘
```

### Parallelization Strategies

#### 1. CLI Multiprocessing (\_\_main\_\_.py)
```python
# Splits data across lat/lon dimensions
# Each worker processes a spatial subset
# Shared memory arrays for inputs/outputs
# Pool size: CPU count - 1
with multiprocessing.Pool(processes=num_workers) as pool:
    pool.map(_apply_along_axis, chunk_params)
```

#### 2. Dask Lazy Evaluation (xarray_adapter.py)
```python
# Dask-backed xarray processing
# Time dimension: single chunk (required)
# Spatial dimensions: chunked (e.g., 50x50 cells)
# Computation triggered by .compute() or .load()
result_da = xr.apply_ufunc(
    spi_computation_fn,
    precip_da,
    dask='parallelized',
    output_dtypes=[float]
)
```

## Testing Architecture

### Test Organization (66 Test Files)
```
tests/
├── conftest.py                      # Session-scoped fixtures
│
├── Core Functionality Tests
│   ├── test_indices.py              # Legacy numpy API tests
│   ├── test_compute.py              # Core algorithm tests
│   ├── test_xarray_adapter.py       # Modern API tests (EXPANDED)
│   ├── test_typed_public_api.py     # Strict typing tests (NEW)
│   └── test_eto.py                  # PET computation tests
│
├── Quality Assurance Tests
│   ├── test_backward_compat.py      # API stability tests
│   ├── test_xarray_equivalence.py   # numpy ↔ xarray parity
│   ├── test_property_based.py       # Hypothesis invariant tests
│   └── test_type_checking.py        # mypy runtime validation
│
├── Validation and Error Handling
│   ├── test_exceptions.py           # Exception hierarchy tests
│   ├── test_input_validation.py     # Input validation tests
│   ├── test_validation.py           # Validation facade tests
│   ├── test_metadata_validation.py  # CF metadata tests
│   ├── test_computation_errors.py   # Error condition tests
│   ├── test_data_quality_warnings.py # Warning behavior tests
│   └── test_input_type_detection.py # Type detection tests
│
├── Observability Tests
│   ├── test_observability.py        # Logging config, lifecycle events
│   └── test_performance_metrics.py  # Performance tracking tests
│
├── Performance Tests (Benchmarks)
│   ├── test_benchmark_overhead.py   # xarray vs numpy overhead
│   ├── test_benchmark_chunked.py    # Dask chunking strategies
│   ├── test_benchmark_fire.py       # Fire-index throughput, sizing, and guards
│   └── test_benchmark_memory.py     # Memory usage profiling
│
└── Regression Tests
    ├── test_palmer.py               # Palmer indices regression
    ├── test_utils.py                # Utility function tests
    └── test_zero_precipitation_fix.py # Specific bug fix test
```

### Test Fixtures (conftest.py)
**Session-scoped fixtures** for performance:
- **Numpy arrays**: `precips_mm_monthly`, `temps_celsius`, `pet_thornthwaite_mm`
- **xarray DataArrays**: `sample_monthly_precip_da`, `gridded_monthly_precip_3d`, `dask_monthly_precip_1d`
- **Edge case fixtures**: `zero_inflated_precip_da`, `leading_nan_block_da`, `non_monotonic_time_da`
- **Benchmark fixtures**: `bench_monthly_precip_np`, `bench_monthly_precip_da`, `bench_gridded_precip_da`
- **Constants**: `_CALIBRATION_YEAR_START_MONTHLY`, `_DATA_YEAR_START_MONTHLY`, `_LATITUDE_DEGREES`

### Test Coverage Targets
- **Overall**: >90% line coverage
- **Critical modules**: `compute.py`, `indices.py`, `xarray_adapter.py` → 100%
- **Exception paths**: All custom exceptions tested with context attributes
- **Property-based**: Mathematical invariants (e.g., SPI mean ≈ 0, std ≈ 1)

### Running Tests
```bash
# All tests (excluding benchmarks)
uv run pytest

# Include benchmarks
uv run pytest -m benchmark

# With coverage
uv run pytest --cov=src --cov-report=html

# Specific test file
uv run pytest tests/test_xarray_adapter.py -v

# Property-based tests only
uv run pytest tests/test_property_based.py
```

## Deployment and CI/CD

### GitHub Actions Workflows

#### 1. unit-tests-workflow.yml
**Trigger**: Pull request and push to main, merge group, weekly schedule, manual dispatch
```yaml
Matrix (core suite):
  test (every event):
    - ubuntu-latest: Python 3.10, 3.11, 3.12, 3.13, 3.14
    - macos-latest: Python 3.14
  test-full (everything except pull requests):
    - macos-latest: Python 3.10
Steps:
  1. Checkout code
  2. Setup Python + uv
  3. uv sync --locked --dev
  4. Run pytest -n auto
Rust jobs (every event; see "Rust in CI" under Optional Rust Backend):
  rust, test-native, native-wheel
```

#### 2. release.yml
**Trigger**: Git tag push (vX.Y.Z)
```yaml
Steps:
  1. Checkout code
  2. Build sdist and wheel (hatchling)
  3. Run twine package checks
  4. Publish to PyPI (trusted publishing via OIDC)
```

#### 3. benchmarks.yml
**Trigger**: Pull request to `main`, manual dispatch
```yaml
Steps:
  1. Checkout code
  2. Setup Python + uv
  3. uv sync --group dev
  4. Run pytest -m benchmark --benchmark-enable --benchmark-json
  5. Upload benchmark artifact
  6. Compare against baseline when available
```

### Docker Container

**Base Image**: `python:3.14-slim`
**Build Strategy**: Multi-stage (builder + production)

```dockerfile
# Builder stage: Install dependencies
FROM python:3.14-slim AS builder
COPY --from=ghcr.io/astral-sh/uv:latest /uv /bin/uv
WORKDIR /app
COPY pyproject.toml uv.lock ./
COPY README.md LICENSE ./
COPY src/ ./src/
RUN uv sync --frozen --no-dev

# Production stage: Copy venv + source
FROM python:3.14-slim
RUN apt-get update && apt-get install -y \
    libhdf5-dev libnetcdf-dev
COPY --from=builder /app/.venv /app/.venv
COPY src/ ./src/
USER climate  # Non-root user
ENTRYPOINT ["python", "-m", "climate_indices"]
```

**Usage**:
```bash
docker build -t climate_indices:X.Y.Z .
docker run -v $(pwd)/data:/data climate_indices:X.Y.Z \
    --index spi --scales 6 --netcdf_precip /data/precip.nc \
    --var_name_precip prcp --output_file_base /data/spi
```

### PyPI Distribution

**Package Name**: `climate_indices`
**Installation**: `pip install climate_indices` or `uv pip install climate_indices`
**Artifacts**:
- **Source distribution** (`climate_indices-X.Y.Z.tar.gz`)
- **Wheel** (`climate_indices-X.Y.Z-py3-none-any.whl`)

**Excludes from package** (`pyproject.toml`):
- `tests/`, `docs/`, `notebooks/`, `assets/`, `.github/`, `.venv/`, cache directories

## Key Design Decisions and Trade-offs

### 1. Dual API (numpy vs xarray)
**Decision**: Maintain both numpy and xarray APIs.

**Rationale**:
- **Numpy**: Minimal dependencies, direct array manipulation, backward compatibility
- **xarray**: Labeled dimensions, CF metadata, Dask integration, modern workflow

**Trade-off**: Code duplication risk mitigated by having xarray API call numpy implementation internally.

See [ADR-0001](docs/adr/0001-dual-numpy-xarray-api.md).

### 2. Multiprocessing vs Dask
**Decision**: Use multiprocessing in CLI, Dask in xarray API.

**Rationale**:
- **Multiprocessing**: Predictable memory usage, no Dask dependency for CLI users
- **Dask**: Lazy evaluation, better integration with xarray ecosystem, dynamic scheduling

**Trade-off**: Separate parallelization logic in CLI and xarray layers.

See [ADR-0002](docs/adr/0002-multiprocessing-cli-dask-xarray.md).

### 3. Time Dimension Chunking Constraint

**Decision**: For fitting and stateful indices, Dask arrays MUST have time as a single chunk. The two PET entry points are the exception: they accept a split `time` and rechunk it internally, because their kernels are period-based arithmetic rather than fits.

**Rationale**: Climate indices require access to full time series for distribution fitting.

**Enforcement**: `validation.validate_dask_chunks()` validates chunking on the adapter path and raises `CoordinateValidationError` if violated; the PET entry points skip it.

See [ADR-0003](docs/adr/0003-dask-time-dimension-single-chunk.md) for the full rationale.

### 4. Exception Hierarchy with Context
**Decision**: Custom exceptions with context attributes instead of plain `ValueError`.

**Rationale**: AI agents and users need structured error information for debugging.

**Example**:
```python
raise DistributionFittingError(
    "Gamma fitting failed",
    distribution_name="gamma",
    input_shape=(480, 5, 6),
    parameters={"alpha": "NaN", "beta": "NaN"},
    suggestion="Try Pearson Type III distribution"
)
```

### 5. Property-Based Testing
**Decision**: Use Hypothesis for mathematical invariant testing.

**Rationale**: Traditional unit tests miss edge cases; property-based tests generate adversarial inputs.

**Example Properties**:
- SPI output has mean ≈ 0, standard deviation ≈ 1 over calibration period
- SPI is monotonic with respect to input precipitation
- PET is always non-negative

## Performance Considerations

### Bottlenecks
1. **Distribution fitting**: scipy.stats fitting functions (CPU-bound)
2. **CDF transformation**: scipy.stats.cdf() calls (CPU-bound)
3. **Temporal scaling**: Rolling sum over large arrays (memory-bound)
4. **I/O**: NetCDF reading for large gridded datasets (I/O-bound)

### Optimization Strategies
1. **Vectorization**: Numpy broadcasting instead of loops
2. **Shared memory**: CLI uses multiprocessing.Array for zero-copy
3. **Lazy evaluation**: Dask defers computation until .compute()
4. **Chunking**: Spatial chunks in Dask, single time chunk for fitting and stateful indices
5. **Caching**: Distribution fitting parameters can be fitted once and passed
   back in through the `fitting_params` argument of `indices.spi()`

### Benchmark Results (Typical)
| Operation | Input Size | Execution Time | Memory |
|-----------|-----------|----------------|--------|
| SPI-6 (numpy, 1D) | 480 months | ~50 ms | <10 MB |
| SPI-6 (xarray, 1D) | 480 months | ~60 ms | <15 MB |
| SPI-6 (xarray, 3D) | 480×20×20 | ~2 sec | ~100 MB |
| SPI-6 (Dask, 3D) | 480×100×100 | ~20 sec | ~500 MB |

## Security Considerations

### Input Validation
- **NetCDF files**: Dimension checks, coordinate validation
- **User inputs**: Scale range [1-72] monthly / [1-2196] daily, year validation
- **Path sanitization**: No path traversal in CLI file arguments

### Dependency Security
- **uv lock**: Pinned dependencies with checksums
- **GitHub Actions**: Pinned action versions with SHA hashes
- **Container**: Non-root user execution

### No Network Communication
Library has no network dependencies; all data is file-based.

---

**Next Steps**: See [development-guide.md](docs/development-guide.md) for setup instructions and [deployment-guide.md](docs/deployment-guide.md) for CI/CD details; the generated API reference is published in the [Sphinx reference page](docs/reference.md).

## docs/development-guide.md

# Development Guide

## Prerequisites

### Required Software
- **Python**: 3.10, 3.11, 3.12, 3.13, or 3.14
- **uv**: Modern Python package installer ([https://github.com/astral-sh/uv](https://github.com/astral-sh/uv))
- **Git**: Version control
- **Make**: Build automation (optional, for Sphinx docs)

### System Dependencies (for NetCDF support)
```bash
# macOS
brew install hdf5 netcdf

# Ubuntu/Debian
sudo apt-get install libhdf5-dev libnetcdf-dev

# Fedora/RHEL
sudo dnf install hdf5-devel netcdf-devel
```

### Optional: the native Rust extension

Everything in the library runs in pure Python, and `uv sync` installs it that way.
The native extension is an optional accelerator: when it is present, the dispatch
described in [architecture.md](docs/architecture.md#optional-rust-backend) sends prepared
float64 arrays to Rust kernels instead of the Python reference. You do **not** need
Rust to install, import, run, or test `climate_indices`; the native suites
(`tests/test_native_*.py`, covering kernel and PET parity plus the extension smoke
checks) skip without it. A change to a ported kernel is not actually compared
against the Rust path until you build the extension, so build it before claiming
parity. No Rust experience is required: [installing the Rust
toolchain](#install-the-rust-toolchain) below is a one-time, two-command step.

## Installation

### 1. Clone Repository
```bash
git clone https://github.com/monocongo/climate_indices.git
cd climate_indices
```

### 2. Install uv (if not already installed)
```bash
# macOS/Linux
curl -LsSf https://astral.sh/uv/install.sh | sh

# Windows
powershell -c "irm https://astral.sh/uv/install.ps1 | iex"

# Or via pip
pip install uv
```

### 3. Create Virtual Environment and Install Dependencies
```bash
# Sync all dependencies (core + dev + test)
uv sync --group dev

# Activate virtual environment
source .venv/bin/activate  # macOS/Linux
.venv\Scripts\activate     # Windows
```

### 4. Install in Development Mode
```bash
# Install package in editable mode
uv pip install -e .
```

## Development Workflow

### Running Tests
```bash
# All tests (excludes benchmarks by default)
uv run pytest

# With verbose output
uv run pytest -v

# Specific test file
uv run pytest tests/test_xarray_adapter.py

# Include benchmarks
uv run pytest -m benchmark

# With coverage report
uv run pytest --cov=src --cov-report=html
# Open htmlcov/index.html in browser
```

### Code Quality Checks

#### Linting with Ruff
```bash
# Check for issues
ruff check src/ tests/

# Auto-fix issues
ruff check --fix src/ tests/

# Format code
ruff format src/ tests/
```

#### Type Checking with Mypy

Run mypy on Python 3.12, the stack CI's lint job uses. The package-wide `strict`
setting in `pyproject.toml` applies automatically, so no `--strict` flag is
needed.

```bash
# Check the package and the type assertions pytest cannot see
uv run --python 3.12 mypy src/ tests/test_type_checking.py

# Check a specific module
uv run --python 3.12 mypy src/climate_indices/typed_public_api.py
```

The `mypy` commands require a Python 3.11+ dependency stack. The Python 3.10
stack pins numpy to 2.2.x, whose stubs raise findings that no single source tree
can satisfy on both stacks, so the release workflow skips its mypy leg on 3.10.

### Building Documentation

#### Sphinx Documentation
```bash
cd docs
make html
# Open _build/html/index.html
```

#### Clean Build
```bash
cd docs
make clean
make html
```

### Building Distribution Packages
```bash
# Build wheel and sdist
uv run python -m build

# Build a local binary wheel for this platform (needs a Rust toolchain). On Linux
# this is not the published manylinux_2_28 wheel: release.yml builds that inside the
# container named in its matrix, which is what sets the glibc floor (ADR-0018)
uv run maturin build --release --out dist

# Output in dist/
# - climate_indices-X.Y.Z-py3-none-any.whl
# - climate_indices-X.Y.Z.tar.gz
# - climate_indices-X.Y.Z-cp310-abi3-<platform>.whl
```

## Project Structure

```
climate_indices/
├── src/climate_indices/     # Source code
├── tests/                    # Test suite
├── docs/                     # Documentation
├── pyproject.toml            # Build config + dependencies
├── uv.lock                   # Dependency lock file
└── .github/workflows/        # CI/CD pipelines
```

## Development Commands Reference

| Task | Command |
|------|---------|
| **Install dependencies** | `uv sync --group dev` |
| **Run tests** | `uv run pytest` |
| **Run benchmarks** | `uv run pytest -m benchmark` |
| **Coverage report** | `uv run pytest --cov=src --cov-report=html` |
| **Lint code** | `ruff check --fix src/ tests/` |
| **Format code** | `ruff format src/ tests/` |
| **Type check** | `uv run --python 3.12 mypy src/ tests/test_type_checking.py` |
| **Build docs** | `cd docs && make html` |
| **Build package** | `uv run python -m build` |
| **Update deps** | `uv lock` |

## Porting a Kernel to Rust

Adding a Rust kernel means adding a second implementation of a numerical block that
already exists in Python. The Python implementation stays, and it stays the oracle:
the port is only finished when the two agree at `rtol = atol = 1e-10` with identical
NaN positions. The architecture is decided in
[ADR 0017](docs/adr/0017-rust-core-acceleration-backend.md); the crate layout, the
already-ported seams, and the dispatch rules are in
[architecture.md](docs/architecture.md#optional-rust-backend). `crates/climate-indices/src/gamma.rs`
and its dispatch in `src/climate_indices/compute.py` are the reference implementation
of every step below.

### Install the Rust toolchain

Rust is distributed through `rustup`, a toolchain manager that plays the same role
for Rust that `uv` plays for Python. You do not need to pick a version: this
repository's `rust-toolchain.toml` pins the compiler (1.99.0) and the two components
the checks use (`rustfmt`, `clippy`), and `rustup` installs that exact toolchain
automatically the first time you run a Rust command inside the checkout. Install
`rustup` once, then let the checkout pick the version.

```bash
# macOS/Linux (official installer — recommended)
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal

# macOS with Homebrew (alternative). rustup is keg-only, so it is installed but
# not put on PATH: resolve the prefix Homebrew actually used (`/opt/homebrew` on
# Apple Silicon, `/usr/local` on Intel) and add its shim directory to your shell's
# startup file (`~/.zshrc` for zsh, `~/.bash_profile` for bash).
brew install rustup
echo "export PATH=\"$(brew --prefix rustup)/bin:\$PATH\"" >> ~/.zshrc

# Windows (PowerShell)
winget install --id Rustlang.Rustup
```

The official installer places the Rust binaries in `~/.cargo/bin` and adds that
directory to your shell's `PATH`; open a new terminal afterwards (or run
`source "$HOME/.cargo/env"`) so the commands below resolve. If you used Homebrew,
the `echo` line does the same for the Homebrew copy, also taking effect in a new
terminal. This is the step people most often miss: a Homebrew rustup without that
line leaves `cargo` and `rustup` off `PATH`, so every Rust command reports
`command not found` even though the toolchain is installed. On Windows, install the
Microsoft C++ build tools (the Visual Studio installer's "Desktop development with
C++" workload) before building: the MSVC toolchain links against them, and the build
fails at the link step without them.

Verify the install:

```bash
rustup --version
rustup show    # the toolchain in effect for the current directory (1.99.0)
cargo --version
```

`rustup show` lists the toolchains installed and names the one in effect. The first
Rust command you run inside the checkout also downloads the pinned 1.99.0 toolchain
if it is missing; to fetch it explicitly, run
`rustup toolchain install 1.99.0 --component rustfmt --component clippy`.

In this project `cargo` is the Rust build tool and package manager, `rustfmt`
formats Rust code, and `clippy` lints it; `maturin` is the Python bridge that wraps
the compiled crates into an importable extension. `maturin` is already a locked dev
dependency, so `uv sync --group dev` installs it and you run it through `uv`.

Build the extension and confirm it imports, from the repository root:

```bash
uv run maturin develop --release
uv run python -c "import climate_indices._native; print(climate_indices._native.__file__)"
```

`maturin develop` compiles the crates in `crates/`, copies the resulting `_native`
extension into `src/climate_indices/`, and installs the package in editable mode, so
the extension stays importable when `uv sync` later reinstalls the project. The first
build downloads the pinned toolchain and the crates it depends on and can take a few
minutes; later builds are incremental. To return to the pure-Python install, delete
the compiled file and start a new interpreter:

```bash
rm -f src/climate_indices/_native.*.so src/climate_indices/_native.*.pyd
```

The everyday Rust commands, for reference:

| Command | What it does |
|---------|--------------|
| `rustup show` | The toolchain this checkout pins (must be 1.99.0) |
| `rustup update` | Updates installed toolchains — do **not** use this to move the pin; change `channel` in `rust-toolchain.toml` deliberately instead |
| `cargo build` / `cargo test --workspace` | Compiles the crates / runs the crate-level Rust tests; set `PYO3_PYTHON` to the project interpreter, as in [step 8](#8-verify), so they link the right Python |
| `cargo fmt --all -- --check` | Formatting; CI fails on any diff |
| `cargo clippy --workspace --all-targets -- -D warnings` | Lints; CI fails on any warning |
| `RUSTDOCFLAGS="-D warnings" cargo doc -p climate-indices --no-deps` | Crate documentation; CI fails on any rustdoc warning |

### Crate quality gates

The published `climate-indices` crate has four required CI gates, one job each, so
each is a required check on the pull request. Run them locally with
`cargo install --locked cargo-deny@0.20.2 cargo-semver-checks@0.51.0` once, then:

| Command | What it proves |
|---------|----------------|
| `cargo package -p climate-indices --locked`, unpack into a directory **outside** this checkout, copy `rust-toolchain.toml` into it, then `cargo clippy --all-targets --all-features --locked -- -D warnings`, `cargo test --all-features --locked`, and `RUSTDOCFLAGS="-D warnings" cargo doc --all-features --no-deps --locked` | The crate that a user receives lints, tests, and documents cleanly, on the pinned compiler. Running it from the unpacked `.crate` is what catches a file the crate needs but does not ship. |
| `cargo "+{msrv}" check -p climate-indices --lib --all-features --locked`, where `{msrv}` is the crate's `rust-version` | The library builds on its declared minimum supported Rust version (1.85). `--lib` leaves out `criterion`, which only the benchmarks use. |
| `cargo semver-checks check-release -p climate-indices` | The public API is compatible with the latest version on crates.io. A deliberate break fails here, cited to the offending item. |
| `cargo deny --locked check -W unmaintained -W unsound` | The licenses, advisories, duplicate versions, and sources permitted by `deny.toml` cover everything that reaches a user — the published crate and the `climate-py` bindings the wheel compiles. Unmaintained and unsound advisories are warnings, so an upstream archive notice does not block an unrelated pull request. |

The gates read the crate's own `rust-version` and version, so moving either needs no
workflow change. `cargo semver-checks` compares against the published baseline. The
`0.0.0` placeholder had no API to compare; since `0.1.0` was published, the check has a
real baseline.

CI runs the native legs with `CLIMATE_INDICES_REQUIRE_NATIVE=1`, which turns a
missing extension into a test failure instead of a skip so a broken build cannot
silently drop the parity suite. Set the same variable locally to confirm a change
really exercises the Rust path.

### 1. Trace the Python seam

Read the Python function end to end and name the block that becomes the kernel, then
name everything that stays: validation, calibration-period resolution, data-quality
and goodness-of-fit warnings, the Pearson-to-gamma fallback, logging, zero placement,
support-limit masks, and output scaling. The seam goes where the prepared float64
values are already valid, so the kernel takes no validation decisions and raises no
`climate_indices` exceptions.

Record the semantics the kernel has to reproduce, because these are what the port is
reviewed against: operation order, what NaN means, how zeros are treated, and what a
degenerate column does (all-missing, all-zero, constant, or no positive value).

### 2. Write the kernel documentation block

The `//!` crate/domain header and the `///` block above each public kernel state:

- the Python source it ports, named by function and by the block inside it;
- inputs and outputs, including the array layout (years by columns) and which
  argument is per-column;
- zero semantics, NaN semantics, invalid-input behaviour, short-series and
  degenerate-column behaviour;
- the numerics: the expressions, the summation order, and where SciPy is being
  reproduced.

### 3. Port the kernel

Put the kernel in `crates/climate-indices/src/`, in the module that owns the index
family, and declare that module in `crates/climate-indices/src/lib.rs`
(`pub mod <family>;`) so `climate-py` can reach it by module path. The crate is pure
Rust: no PyO3, no NumPy bindings, no Python exceptions. Where SciPy evaluates a
special function, port the routine SciPy evaluates into `special/` (see the Cephes
ports for `igam`, `igamc`, `ndtri`, `ndtr`, and `lgam`) rather than calling a
general crate — a generic implementation does not hold parity in the tails.
Reproduce the Python arithmetic rather than improving it, and do not add a
fast-math or fused kernel: `[profile.release]` keeps IEEE semantics on purpose.

### 4. Add the binding and the stub

Add the `#[pyfunction]` to `crates/climate-py/src/lib.rs`, converting arrays at the
boundary and returning only `climate-indices` types across it; that crate holds no
algorithm. Register it in the `#[pymodule]` in that file with
`m.add_function(wrap_pyfunction!(<kernel>, m)?)?`: a `#[pyfunction]` that is not
registered is never exposed to Python. Add the matching signature to
`src/climate_indices/_native.pyi`, annotating each array argument with its dtype:
float64 values, bool validity masks, int64 indices. Existing bindings call
`checked_copy`, which rejects unaligned arrays, copies empty and caller-owned
storage before `py.detach`, and wraps the view conversion. Simple elementwise
stages such as PNP percentages and Palmer CAFEC use `checked_view` without releasing
the GIL, avoiding input-sized copies.

### 5. Add dispatch and routing

The index dispatch helper is named `_native_<kernel>`, lives next to the Python
implementation it ports in `compute.py`, and is reached by its public callers in
`indices.py` as `compute._native_<kernel>`. The fire recurrences use a different
pattern: `fire/_native.py` defines `moisture_code_recurrence` and `kbdi_recurrence`,
which return a callable for the shared daily runner rather than a computed result.
Two caller patterns are in use: a helper returns `None` when a
routing condition fails and its caller runs the Python path (as
`_native_gamma_probabilities` does), or the helper carries no guard of its own and is
called only after the caller's `_native_float64` guard has passed (as
`_native_gamma_parameters` does). Dispatch takes the Rust path only when every
routing condition holds:

- the extension imported (`_native is not None`) and actually exporting the kernel —
the two EDDI helpers probe with `hasattr(_native, "<kernel>")`, so an extension
built before the port falls back to Python;
- a plain, aligned float64 `ndarray` whose layout the kernel supports — use the
existing `_native_float64` / `_native_lmoment_input` guards rather than writing new
ones;
- prepared arguments the kernel accepts, e.g. parameters that are one per calendar
step (and cell); caller-supplied parameters that vary by year stay in Python;
- NumPy floating-point errors ignored (`np.errstate(all="ignore")`): warnings,
exceptions, callbacks, logging, or printing keep the Python path, because the
kernels do not implement NumPy's reporting policies.

A kernel is never retried in Python after a Rust runtime error: the exception
propagates. A missing extension is not an error, it is the pure-Python install.

### 6. Write the parity tests

Extend the suite that covers the index family — `tests/test_native_parity.py` for the
SPI/EDDI kernels, `tests/test_native_parity_distributions.py` for the distribution
fits, `tests/test_native_parity_fire.py` for the fire recurrences,
`tests/test_native_parity_palmer.py` for the Palmer family — running the same
computation twice through the public or compute-level API: once with `_native`
replaced by the `_Recorder` around the extension, once with `_native` set to `None`,
comparing at `rtol = atol = 1e-10` with matching NaN positions. The compute suites
patch `compute._native`; the fire suite patches `climate_indices.fire._native._native`,
and the Palmer suite patches `palmer._native`.
Tests that establish invocation assert the recorded call set, so a run that never
reached Rust fails instead of passing Python against Python; result-only comparisons
discard the call set, and intentional-fallback cases assert an empty one. The EDDI
guards probe with `hasattr`, which the compute recorder also records, so a stale
extension that imports but predates an EDDI export records the probe, falls back to
Python, and still satisfies the expected EDDI call set: the probe is never evidence
of invocation, and the extension must be rebuilt after adding a kernel. Cover the edge
cases the doc block promises: all-missing and all-zero columns, a constant column,
supplied parameters, and masked or unaligned inputs that must fall back.

A test that targets the Python reference itself uses the `python_backend` fixture
from `tests/conftest.py`, which sets the `_native` attribute of every dispatch module
(`compute`, `eto`, `pm_eto`, `fire`, `flood`, and `palmer`) to `None` for the test; a new
dispatch module is added to that fixture in the same change. Loosening a tolerance needs a measured justification in the
ticket (maximum absolute and relative error, where it occurs, which primitive
diverges, and whether reproducing the Python method closes the gap); never edit a
fixture or a reference test to match Rust.

### 7. Regenerate the documentation bundles

If the change touched a file listed in `SUMMARY_FILES` or `FULL_FILES` in
`scripts/generate_llms_txt.py` — this guide and `docs/architecture.md` among them —
regenerate and commit the bundles before verifying, because the verification run
below includes the bundle-sync test:

```bash
uv run scripts/generate_llms_txt.py
uv run pytest tests/test_review_scripts.py
```

### 8. Verify

```bash
cargo fmt --all -- --check
cargo clippy --workspace --all-targets -- -D warnings
RUSTDOCFLAGS="-D warnings" cargo doc -p climate-indices --no-deps
PYO3_PYTHON="$(uv run python -c 'import sys; print(sys.executable)')" cargo test --workspace
uv run maturin develop --release
uv run python -c "import climate_indices._native"
uv run ruff check src/ tests/
uv run ruff format --check src/ tests/
uv run --python 3.12 mypy src/ tests/test_type_checking.py
CLIMATE_INDICES_REQUIRE_NATIVE=1 uv run pytest -n auto   # with the extension built
uv run pytest -m validation -n auto  # external reference suites, per the RUST ticket
rm -f src/climate_indices/_native.*.so src/climate_indices/_native.*.pyd && uv run pytest -n auto   # pure-Python fallback
uv run --extra docs sphinx-build -E -b html -W --keep-going docs docs/_build/html
uv run --extra docs sphinx-build -E -b doctest docs docs/_build/doctest
```

The fallback run matters as much as the native one: nothing user-facing may require
the extension. The native run above sets `CLIMATE_INDICES_REQUIRE_NATIVE=1`, as CI's
native legs do, so a missing or stale extension fails instead of silently skipping
every parity test.

## Coding Standards

### Python Style
- **Line length**: 120 characters
- **Python version**: 3.10 minimum (avoid syntax and stdlib added after 3.10)
- **Type hints**: Required for all functions
- **Docstrings**: Google-style for all public functions
- **Imports**: stdlib → third-party → local (enforced by ruff)

### Import Order Example
```python
# Standard library
from __future__ import annotations
import os
from typing import Optional

# Third-party
import numpy as np
import xarray as xr

# Local
from climate_indices import compute, exceptions
```

### Testing Standards
- **Framework**: pytest with fixtures
- **Coverage target**: >90%
- **Property-based tests**: Use hypothesis for invariants
- **Markers**: Use `@pytest.mark.benchmark` for performance tests

### Icechunk I/O tests

The optional Icechunk integration tests live in `tests/test_icechunk_io.py` and are
marked `icechunk`. They cover round-trip fidelity, SPI and SPEI parity across Icechunk,
NetCDF, and Zarr inputs, session and snapshot semantics, and eager, threaded, and
distributed writes. Install the extra and run just that suite:

```bash
uv sync --extra icechunk   # Python 3.11 or later; `--dev` installs it too
uv run pytest -m icechunk
```

Without the extra (Python 3.10, or a base install) the module skips. Every CI matrix leg
on Python 3.11 and later installs the `dev` group, so those legs run the suite already;
the marker is for selecting it locally.

One test is opt-in in a second way: `test_r2_round_trip` writes to the shared Cloudflare
R2 bucket, so it is marked `validation` and excluded from the default selection. It needs
both an explicit switch and credentials in the environment:

```bash
R2_TEST=1 R2_ACCOUNT_ID=<account> R2_ACCESS_KEY_ID=<key> R2_SECRET_ACCESS_KEY=<secret> \
  uv run pytest -m validation tests/test_icechunk_io.py
```

It creates a fresh `R2_PREFIX/pytest-<id>` repository (default prefix `examples/icechunk`,
default bucket `climate-indices`) and leaves that snapshot for inspection, so it never
touches another object in the bucket. Delete the prefix when it is no longer needed.

### Git Workflow
1. Create feature branch from `main`
2. Make changes with descriptive commits
3. Run tests locally (`uv run pytest`)
4. Push and create pull request
5. CI runs full test matrix
6. Merge after review

## Troubleshooting

### Common Issues

#### NetCDF Import Errors
```bash
# Install system dependencies
brew install hdf5 netcdf  # macOS
sudo apt-get install libhdf5-dev libnetcdf-dev  # Ubuntu
```

#### uv Not Found
```bash
# Add to PATH
export PATH="$HOME/.cargo/bin:$PATH"  # Add to ~/.bashrc or ~/.zshrc
```

#### Test Failures
```bash
# Clear pytest cache
pytest --cache-clear

# Reinstall dependencies
uv sync --group dev --force
```

---

See [CONTRIBUTING.md](https://github.com/monocongo/climate_indices/blob/main/CONTRIBUTING.md) for detailed contribution guidelines.

## docs/contribution-guide.md

# Contribution Guide

## Getting Started

Thank you for your interest in contributing to `climate_indices`! This guide will help you get started with the contribution process.

`CONTRIBUTING.md` at the repository root is the canonical contributor process;
this guide is the site's walkthrough of it. Where the two differ, the root
document wins.

### Prerequisites
1. **Python 3.10+** installed
2. **uv** package manager installed
3. **Git** for version control
4. Familiarity with pytest and type hints

## Development Setup

### 1. Fork and Clone
```bash
# Fork the repository on GitHub, then clone your fork
git clone https://github.com/YOUR_USERNAME/climate_indices.git
cd climate_indices
```

### 2. Set Up Upstream Remote
```bash
git remote add upstream https://github.com/monocongo/climate_indices.git
git fetch upstream
```

### 3. Install Dependencies
```bash
uv sync --group dev
source .venv/bin/activate
```

## Coding Conventions

### Style Guidelines
- **Indentation**: 4 spaces (no tabs)
- **Line length**: 120 characters
- **Whitespace**: After list items, around operators, around method parameters
- **Comments**: Above-line comments only, liberal use encouraged
- **Naming**: `snake_case` for functions/variables, NO `camelCase`
- **Descriptive names**: Avoid `data`, `df2`, `temp` - use meaningful names

### Type Hints
- **Required**: All function arguments and return types MUST be type-hinted
- **Format**:
  ```python
  def compute_index(
      values: np.ndarray,
      scale: int,
      *,  # Keyword-only separator
      distribution: Distribution = Distribution.gamma,
  ) -> np.ndarray:
      """Compute climate index."""
  ```

### Docstrings
- **Format**: Google-style docstrings
- **Required**: All public functions must have docstrings
- **Example**:
  ```python
  def spi(
      values: np.ndarray,
      scale: int,
      distribution: Distribution,
  ) -> np.ndarray:
      """Compute Standardized Precipitation Index.

      Args:
          values: Precipitation time series in mm
          scale: Temporal scale in months (1-72) or days (1-2196)
          distribution: Fitting distribution (gamma or pearson)

      Returns:
          SPI values as numpy array (same shape as input)

      Raises:
          InsufficientDataError: If fewer than 4 non-zero values for Pearson Type III fitting
          InvalidArgumentError: If scale is out of range

      Example:
          >>> precip = np.array([...])
          >>> spi_6 = spi(precip, scale=6, distribution=Distribution.gamma)
      """
  ```

### Import Organization
```python
# 1. Standard library
from __future__ import annotations
import os
from typing import Optional

# 2. Third-party packages
import numpy as np
import xarray as xr

# 3. Local modules
from climate_indices import compute, exceptions
```

**Forbidden**: Wildcard imports (`from module import *`)

## Testing Requirements

### Test Coverage
- **Minimum**: >90% line coverage
- **New features**: MUST include tests
- **Bug fixes**: MUST include regression test

### Test Structure
```python
# tests/test_new_feature.py
import pytest
import numpy as np
from climate_indices import new_feature

def test_new_feature_basic():
    """Test basic functionality."""
    result = new_feature(input_data)
    assert result.shape == expected_shape
    np.testing.assert_allclose(result, expected_output)

def test_new_feature_edge_case():
    """Test edge case handling."""
    with pytest.raises(InvalidArgumentError):
        new_feature(invalid_input)
```

### Running Tests Before Submission
```bash
# All tests
uv run pytest

# With coverage
uv run pytest --cov=src --cov-report=term

# Code quality
ruff check --fix src/ tests/
ruff format src/ tests/
mypy src/
```

## Contribution Workflow

### 1. Create Feature Branch
```bash
git switch main
git pull --ff-only upstream main
git switch -c fix/bug-description
# Use descriptive branch names:
# - feature/add-new-index
# - fix/handle-nan-values
# - docs/improve-api-examples
```

### 2. Make Changes
- Follow coding conventions above
- Write tests for new functionality
- Update documentation if needed
- Keep commits focused and atomic

### 3. Commit Messages
```bash
# Good: Clear, concise, imperative mood
git commit -m "Add validation for empty latitude arrays"
git commit -m "Fix NaN handling in temporal scaling"

# Bad: Vague, past tense
git commit -m "Fixed stuff"
git commit -m "Changed things"
```

**Format for larger commits**:
```
Add Hargreaves PET method

- Implement daily PET computation
- Add validation for tmin/tmax inputs
- Include comprehensive test suite
- Update documentation with examples
```

### 4. Pull Upstream Changes Regularly
```bash
git fetch upstream
git merge upstream/main
# Or rebase if preferred:
git rebase upstream/main
```

### 5. Push and Create Pull Request
```bash
git push origin fix/bug-description
```

Then on GitHub:
1. Click "New Pull Request"
2. Provide clear description of changes
3. Reference any related issues (`Fixes #123`)
4. Wait for CI to pass
5. Address review comments

## Pull Request Checklist

Before submitting:
- [ ] Tests pass locally (`uv run pytest`)
- [ ] Code is formatted (`ruff format`)
- [ ] Linting passes (`ruff check --fix`)
- [ ] Type checking passes (`mypy src/`)
- [ ] Documentation is updated
- [ ] Commit messages are clear
- [ ] No unrelated changes included

## Code Review Process

### What Reviewers Check
1. **Functionality**: Does the code work as intended?
2. **Tests**: Are there sufficient tests with good coverage?
3. **Style**: Does it follow project conventions?
4. **Documentation**: Are docstrings and comments clear?
5. **Backward compatibility**: Does it break existing APIs?

### Responding to Reviews
- Address all comments
- Ask questions if feedback is unclear
- Make requested changes in new commits
- Mark conversations as resolved when addressed

## Best Practices

### DO:
✅ Write clear, focused pull requests
✅ Include tests for all new code
✅ Follow the existing code style
✅ Update documentation when needed
✅ Respond promptly to review feedback
✅ Keep commits atomic and well-described

### DON'T:
❌ Mix unrelated changes in one PR
❌ Submit without running tests locally
❌ Ignore linting/formatting errors
❌ Break backward compatibility without discussion
❌ Include commented-out code or debug prints
❌ Use magic numbers without constants

## Specific Contribution Areas

### Adding a New Climate Index
1. Implement core algorithm in `compute.py`
2. Add high-level function to `indices.py` (numpy API)
3. Add wrapper to `xarray_adapter.py` (xarray API)
4. Add tests in `tests/test_compute.py`, `test_indices.py`, `test_xarray_adapter.py`
5. Update documentation in `docs/reference.md` and any relevant `docs/*.md` guides

### Fixing a Bug
1. Write a failing test that reproduces the bug
2. Fix the bug in the appropriate module
3. Verify the test now passes
4. Add regression test to prevent recurrence

### Improving Documentation
1. Update the relevant Markdown page in `docs/` for Sphinx docs (`docs/index.md` is the Sphinx root document whose toctrees include them)
2. Update docstrings in source code
3. Build docs locally: `cd docs && make html`
4. Check for broken links and formatting issues

## Getting Help

- **Questions**: Open a GitHub issue with `question` label
- **Bugs**: Open a GitHub issue with `bug` label
- **Features**: Open a GitHub issue with `enhancement` label
- **Discussion**: Use GitHub Discussions

## Recognition

Contributors are acknowledged in:
- Git commit history
- GitHub contributor graph
- Release notes (for significant contributions)

---

Thank you for contributing to `climate_indices`! Your efforts help advance open science and climate research.

## docs/ai-assisted-development.md

# How this project is developed

`climate_indices` is developed by one maintainer directing several concurrent AI
coding agents. The agents carry the volume; the maintainer owns every
irreversible decision. This page describes that workflow so another project can
copy it, and states where it does not help.

## What stays human

Merge authority, release tags, PyPI publication, and scientific sign-off. These
are policy, not etiquette:

- `AGENTS.md` says agents never merge pull requests, push `main`, or create or
  push release tags, and that no plan, handoff, issue text, or prior session's
  notes can authorize a merge.
- `CONTRIBUTING.md` gives contributors the same rule: open a PR, and the
  maintainer merges it after review and passing CI.
- Releases are maintainer-owned and tag-based from `main`, per
  `docs/release-process.md`.
- Scientific claims are bounded by `VALIDATION.md`, which records a tolerance
  and a known gap for every index. Accepting a result the table does not yet
  support stays a human decision.

## Wayfinding before code

New work is shaped before it is implemented. Large efforts are mapped as a set
of decision tickets and resolved one at a time: `wayfinder` for the map,
`grilling` and `domain-modeling` to settle hard-to-reverse choices and pin down
terminology, `to-issues` to slice the result into independently grabbable
issues, and `triage` to move them through the workflow. The GitHub operations
behind those skills are written down in the repository's maintainer notes,
`docs/agent/issue-tracker.md`, which are not published here.

The output is an ordinary GitHub milestone, project board, and numbered issues.
Labels carry the state: `ready-for-agent` opens an issue to parallel agent
sessions, `ready-for-human` keeps it with the maintainer, and `status:*` marks
progress. Large subsystems get an epic — the fire index family shipped as
FIRE-01 through FIRE-21. Branches carry the type prefix and, for issue-backed
work, the issue number, such as `perf/923-vectorize-apply-ufunc`; the prefix
list is in `CONTRIBUTING.md`.

## Parallel execution

Implementation runs in the `herdr` harness, which keeps 8–12 agent sessions
alive at once. Sessions run under the `pi` coding agent, mostly driving DeepSeek
models, with Claude and GPT models used for particular tasks. The tooling is the
maintainer's own choice; nothing in the package depends on it.

The tooling matters less than the isolation rule: every session gets its own
`git worktree` off `origin/main`, documented in `AGENTS.md`. Sessions that share
a checkout overwrite each other's working tree and stage each other's hunks, so
the rule is not optional. A session that must hand off writes a handover note
rather than guessing from history.

## Verification as the design constraint

The working rule is one logical change per pull request, reviewed the same day.
Commit messages conventionally follow Conventional Commits,
`type(scope): summary`. Review findings come back as `fix(review):` commits
rather than silent edits, so the correction stays visible in history; 13 such
commits landed between `v2.4.0` and this document.

Evidence is a merge condition, not a description:

- Performance claims ship a benchmark script and a committed result file.
- Scientific claims ship a reference fixture with a stated tolerance and a
  stated gap; `VALIDATION.md` is the index of those.
- Architectural claims resting on a non-obvious, hard-to-reverse decision ship
  an ADR under `docs/adr/`.

The local gates are the ones `AGENTS.md` lists under *Validate source or test
changes*, and they are what CI runs: lint, format, types, and tests for every
change, the published-docs build with warnings as errors and its doctests for
documentation, and `tests/test_release_integrity.py` for anything touching
packaging or release. The repository also bans AI/tool attribution in commits
and PR descriptions (`CONTRIBUTING.md`): authorship belongs to the human author
only.

## What the workflow produced

Pull requests are the honest unit of output; commit counts double-count merges.
The window below runs from the `v2.4.0` tag (2026-04-06) to commit `3bf9cafc`
(2026-09-17). Every number is recomputed by the commands that follow it, which
use the same window.

| Metric | Value |
| --- | --- |
| Merged pull requests (trunk) | 191 |
| Issues closed | 117 |
| Median PR open-to-merge | about 68 minutes |
| Source changes (`src/`) | +11,723 / −3,881 lines |
| Test changes (`tests/`) | +29,991 / −4,744 lines |
| Documentation changes (`docs/`) | +3,498 / −3,297 lines |
| ADRs | 11 |
| Tests | 2,566 collected; 2,143 in the default gate |

```bash
R=v2.4.0
E=3bf9cafc

git log $R..$E --first-parent --merges --format='%s' \
  | grep -oE '^Merge pull request #[0-9]+' | sort -u | wc -l
gh api -X GET search/issues \
  -f q='repo:monocongo/climate_indices is:issue is:closed closed:2026-04-06..2026-09-17' --jq .total_count
gh pr list --state merged --limit 1000 --json createdAt,mergedAt \
  | jq '[.[] | select(.mergedAt >= "2026-04-06" and .mergedAt <= "2026-09-17T23:59:59Z") | (((.mergedAt|fromdateiso8601) - (.createdAt|fromdateiso8601))/60)] | sort | (if length % 2 == 1 then .[length/2|floor] else ((.[length/2-1] + .[length/2]) / 2) end)'

git diff --shortstat $R $E -- src
git diff --shortstat $R $E -- tests
git diff --shortstat $R $E -- docs
git log $R..$E --format='%s' | grep -c '^fix(review)'

# The ADR and test counts are snapshots at the window end, not at HEAD.
git checkout --detach $E
ls docs/adr/*.md | wc -l
uv run pytest --collect-only -q -o addopts=   # 2,566 collected
uv run pytest --collect-only -q               # 2,143 in the default gate
```

## Limits

Merge, tag, and scientific sign-off stay human by design: a maintainer
authorizes the release, pushes the tag, and approves the release environment.
After those gates, the release workflow publishes the built artifacts to PyPI
and creates the GitHub Release, so the pipeline carries out a release but never
decides one. Volume can outpace
verification, which is why `VALIDATION.md` records gaps instead of a green
checkmark and why the release is gated on a human approving an environment
rather than on a green pipeline. The likeliest failure mode is not wrong code
shipping silently; it is review findings arriving faster than they can be
triaged. This page describes how the chosen work gets verified and merged, not
how to choose what to build.

## docs/xarray_migration.md

# xarray Migration Guide

```{contents}
:backlinks: none
:local: true
```

## Introduction

Starting with version 2.3.0, the `climate_indices` library provides native support for xarray DataArrays alongside the traditional NumPy array API. This migration guide helps existing users understand the benefits of adopting the xarray API and provides practical examples for transitioning from NumPy-based workflows to xarray-based workflows.

:::{note}
**Stable index API from 3.1.0**

The DataArray index APIs in this guide are **Stable** from 3.1.0
([ADR-0012](docs/adr/0012-xarray-api-stays-beta-through-3.0.0.md)).
Parameter inference, metadata attributes, and coordinate handling have the same
no-breaking-changes-in-minor-versions guarantee as the NumPy API.
The promotion changes no numerical algorithm, calendar conversion, or Dask chunk requirement.
It includes index and diagnostic functions that accept DataArrays and return Datasets.

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`.
The warning remains publicly importable from `climate_indices` and
`climate_indices.exceptions`; existing filters need no changes.
Stable DataArray index APIs do not emit it.
No warning import or filter migration is required.
See the [API change record](docs/deprecations/api-changes.md#dataarray-index-api-promotion-310)
for the promotion details and the
[compatibility matrix](docs/xarray_compatibility.md) for each route's constraints.
:::

:::{warning}
**Calendar-alignment fix changes daily results (3.0.0)**

Before 3.0.0, the xarray adapter passed daily Gregorian values straight into
the NumPy core without accounting for its 366-day-per-year layout, silently
shifting every value after February 28 in a non-leap year. This is now
corrected: {func}`climate_indices.spi`, {func}`climate_indices.spei`,
{func}`climate_indices.eddi`, {func}`climate_indices.percentage_of_normal`, and
{func}`climate_indices.xarray_adapter.pet_hargreaves` may return **different,
corrected** daily values for any series spanning a non-leap year. For inputs
using a supported Gregorian calendar (`standard`, `gregorian`, or
`proleptic_gregorian`) with daily coordinates that begin on January 1,
this happens with **no error raised** — the input still succeeds, just with
different numbers. Other previously-successful inputs (unsupported
calendars, or daily coordinates that do not begin on January 1) now raise a
validation error instead of being silently misinterpreted. The NumPy array
API is unaffected. See
`docs/adr/0004-xarray-calendar-semantics.md` and the "Unsupported calendar or
calendar origin" pitfall below for the new, separately-raised validation
errors this release also adds.
:::

**Who this guide is for:** Users currently working with NumPy arrays who want to leverage xarray's labeled dimensions, automatic metadata handling, and coordinate-aware operations.

**What changed:** All primary index functions (`spi()`, `spei()`, `eddi()`, `percentage_of_normal()`, `pet_thornthwaite()`, `pet_hargreaves()`, `pdsi()`) and the audit surface `fit_diagnostics()` now accept both `np.ndarray` and `xr.DataArray` inputs. When you pass an xarray DataArray, the library automatically:

- Infers temporal parameters from time coordinates
- Preserves all coordinate information in outputs
- Adds CF-compliant metadata attributes
- Aligns multi-input arrays automatically
- Supports Dask for out-of-core computation

## Why Migrate to xarray?

The xarray API offers several key advantages over the traditional NumPy API:

**Automatic parameter inference**

: Parameters like `data_start_year`, `periodicity`, and calibration bounds are automatically inferred from the time coordinate, reducing boilerplate and eliminating a common source of errors.

**Coordinate preservation**

: Output DataArrays retain all input coordinates (time, lat, lon, etc.) and their attributes, making it easy to track spatial and temporal metadata through your analysis pipeline.

**CF-compliant metadata**

: Outputs automatically include standardized attributes (`long_name`, `units`, `references`) and a provenance history trail following CF Convention standards, improving reproducibility and dataset documentation.

**Multi-input alignment**

: Functions like `spei()` automatically align precipitation and PET arrays using `xr.align(join="inner")`, handling mismatched time ranges transparently and warning when data is dropped.

**Dask support**

: Pass chunked DataArrays to enable lazy, out-of-core computation on datasets larger than available RAM, with automatic parallelization across spatial dimensions.

**Type safety**

: IDE autocomplete and type checkers like mypy can narrow return types based on input types, catching potential errors at development time rather than runtime.

## Side-by-Side Comparison

This section shows how to convert common NumPy workflows to the equivalent xarray patterns.

### 2.1 SPI — Standardized Precipitation Index

**Before (NumPy):**

```python
import numpy as np
from climate_indices import indices
from climate_indices.compute import Periodicity
from climate_indices.indices import Distribution

# 40 years × 12 months/year = 480 monthly values
precip_mm = np.load("monthly_precip.npy")  # shape: (40, 12)

# compute SPI-6 with gamma distribution
spi_values = indices.spi(
    values=precip_mm,
    scale=6,
    distribution=Distribution.gamma,
    data_start_year=1980,                    # required
    calibration_year_initial=1981,           # required
    calibration_year_final=2010,             # required
    periodicity=Periodicity.monthly,         # required
)
# result: np.ndarray, shape (40, 12)
```

**After (xarray):**

```python
import pandas as pd
import xarray as xr
from climate_indices import spi
from climate_indices.indices import Distribution

# load data as xarray DataArray with time coordinate
time = pd.date_range("1980-01-01", periods=480, freq="MS")
precip_da = xr.DataArray(
    precip_mm.flatten(),
    coords={"time": time},
    dims=["time"],
    attrs={"units": "mm", "long_name": "Monthly Precipitation"},
)

# compute SPI-6 with gamma distribution
spi_result = spi(
    values=precip_da,
    scale=6,
    distribution=Distribution.gamma,
    # data_start_year, calibration bounds, and periodicity
    # are all inferred from the time coordinate!
)
# result: xr.DataArray with time coord and CF metadata
```

**Key changes:**

- Import from the package root (`from climate_indices import spi`) instead of `climate_indices.indices`
- Data shape changes from 2D (years, months) to 1D with labeled time dimension
- Only 2 required parameters (`scale`, `distribution`) instead of 6
- All temporal parameters (`data_start_year`, `calibration_year_initial`, `calibration_year_final`, `periodicity`) are optional and inferred automatically
- Output includes coordinate metadata and CF attributes

### 2.2 SPEI — Standardized Precipitation Evapotranspiration Index

**Before (NumPy):**

```python
import numpy as np
from climate_indices import indices
from climate_indices.compute import Periodicity
from climate_indices.indices import Distribution

# load precipitation and PET arrays
precip_mm = np.load("monthly_precip.npy")  # shape: (40, 12)
pet_mm = np.load("monthly_pet.npy")        # shape: (40, 12)

# arrays must be pre-aligned manually!
assert precip_mm.shape == pet_mm.shape

# compute SPEI-3
spei_values = indices.spei(
    precips_mm=precip_mm,
    pet_mm=pet_mm,
    scale=3,
    distribution=Distribution.gamma,
    periodicity=Periodicity.monthly,
    data_start_year=1980,
    calibration_year_initial=1981,
    calibration_year_final=2010,
)
```

**After (xarray):**

```python
import pandas as pd
import xarray as xr
from climate_indices import spei
from climate_indices.indices import Distribution

# load precipitation and PET as DataArrays
# they can have different time ranges!
time_precip = pd.date_range("1980-01-01", periods=480, freq="MS")
time_pet = pd.date_range("1985-01-01", periods=420, freq="MS")

precip_da = xr.DataArray(
    precip_values,
    coords={"time": time_precip},
    dims=["time"],
)
pet_da = xr.DataArray(
    pet_values,
    coords={"time": time_pet},
    dims=["time"],
)

# compute SPEI-3 — automatic inner join alignment!
spei_result = spei(
    precips_mm=precip_da,
    pet_mm=pet_da,
    scale=3,
    distribution=Distribution.gamma,
    # alignment, calibration, and periodicity all handled automatically
)
# result covers 1985-2019 (intersection of input time ranges)
```

**Key changes:**

- Import from `climate_indices.spei` instead of `climate_indices.indices.spei`
- Inputs no longer need to be pre-aligned — `xr.align(join="inner")` handles mismatched time ranges
- A warning is emitted if alignment drops timesteps, helping detect data issues
- Only 3 required parameters (`scale`, `distribution`) instead of 7

### 2.3 PET Thornthwaite

**Before (NumPy):**

```python
from climate_indices import indices

# 1D array of 480 monthly temperatures
temp_celsius = np.load("monthly_temp.npy")  # shape: (480,)

# compute PET using Thornthwaite method
pet_mm = indices.pet(
    temperature_celsius=temp_celsius,
    latitude_degrees=40.0,
    data_start_year=1980,  # required
)
# result: np.ndarray, shape (480,)
```

**After (xarray):**

```python
from climate_indices.xarray_adapter import pet_thornthwaite

# load as DataArray
time = pd.date_range("1980-01-01", periods=480, freq="MS")
temp_da = xr.DataArray(
    temp_values,
    coords={"time": time},
    dims=["time"],
    attrs={"units": "degC"},
)

# compute PET — data_start_year inferred from time coord
pet_result = pet_thornthwaite(
    temperature=temp_da,
    latitude=40.0,
)
# result: xr.DataArray with CF metadata
```

**Gridded example with spatial broadcasting:**

```python
# load gridded temperature (time, lat, lon)
temp_grid = xr.DataArray(
    temp_3d_values,  # shape: (480, 5, 6)
    coords={
        "time": pd.date_range("1980-01-01", periods=480, freq="MS"),
        "lat": [30, 35, 40, 45, 50],
        "lon": [-120, -110, -100, -90, -80, -70],
    },
    dims=["time", "lat", "lon"],
)

# latitude as DataArray for spatial broadcasting
lat_array = xr.DataArray([30, 35, 40, 45, 50], dims=["lat"])

# automatically broadcasts latitude across grid
pet_grid = pet_thornthwaite(temp_grid, latitude=lat_array)
# result: xr.DataArray, shape (480, 5, 6)
```

**Key changes:**

- Import from `climate_indices.xarray_adapter` instead of `climate_indices.indices`
- `data_start_year` is optional and inferred from time coordinate
- Latitude can be a scalar or DataArray for spatial broadcasting
- Function name changes to `pet_thornthwaite()` for clarity

### 2.4 PET Hargreaves

**Before (NumPy):**

```python
from climate_indices.eto import eto_hargreaves

# 1D arrays of daily min/max temperatures
tmin = np.load("daily_tmin.npy")  # shape: (1825,)
tmax = np.load("daily_tmax.npy")  # shape: (1825,)

# manually derive mean temperature
tmean = (tmin + tmax) / 2.0

# compute PET using Hargreaves method
pet_mm = eto_hargreaves(
    daily_tmin_celsius=tmin,
    daily_tmax_celsius=tmax,
    daily_tmean_celsius=tmean,  # must provide manually
    latitude_degrees=40.0,
)
```

**After (xarray):**

```python
from climate_indices.xarray_adapter import pet_hargreaves

# load as DataArrays
time = pd.date_range("2015-01-01", periods=1825, freq="D")
tmin_da = xr.DataArray(tmin_values, coords={"time": time}, dims=["time"])
tmax_da = xr.DataArray(tmax_values, coords={"time": time}, dims=["time"])

# compute PET — tmean derived automatically!
pet_result = pet_hargreaves(
    daily_tmin_celsius=tmin_da,
    daily_tmax_celsius=tmax_da,
    latitude=40.0,
)
# result: xr.DataArray with CF metadata
```

**Key changes:**

- Import from `climate_indices.xarray_adapter`
- Mean temperature is automatically derived as `(tmin + tmax) / 2`
- Inputs are automatically aligned via `xr.align(join="inner")`
- Supports spatial broadcasting with latitude as DataArray

## Understanding the Output Metadata

When you use the xarray API, outputs include rich metadata that improves reproducibility and interoperability with other tools in the scientific Python ecosystem.

**Example SPI output attributes:**

```python
result = spi(precip_da, scale=6, distribution=Distribution.gamma)
print(result.attrs)

# {
#     'long_name': 'Standardized Precipitation Index',
#     'units': 'dimensionless',
#     'references': 'McKee, T. B., Doesken, N. J., & Kleist, J. (1993). ...',
#     'scale': 6,
#     'distribution': 'gamma',
#     'calibration_year_initial': 1980,
#     'calibration_year_final': 2019,
#     'climate_indices_version': '2.3.0',
#     'history': '2026-02-09T14:23:45Z: SPI-6 calculated using gamma distribution (climate_indices v2.3.0)',
# }
```

### CF Convention Compliance

The library follows CF (Climate and Forecast) Convention standards for metadata:

**CF metadata fields:**
: - `long_name`: Human-readable description of the variable
  - `units`: Physical units (e.g., "dimensionless", "mm/month")
  - `references`: Academic citation for the algorithm
  - `standard_name`: CF standard name (when officially defined)

**Calculation metadata:**

: Parameters used in the computation (`scale`, `distribution`, `calibration_year_initial`, `calibration_year_final`) are stored as attributes for reproducibility.

**Provenance tracking:**

: The `history` attribute follows CF Convention format with ISO 8601 timestamps. When you chain operations, new history entries are appended with newline separators:

  ```python
  # original history from upstream processing
  precip_da.attrs["history"] = "2026-02-08T10:00:00Z: Data downloaded from NCAR"

  # compute SPI
  result = spi(precip_da, scale=3, distribution=Distribution.gamma)

  # history is appended, not replaced
  print(result.attrs["history"])
  # 2026-02-08T10:00:00Z: Data downloaded from NCAR
  # 2026-02-09T14:23:45Z: SPI-3 calculated using gamma distribution (climate_indices v2.3.0)
  ```

### Coordinate Preservation

All coordinates and their attributes are deep-copied from input to output, preserving spatial and temporal metadata:

**Dimension coordinates:**

: Primary axis coordinates (time, lat, lon) with their attributes

**Auxiliary coordinates:**

: Non-dimension coordinates like month-of-year or station ID

**Scalar coordinates:**

: Zero-dimensional coordinates like ensemble member or experiment name

```python
# input with rich coordinate metadata
precip_da.coords["time"].attrs = {
    "axis": "T",
    "calendar": "standard",
    "long_name": "time",
}

# compute SPI
result = spi(precip_da, scale=6, distribution=Distribution.gamma)

# all coordinate attributes are preserved
assert result.coords["time"].attrs == precip_da.coords["time"].attrs
```

### Direct NetCDF export

Because outputs are standard xarray DataArrays with CF-compliant metadata, you can write them directly to NetCDF files without additional metadata preparation:

```python
# compute index
spi_result = spi(precip_da, scale=12, distribution=Distribution.gamma)

# write to NetCDF with full metadata
spi_result.to_netcdf("spi_12month_1980_2019.nc")

# later, reload with metadata intact
spi_loaded = xr.open_dataarray("spi_12month_1980_2019.nc")
print(spi_loaded.attrs["references"])  # citation preserved
```

## Common Pitfalls and Solutions

### 1. Time dimension must be named "time"

:::{warning}
**Error:** `DimensionMismatchError: Time dimension 'time' not found in input. Available dimensions: ['date']. Use time_dim parameter to specify custom name.`

**Cause:** Your DataArray uses a different dimension name (e.g., `date`, `timestamp`, `t`).

**Solution:** Rename the dimension to `time` before calling the function:

```python
# if your data uses "date" instead of "time"
da_renamed = da.rename({"date": "time"})
result = spi(da_renamed, scale=6, distribution=Distribution.gamma)
```
:::

### 2. Non-monotonic time coordinate

:::{warning}
**Error:** `CoordinateValidationError: Time coordinate is not monotonically increasing. Sort the data using data.sortby('time') before processing.`

**Cause:** Time values are not in chronological order (e.g., data was shuffled or merged from multiple sources).

**Solution:** Sort by time coordinate before processing:

```python
da_sorted = da.sortby("time")
result = spi(da_sorted, scale=6, distribution=Distribution.gamma)
```
:::

### 3. Unsupported frequency

:::{warning}
**Error:** `CoordinateValidationError: Unsupported frequency 'W' - only 'MS'/'ME'/'M' (monthly) and 'D' (daily) supported`

**Cause:** Your data has weekly, quarterly, or annual frequency.

**Solution:** The library only supports monthly (MS/ME/M) and daily (D) data. Resample your data to a supported frequency:

```python
# resample weekly to monthly using sum aggregation
da_monthly = da.resample(time="MS").sum()
result = spi(da_monthly, scale=6, distribution=Distribution.gamma)
```
:::

### 4. Unsupported calendar or calendar origin

:::{warning}
**Error:** `CoordinateValidationError: Unsupported calendar` or a message
that monthly input must begin in January or daily input must begin on January 1.

**Cause:** SPI, SPEI, EDDI, and PNP group monthly values from January and
daily values into 366 calendar positions. `pet_thornthwaite` (monthly) and
`pet_hargreaves` (daily) enforce the same contract. `cftime` calendars and
non-January input origins cannot be represented without changing those
semantics.

**Solution:** For those index computations, use standard, gregorian, or
proleptic_gregorian datetime coordinates. Start monthly data in January and
daily data on January 1. A final daily year may end on any calendar date.
For non-leap years, the xarray adapter synthesizes February 29 from February
28 and March 1 before computation, then restores the original Gregorian
coordinate.
:::

### 5. Multi-chunked time dimension (Dask)

:::{warning}
**Error:** `CoordinateValidationError: Dimension 'time' is split across 4 chunks. Climate index computation requires this dimension in a single chunk. Rechunk using: data = data.chunk({'time': -1})`

**Cause:** Your Dask-backed DataArray has the time dimension split into multiple chunks. SPI/SPEI require the full time series for fitting distributions.

**Solution:** Rechunk time dimension to a single chunk:

```python
# rechunk time to single chunk, keep spatial chunks
da_rechunked = da.chunk({"time": -1, "lat": 10, "lon": 10})
result = spi(da_rechunked, scale=6, distribution=Distribution.gamma)
```
:::

### 6. Passing Dataset instead of DataArray

:::{warning}
**Error:** `InputTypeError: Unsupported input type: xarray.core.dataset.Dataset. Accepted types: np.ndarray, list, tuple, int, float, np.integer, np.floating, xr.DataArray. xr.Dataset detected: Use ds['variable_name'] to select a DataArray`

**Cause:** You passed an `xr.Dataset` (collection of variables) instead of an `xr.DataArray` (single variable).

**Solution:** Select a specific variable from the Dataset:

```python
# if you loaded a NetCDF file with xr.open_dataset()
ds = xr.open_dataset("precipitation.nc")

# extract the precipitation variable as a DataArray
precip_da = ds["precipitation"]
result = spi(precip_da, scale=6, distribution=Distribution.gamma)
```
:::

### 7. Insufficient calibration data with NaN

:::{warning}
**Error:** `InsufficientDataError: Insufficient non-NaN data in calibration period (1980-2010). Found 245 non-NaN values (20.4 effective years), but at least 30 years of non-NaN data required for reliable distribution fitting.`

**Cause:** Your data contains too many NaN values within the calibration period, leaving fewer than 30 years of valid observations. This check applies to a single-series (1-D) in-memory input; gridded in-memory and Dask-backed inputs are not preflighted.

**Solution:** Either extend the calibration period or use a different time range with denser data:

```python
# option 1: extend calibration period if you have more data
result = spi(
    precip_da,
    scale=6,
    distribution=Distribution.gamma,
    calibration_year_initial=1950,  # earlier start
    calibration_year_final=2020,    # later end
)

# option 2: use only the dense portion of your dataset
dense_period = precip_da.sel(time=slice("1990", "2020"))
result = spi(dense_period, scale=6, distribution=Distribution.gamma)
```
:::

## Performance Considerations

### When to use Dask chunking

Use Dask-backed DataArrays when:

- **Gridded datasets exceed available RAM** (e.g., global climate model output)
- **Spatial dimensions are large** (e.g., 720×360 global grid)
- **Processing multiple realizations** (ensemble forecasts with time×lat×lon×member dimensions)

Do NOT use Dask for:

- **1D time series** — overhead exceeds benefit for small arrays
- **Small grids** (< 100 spatial points) — eager computation is faster

```python
# good use case: large gridded dataset
precip = xr.open_dataarray(
    "global_precip.nc",
    chunks={"time": -1, "lat": 50, "lon": 50},  # single time chunk
)
result = spi(precip, scale=6, distribution=Distribution.gamma)
result.to_netcdf("spi_global.nc")  # triggers computation
```

### Optimal chunking strategy

For climate index calculations:

1. **Time dimension:** Must be a single chunk (`{"time": -1}`)

2. **Spatial dimensions:** Chunk to balance parallelism vs. overhead

   - Too large: Limited parallelism, potential memory issues
   - Too small: Excessive overhead from task scheduling

   A good starting point: **10-20 chunks per spatial dimension**

```python
# example: 40 years × 180 lat × 360 lon global grid
da = xr.open_dataarray(
    "global_data.nc",
    chunks={
        "time": -1,      # single chunk (required)
        "lat": 18,       # 10 chunks (180/18)
        "lon": 36,       # 10 chunks (360/36)
    },
)
```

### Overhead

**1D time series:** Xarray overhead is minimal (\<5%) compared to NumPy path. The convenience and metadata benefits outweigh the small performance cost.

**Gridded data:** Performance is computation-dominated. Overhead from xarray/Dask is negligible relative to distribution fitting time. Spatial parallelism via Dask provides significant speedup on multi-core systems.

**Benchmark results** (40-year monthly data, SPI-6):

| Input type | Shape         | Time                 |
| ---------- | ------------- | -------------------- |
| NumPy      | (40, 12)      | 12.3 ms              |
| xarray 1D  | (480,)        | 12.8 ms (+4%)        |
| xarray 3D  | (480, 20, 20) | 1.42 s               |
| Dask 3D    | (480, 20, 20) | 0.58 s (2.4× faster) |

### Lazy evaluation

Dask-backed DataArrays use lazy evaluation — computations build a task graph but don't execute until you call `.compute()` or trigger output (e.g., `.to_netcdf()`). Invalid SPI/SPEI options (including timescale, distribution, and output scale) raise when you call the function, before a graph is built:

```python
# open dataset with chunking (no data loaded yet)
precip = xr.open_dataarray("large_file.nc", chunks={"time": -1, "lat": 20})

# compute SPI (builds task graph, no computation yet)
spi_result = spi(precip, scale=6, distribution=Distribution.gamma)

# trigger computation and load result into memory
spi_eager = spi_result.compute()

# OR write directly to disk without loading into memory
spi_result.to_netcdf("output.nc")  # triggers compute + streams to disk
```

## Real-World Examples

### Example 1: SPI from NetCDF

```python
import xarray as xr
from climate_indices import spi
from climate_indices.indices import Distribution

# load monthly precipitation from NetCDF
ds = xr.open_dataset("PRISM_precip_monthly_1981_2020.nc")
precip = ds["ppt"]  # extract DataArray

# compute SPI-12 using full time range for calibration
spi_12 = spi(
    values=precip,
    scale=12,
    distribution=Distribution.gamma,
)

# add descriptive metadata
spi_12.attrs["title"] = "SPI-12 for CONUS 1981-2020"
spi_12.attrs["source"] = "PRISM Climate Group"

# write to NetCDF with full provenance
spi_12.to_netcdf("spi_12_CONUS_1981_2020.nc")
```

### Example 2: SPEI from separate files

```python
import xarray as xr
from climate_indices import spei
from climate_indices.indices import Distribution

# load precipitation and PET from different files
# they may have different time ranges!
precip = xr.open_dataarray("precip_1950_2020.nc")
pet = xr.open_dataarray("pet_1980_2020.nc")

# compute SPEI-6 — automatic inner join alignment
spei_6 = spei(
    precips_mm=precip,
    pet_mm=pet,
    scale=6,
    distribution=Distribution.gamma,
)
# result covers 1980-2020 (intersection of input ranges)

# explicit calibration period (optional)
spei_6_custom = spei(
    precips_mm=precip,
    pet_mm=pet,
    scale=6,
    distribution=Distribution.gamma,
    calibration_year_initial=1991,
    calibration_year_final=2020,
)
```

### Example 3: Gridded PET with spatial latitude broadcasting

```python
import pandas as pd
import xarray as xr
from climate_indices.xarray_adapter import pet_thornthwaite

# load gridded monthly temperature
temp_grid = xr.open_dataarray("temp_monthly_grid.nc")
# shape: (480 time, 50 lat, 60 lon)

# create latitude DataArray matching the grid
lat_array = xr.DataArray(
    temp_grid.coords["lat"].values,
    dims=["lat"],
)

# compute PET — automatically broadcasts latitude across grid
pet_grid = pet_thornthwaite(
    temperature=temp_grid,
    latitude=lat_array,
)
# result shape: (480 time, 50 lat, 60 lon)

# write to NetCDF
pet_grid.to_netcdf("pet_thornthwaite_grid.nc")
```

### Example 4: Out-of-core SPI with Dask

```python
import xarray as xr
from climate_indices import spi
from climate_indices.indices import Distribution

# open large global dataset with Dask chunking
# (data not loaded into memory yet)
precip_global = xr.open_dataarray(
    "CMIP6_pr_global_1850_2100.nc",
    chunks={
        "time": -1,     # single time chunk (required for SPI)
        "lat": 30,      # 6 chunks along latitude
        "lon": 60,      # 6 chunks along longitude
    },
)

# compute SPI-12 (builds task graph, no computation yet)
spi_12_global = spi(
    values=precip_global,
    scale=12,
    distribution=Distribution.gamma,
    calibration_year_initial=1981,
    calibration_year_final=2010,
)

# write to NetCDF — triggers Dask computation under a distributed client, so
# the blocks are fitted in worker processes and streamed to disk. Without a
# client, the default threaded scheduler runs the Python-level fitting work
# serially. (Install the client with `pip install distributed`.)
from dask.distributed import Client

with Client(n_workers=4, threads_per_worker=1) as client:
    spi_12_global.to_netcdf("spi_12_global_cmip6.nc")
```

## API Reference Links

**Index functions:**

- {func}`climate_indices.spi` — Standardized Precipitation Index
- {func}`climate_indices.spei` — Standardized Precipitation Evapotranspiration Index
- {func}`climate_indices.xarray_adapter.pet_thornthwaite` — Potential Evapotranspiration (Thornthwaite)
- {func}`climate_indices.xarray_adapter.pet_hargreaves` — Potential Evapotranspiration (Hargreaves)
- {func}`climate_indices.eddi` — Evaporative Demand Drought Index
- {func}`climate_indices.percentage_of_normal` — Percent of Normal Precipitation
- {func}`climate_indices.pdsi` — Palmer Drought Severity Index family (Dataset output)
- {func}`climate_indices.fit_diagnostics` — Distribution fit diagnostics (Dataset output, stable from 3.1.0)

**Supporting classes:**

- {class}`climate_indices.indices.Distribution` — Distribution types for fitting
- {class}`climate_indices.compute.Periodicity` — Time series periodicity (monthly/daily)

**Exceptions:**

- {exc}`climate_indices.exceptions.CoordinateValidationError` — Invalid or missing coordinates
- {exc}`climate_indices.exceptions.InputTypeError` — Unsupported input type
- {exc}`climate_indices.exceptions.InsufficientDataError` — Not enough data for calibration
- {exc}`climate_indices.exceptions.InputAlignmentWarning` — Data dropped during multi-input alignment

**Full API reference:**

- {doc}`reference` — Complete API documentation for all modules

**External documentation:**

- [xarray documentation](https://docs.xarray.dev/)
- [NumPy documentation](https://numpy.org/doc/stable/)
- [Dask documentation](https://docs.dask.org/)
- [CF Conventions](http://cfconventions.org/)

## docs/performance.md

# Gridded Performance

:::{note}
The gridded DataArray index API is stable from 3.1.0, as is the NumPy API.
This page covers the Dask execution path for the index adapters described in
the [xarray compatibility matrix](docs/xarray_compatibility.md).
Compound Dataset routes and eager event catalogs remain beta.
:::

The gridded index kernels are vectorized once per spatial block: a Dask-backed
`(time, lat, lon)` input reaches the NumPy core once per block instead of once
per grid cell. Chunking, memory, and worker processes are then Dask's job, and
the choices that matter are the chunk shape and the scheduler. This page gives a
complete runnable example and the measured numbers behind the recommendations in
[Operational Guidance](docs/xarray_compatibility.md#operational-guidance) and
[Chunking guidance for gridded indices](docs/xarray_compatibility.md#chunking-guidance-for-gridded-indices).

## Runnable example: gridded SPI and SPEI

This computes 3-month SPI and SPEI for every cell of a synthetic monthly grid,
with `time` in one chunk and spatial blocks driving the parallelism. Both
indices are materialized in one `dask.compute` call, so they share the input
blocks and one worker pool.

```python
import dask
import numpy as np
import pandas as pd
import xarray as xr

from climate_indices import spei, spi
from climate_indices.indices import Distribution

def main() -> None:
    time = pd.date_range("1981-01-01", periods=40 * 12, freq="MS")
    lat = np.linspace(25.0, 49.0, 25)
    lon = np.linspace(-125.0, -101.0, 25)
    shape = (time.size, lat.size, lon.size)
    rng = np.random.default_rng(42)

    precip = xr.DataArray(
        rng.gamma(shape=2.0, scale=15.0, size=shape),
        coords={"time": time, "lat": lat, "lon": lon},
        dims=["time", "lat", "lon"],
        attrs={"units": "mm"},
    )
    pet = xr.DataArray(
        np.full(shape, 60.0),
        coords={"time": time, "lat": lat, "lon": lon},
        dims=["time", "lat", "lon"],
        attrs={"units": "mm"},
    )

    # Time in one chunk is required; spatial chunks are the parallelism and memory lever.
    chunks = {"time": -1, "lat": 10, "lon": 10}
    precip = precip.chunk(chunks)
    pet = pet.chunk(chunks)

    # No Python loop over cells: each index runs once per spatial block.
    spi_lazy = spi(
        values=precip,
        scale=3,
        distribution=Distribution.gamma,
        calibration_year_initial=1981,
        calibration_year_final=2010,
    )
    spei_lazy = spei(
        precips_mm=precip,
        pet_mm=pet,
        scale=3,
        distribution=Distribution.gamma,
        calibration_year_initial=1981,
        calibration_year_final=2010,
    )

    # One process pool and one shared input graph for both indices. The
    # `__main__` guard keeps spawned workers from re-running this block.
    spi_grid, spei_grid = dask.compute(spi_lazy, spei_lazy, scheduler="processes")


if __name__ == "__main__":
    main()
```

`spi_grid` and `spei_grid` are `DataArray`s with the input coordinates, one
value per cell and month, computed without a per-cell Python loop. With
synthetic values the gamma goodness-of-fit check can report
`GoodnessOfFitWarning` for a few cell-months; that is the fit report, not a
failure, and under `scheduler="processes"` it is raised in the workers (their
stderr), not catchable in the caller. Give multi-input indices such as SPEI the
same chunking on the dimensions they share with the other input: each input's
`time` chunking is validated independently, and a mismatch on `lat`/`lon`
survives into the compute as a Dask `rechunk-merge` copy.

## Numerical equivalence

The Spatial Kernel runs the same NumPy core as the serial API, so on a fully populated
monthly grid a gridded result is the serial result: SPI and SPEI with the gamma
distribution, EDDI, and percentage of normal are bit-for-bit identical to calling the
NumPy API once per cell. Two paths differ by a few float64 ULP instead, on CPUs where
scalar and broadcast evaluation of the same library call round differently: Thornthwaite
PET, whose Spatial Block form reorders the same arithmetic (1.14e-13 measured on the
`5 x 6` test grid), and the Pearson L-moment fit used by `spi`/`spei` with
`Distribution.pearson` (4.9e-15 on that grid under NumPy 2.4/scipy 1.17, 1.4e-14 under the
minimum dependencies). Both are asserted at `atol=1e-12`.

A Spatial Block runs one Pearson fit call over its cells; when that fit fails, `spi`
falls back to gamma for the whole block rather than per cell. That block semantics is
documented in [ADR-0009](docs/adr/0009-spatial-block-declaration.md) and is not per-cell
equivalence.

Chunk layout does not change the numbers either: the same grid computed from a
different chunk shape matches the in-memory result, bit for bit for gamma SPI, EDDI, and
percentage of normal, and within the PET tolerance above. Because the Pearson
gamma fallback is applied per Spatial Block, a chunk shape that isolates different cells
can change which cells fall back, so the chunk-layout bounds cover the gamma path only.
`tests/test_numerical_equivalence.py` enforces all of these bounds, so divergence fails
`pytest` instead of silently changing an index value. Daily xarray grids are
calendar-adapted to an all-leap 366-day series before the NumPy core and converted back
afterward, so their serial counterpart is that adapted series rather than the raw Gregorian
input. Missing-data, partial-year, and daily-grid shapes keep the looser bounds asserted
in `tests/test_spatial_kernel.py`.

## Measured speedup

Vectorization is what removes the per-cell Python loop, so its effect is
measured against the serial in-memory call. On the #893 reference grid
(38 x 87 cells, 40 years monthly, scale 3, fastest of three runs after a warm-up
with logging quiet so the fit is isolated; the speedups compare quiet-logging
serial runs and are ±15% rather than exact; Python 3.14.7, 10-core macOS arm64):

| index | serial before | serial after | vectorization speedup |
| --- | ---: | ---: | ---: |
| SPI | 0.898 s | 0.204 s | 4.4x |
| SPEI | 0.783 s | 0.220 s | 3.6x |
| Thornthwaite PET | 1.063 s | 0.020 s | 53x |
| EDDI | 14.717 s | 0.043 s | 342x |

The "serial before" column is the last commit before the spatial-block
conversion; the per-cell logging volume disappeared with the loop, collapsing
the INFO-logged figures as well. The raw runs, the full command, and the
end-to-end Dask numbers are in
[the #929 before/after table](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md#beforeafter-on-the-reference-grid-929).

On this grid the Dask path does not beat the serial call for SPI, SPEI, or PET:
after vectorization each finishes in about 0.2 s or less, while a fresh
`processes` pool plus result transfer adds roughly 0.7 s at one worker. EDDI,
whose pre-vectorization call took 14.7 s, is the only index whose end-to-end Dask
figure still clears 10x (0.699 s, ~21x, and at a single worker): the win is the
vectorization surviving the fixed pool cost, not the parallelism. Parallelism
pays when the work per block exceeds that fixed cost -- larger grids, longer
series, or a worker pool that outlives one compute call.

The real-grid case is that boundary measured (#1097): on the Morocco CHIRPS v3
grid (40,260 cells, 528 months, SPI-6 gamma) -- after the numerical-equivalence
gate passed (gamma bit-for-bit against the per-cell NumPy API, and the Dask grid
bit-for-bit against the eager grid under every worker count the benchmark times)
-- `main` runs the eager in-memory
call in 1.22 s while the `v2.4.0` release took 34.35 s, and the release's Dask
path needed 8 workers (11.70 s) to come within 10x of the new eager figure. On
`main` the one-worker Dask call is 1.32 s above the eager call (2.55 s vs
1.22 s), and the best configuration (4 workers) is still 1.4x slower than the
same call in memory. Interpretation: that 1.32 s is the total process-scheduler
overhead, and it dominates because the whole grid is only 1.2 s of serial work.
The numbers, method, and retained raw output are in
[the #1097 section of benchmarks/README.md](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md#real-grid-spi-benchmark-v240-vs-main-on-the-morocco-chirps-case-1097).

#1097 compares eager against Dask -- both xarray/Dask API. #1121 puts the
CLI's `multiprocessing.Pool` path (ADR-0002) in the same comparison, at CONUS
scale: on a full nClimGrid-Monthly precipitation grid (825,460 cells, 469,758
land, 528 months, SPI-6 gamma), the CLI's 9-worker Pool path (28.68 s compute)
is faster than both the xarray eager call (63.16 s) and the best Dask
configuration (8 workers, 65.00 s) -- 2.2x and 2.3x respectively -- with
output equivalence asserted between the CLI and xarray/Dask paths before that
comparison. Dask does not beat eager at this scale either (8 workers, 65.00 s,
is slightly slower than eager's 63.16 s), consistent with #1097's finding that
Dask's process-scheduler overhead dominates once the serial call is fast.
Interpretation: the CLI's `multiprocessing.Array` shared memory is a
zero-copy buffer every worker indexes into directly, while Dask's
`scheduler="processes"` pickles each block's input and result across process
boundaries -- a cost that appears to outweigh Dask's parallelism gain at this
grid's land-cell count. The numbers, method, and retained raw output are in
[the #1121 section of benchmarks/README.md](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md#legacy-cli-multiprocessing-vs-xarraydask-on-conus-nclimgrid-1121).

## Scheduler

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)):

```python
spi_grid = spi_lazy.compute(scheduler="processes")
```

- Use `scheduler="processes"` for the CPU-bound fitting and ranking kernels: the
  default threaded scheduler does not run their Python-level portion in parallel.
- Every `.compute(scheduler="processes")` builds and tears down its own pool, so
  short or repeated computations spend most of their wall clock there. Keep a
  `dask.distributed.Client` alive across calls (`distributed` ships with the
  `dev` extra), pass a pre-created pool as
  `dask.compute(..., scheduler="processes", pool=pool)`, or compute several lazy
  results in one `dask.compute` call as in the example above.
- Dask's processes scheduler batches up to six ready tasks per submission by
  default, which can flatten scaling on large grids; pass `chunksize=1` to
  `dask.compute` (or set `dask.config.set({"chunksize": 1})`) so each ready
  block is submitted as soon as a worker is free.
- The threaded scheduler remains the right choice for I/O-bound work, reductions
  such as `.mean("time")`, and small in-memory grids where serializing each block
  to a worker costs more than the parallelism saves.

The full scheduler guidance, including Zarr and NetCDF write behavior, is in
[Operational Guidance](docs/xarray_compatibility.md#operational-guidance).

## Running on a cluster

A `dask.distributed` client keeps the worker pool across calls: it outlives one
`.compute()`, can stream a Zarr write, and lets a later diagnostic reuse persisted
blocks. Size it from the block budget rather than from the grid — one worker
process per CPU with one thread each, since the kernels are Python/scipy work, and
memory per worker for the blocks it runs at once (the measured per-block figure is
in [Chunking](#chunking)) — and leave a few blocks per worker.

The [notebook's client cell](https://github.com/monocongo/climate_indices/blob/main/notebooks/zarr_dask_spi_spei.ipynb)
creates one from the host's physical memory, uses it as the default scheduler for
the calculation and the Zarr write, and the write cell closes it in a `finally`
block on both the success and failure paths. The measured fresh-pool versus
warm-pool cost is in
[Operational Guidance](docs/xarray_compatibility.md#operational-guidance).

## Chunking

Keep `time` in a single chunk and size spatial blocks by cell count, not shape.
The PET adapters (`pet_thornthwaite`, `pet_hargreaves`) are the exception: they
accept a split `time` and rechunk it internally, paying that copy inside every
call. The measured working set of a monthly gamma fit is about 60 KB per cell,
so a 100 MB per-block budget is roughly 1,600 cells. 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 reference grid measures well under that at `10 x 10` to
`20 x 20`. Daily grids carry up to 366 steps per cell-year instead of 12 and
want blocks near `7 x 7`. Chunk both spatial dimensions rather than long rows,
and leave at least a few blocks per worker.

Rechunk once, at read or prepare time: `.chunk(...)` only sets the layout of a
lazy graph, so persist the rechunked array (or write the layout back to the
store) before the index calls to pay the copy once. The measured block table and
the daily-grid caveat are in
[Chunking guidance for gridded indices](docs/xarray_compatibility.md#chunking-guidance-for-gridded-indices).

## Fit diagnostics

`indices.fit_diagnostics()` returns the fitted parameters, the probability of zero,
the valid calibration count, and the Kolmogorov-Smirnov D statistic and exact
p-value for every calendar step with a valid sample and usable fitted parameters.
It is for auditing a fit, not for operational grids: the index path skips the
exact p-value whenever a fit is clearly acceptable, while diagnostics computes it
for each such step and cell, and each p-value is a `scipy.stats.kstest` call. On a
gridded block that is the dominant cost by a wide margin. Run it on a calibration
series, or on the cells you are auditing, rather than on a full production grid.
The same NumPy core backs the xarray `fit_diagnostics()` surface, which returns the
diagnostics as an `xr.Dataset` over a `month`/`dayofyear` dimension plus the input's
cell dimensions; it re-fits per Dask spatial chunk, so the same block cost applies.

## Measuring your own

Re-run the reference harness on your own machine and dependency set rather than
trusting these numbers. The scaling harness covers SPI, SPEI, PET, and EDDI from
one to N workers:

```bash
uv run benchmarks/parallel_scaling.py --indices spi,spei,pet,eddi --repeat 3
```

`benchmarks/profile_gridded_spi.py` is the `cProfile` entry point for
attributing a slow SPI pass, and
[`benchmarks/README.md`](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md)
documents both harnesses and the committed raw output.

## Rust backend

The optional Rust backend (`climate_indices._native`, built from `crates/`) is
transparent to this page: the xarray adapters call the same Python entry points, and
dispatch happens inside those modules, so the adapter's chunking and scheduler choices
are unchanged. Native kernels can copy inputs across the PyO3 boundary, adding per-block
allocations and peak worker memory; account for those copies when sizing chunks.

RUST-011 measured every ported kernel against the Python implementation it replaces,
from one registry rather than one family at a time. Two results bear on the guidance
above:

- The per-cell kernels are 34x to 219x faster (Thornthwaite, the Palmer recursions, the
  fire recurrences, Hargreaves, the Antecedent Precipitation Index), and the fitting-based
  indices gain 1.1x to 1.8x: their fit and transform are ported and dispatch natively on
  the Rust path, but the Python path they replace already calls SciPy's compiled
  routines, so little time moves. Importing the package costs about 0.7 seconds in this
  fresh-interpreter measurement, the extension included.
- The dispatched kernels release the GIL, so four threads across four 1024-cell SPI blocks finish
  in 0.096 s against 0.278 s single-threaded with no Rust-side parallelism. Rayon is
  therefore not adopted; an outer pool already parallelizes cell blocks.
- At CONUS scale (an opt-in run of the same harness on nClimGrid-Monthly, 469,758
  land cells x 528 months), one eager call is 1.33x to 1.70x faster through Rust for
  SPI, SPEI, and EDDI and 2.8x for PDSI, but most of the SPI and SPEI gain is the
  Python path's cost on whole-grid arrays: in blocks of about 58,700 cells on one
  thread they are at parity (0.98 to 1.07), and Thornthwaite is at parity to a slight Rust
  lead for eager, one, and two threads (1.03 to 1.12) and 1.17-1.68x ahead at four and
  eight, because its spatial Python path is already vectorized. PDSI (6.2x to 9.4x in
  blocks), EDDI (1.7x to 2.0x), and SPI Pearson (1.1x to 2.0x) keep or gain, and the
  fitting-based indices pull ahead at eight threads (1.5x to 2.0x). Time inside the
  extension -- binding, copy-in, and kernel, not kernel alone -- is most of only the
  Thornthwaite (97%), EDDI (80%), and SPI gamma/PDSI (58%) calls, and an outer pool splits
  those, so the Rayon decision stands. Either backend computes a full grid faster in
  chunks than in one call, except Rust EDDI, where eager and chunked tie (3.93 s against
  4.00 s).

The dispatch guard requires NumPy's floating-point errors ignored and no warning-as-error
filters; NumPy's error state is thread-local. The harness explicitly initializes thread
workers with each requested policy: only its all-ignore row reaches Rust. Process rows
set the caller's policy only and report native reachability as `n/a`, because the recorder
cannot observe another process. A separate default spawned-worker probe reports a mixed
policy, consistent with Python dispatch in that configuration, not every process pool.
A caller-provided pool initializer can change worker policies. The process-versus-thread
ignore timing gap therefore combines backend and scheduler/serialization effects, not
scheduler overhead alone; the probe does not verify routing in the timed process workers.
The guidance above and in [Operational Guidance](docs/xarray_compatibility.md#operational-guidance)
stands.

The full tables -- every entry, the fixed and per-cell fit, thread scaling, and the
Dask policy rows -- are in
[`benchmarks/README.md`](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md#rust-kernels-vs-the-python-reference-across-the-parity-registry-rust-011),
and the harness that produced them is `benchmarks/rust_vs_python.py` (the full-grid
run is its `--netcdf` mode):

```bash
uv run benchmarks/rust_vs_python.py --repeat 5 --write
```

### Dispatch optimization measurements (#1281)

2026-10-10, Apple M5/macOS arm64, Python 3.14.7, NumPy 2.4.2, SciPy 1.17.0,
Rust 1.99.0. Baseline `origin/main` revision `74fef017`; release builds on both
sides, 30 alternating timed calls after warm-up, default console INFO logging.
These are station/small-block measurements, not grid or Xeon results. IQR is
Q3−Q1; it is not a confidence interval. Untouched rows vary with host/run noise.

```bash
# In a baseline worktree at 74fef017, copy the harness and select its original registry:
git show 846fa230:benchmarks/profile_dispatch_overhead.py > benchmarks/profile_dispatch_overhead.py
uv sync && uv run maturin develop --release
uv run python benchmarks/profile_dispatch_overhead.py --label before-full \
  --entries $(uv run python -c 'from tests.parity_registry import ENTRIES_BY_NAME; print(*ENTRIES_BY_NAME)')
# In the optimized worktree, use its shipped harness and rebuild:
uv sync && uv run maturin develop --release
uv run python benchmarks/profile_dispatch_overhead.py --label after-final --profile
```

Native public medians in microseconds (IQR); ratios are Python/native. The
[full table](https://github.com/monocongo/climate_indices/blob/main/benchmarks/README.md#dispatch-and-pm-optimization-1281--1324)
includes both Python timing distributions, raw/prep accounting, realistic and
monotonic day sweeps, profiling ceilings, and source/extension hashes.
Raw samples: `benchmarks/results/before-full-20261010T004923Z.json` and
`benchmarks/results/after-final-20261010T011635Z.json`.
The baseline command excludes supplemental SPEI Pearson, which requires the new KS
seam. These retained measurements predate review fixes to PM cache allocation,
operand validation, alignment, and SciPy warning routing; they have not been rerun
against the final reviewed tree.

| Entry | Native before µs (IQR) | Native after µs (IQR) | Python/native before→after |
| --- | ---: | ---: | ---: |
| `spi_gamma` | 303.21 (11.81) | 318.77 (15.70) | 1.14 → 1.16 |
| `spi_gamma_mean_zero` | 349.19 (16.51) | 367.98 (17.46) | 1.11 → 1.12 |
| `spi_pearson` | 848.31 (24.55) | 329.50 (11.39) | 1.23 → 3.34 |
| `spei_loglogistic` | 216.75 (6.68) | 227.71 (6.43) | 1.57 → 1.56 |
| `spi_gamma_spatial_block` | 377.31 (13.90) | 394.42 (20.36) | 1.16 → 1.18 |
| `spei_gamma` | 306.23 (10.05) | 321.58 (7.51) | 1.17 → 1.17 |
| `percentage_of_normal` | 51.31 (2.89) | 53.54 (3.39) | 1.08 → 1.12 |
| `eddi` | 168.27 (4.80) | 153.31 (6.86) | 1.56 → 1.81 |
| `eddi_spatial_block` | 211.52 (5.56) | 197.29 (4.73) | 1.83 → 2.03 |
| `pci` | 51.48 (2.73) | 51.96 (2.41) | 1.01 → 1.04 |
| `fit_diagnostics` | 1893.42 (73.11) | 1976.29 (27.54) | 1.00 → 1.01 |
| `thornthwaite` | 29.71 (4.69) | 30.33 (5.53) | 220.02 → 223.27 |
| `hargreaves` | 105.06 (7.03) | 107.19 (9.15) | 34.58 → 35.38 |
| `penman_monteith` | 133.12 (4.43) | 100.65 (3.47) | 0.86 → 1.19 |
| `pm_eto_intermediates` | 17.60 (0.55) | 10.08 (0.38) | 0.45 → 0.84 |
| `fire_ffmc` | 217.02 (15.73) | 225.12 (11.57) | 45.02 → 45.30 |
| `fire_duff_moisture_code` | 234.90 (9.81) | 243.75 (11.79) | 36.98 → 37.34 |
| `fire_drought_code` | 221.29 (12.71) | 229.62 (8.28) | 37.51 → 37.78 |
| `fire_kbdi` | 200.00 (13.94) | 204.25 (11.85) | 43.10 → 44.02 |
| `flood_pe` | 21.00 (1.12) | 20.85 (0.61) | 2.33 → 2.48 |
| `flood_edi` | 33.90 (3.75) | 35.44 (1.98) | 2.10 → 2.14 |
| `flood_flood_index` | 34.10 (2.15) | 34.69 (1.83) | 2.21 → 2.21 |
| `flood_api` | 167.08 (8.99) | 177.25 (21.19) | 44.95 → 45.05 |
| `palmer_pdsi` | 170.67 (23.08) | 190.19 (72.61) | 202.50 → 195.18 |
| `palmer_scpdsi` | 244.44 (13.15) | 250.75 (38.44) | 64.66 → 65.70 |

PM prep drops from 39.29 (3.47) to 12.71 (2.57) µs. Scalar operands allocate
no element buffer: at 10⁷ elements the fresh-process prep RSS delta is zero;
whole-call peak is 827.5 MB, including real-array copies and output. The separate
cycling-day PM case reaches 1.84× at 146,000 and 1.88× at 10⁶ elements; the
unchanged monotonic-day case remains ≈1.0×. Eq 6 still loses to NumPy on M5.
No performance threshold or backend exception was inferred from one host.

Before the Pearson port, even an optimistic zero-time entire-GoF ceiling was
1.28× for SPI/SPEI Gamma, 1.01× for log-logistic/PNP, but 2.76–2.78× for Pearson.
Only Pearson calibration CDF→KS D was batched. Python still computes candidate
exact p-values and emits unchanged warnings; near-threshold D values keep the
SciPy oracle. Supplemental SPEI Pearson improves from 847.31 (38.08) to
333.31 (9.61) µs native, with Python/native 1.22→3.28. Spatial GoF and further
fit/transform fusion remain unported. Parity remains `rtol = atol = 1e-10`;
new KS maxima are 3.331e-16 absolute and 1.830e-15 relative.

PCI and `fit_diagnostics` are expected ≈1.0×: validation-bound and SciPy-KS-bound,
respectively. Keep their existing dispatch for consistency; do not promise a
standalone gain or port more trivial work. SIMD was deferred because realistic
cycling PM already exceeds 1.5× at n≥10⁵ here; monotonic PM remains an unfavorable
case. #1281 needs a Xeon repeat and an Eq 6 routing/crossover investigation;
#1324 needs isolated-host repeats of the now-persisted raw-sample harness.

## docs/algorithms.md

# Climate Index Background and Selection Guidance

What each index measures, when to use it, and how to choose parameters. For the
computation steps, formulas, and implementation details, see
{doc}`algorithm-reference`; the scientific literature is listed in that page's
bibliography.

```{contents} Table of Contents
:backlinks: none
:local: true
```

## Overview

The `climate_indices` library implements algorithms for computing drought and climate variability indices. The implementations follow peer-reviewed methodologies, but validation evidence varies by index: some have been validated against reference datasets, while the flood indices have regression and specification checks only, with no adopted external numeric oracle. The algorithms are designed for operational climate monitoring, research applications, and climate impact assessments.

## Standardized Precipitation Index (SPI)

### Overview

The Standardized Precipitation Index (SPI) is a widely used indicator for characterizing meteorological drought on multiple timescales. Developed by McKee, Doesken, and Kleist (1993), SPI transforms precipitation data into standardized units that represent the probability of occurrence relative to the long-term climatological distribution.

SPI values are dimensionless and follow a normal distribution with mean 0 and standard deviation 1, making them comparable across different locations and climate regimes.

### When to Use SPI

**Ideal applications:**

- Precipitation-only drought monitoring
- Locations where temperature data is unavailable or unreliable
- Comparing drought conditions across different climate zones
- Operational drought monitoring systems
- Agricultural drought assessment (short timescales: 1-3 months)
- Hydrological drought assessment (long timescales: 12-24 months)

**Limitations:**

- Does not account for evapotranspiration losses
- May underestimate drought severity in warming climates
- Not suitable for studying temperature-driven drought
- Limited applicability in arid regions with many zero-precipitation months

See {doc}`algorithm-reference` for the SPI computation steps and distribution
parameters.

### Parameter Selection

#### Scale (Timescale)

The scale parameter determines the temporal aggregation period: 1-72 months for
monthly data, 1-2196 days for daily data.

**Common scales and applications:**

- **1-month**: Short-term precipitation anomalies, agricultural impacts
- **3-month**: Seasonal drought, soil moisture deficits
- **6-month**: Medium-term drought, streamflow impacts
- **12-month**: Annual drought, reservoir management
- **24-month**: Long-term drought, groundwater depletion

**Selection guidance:**

- Match scale to the impact domain of interest
- Agricultural applications: 1-6 months
- Hydrological applications: 6-24 months
- Water resource management: 12-48 months

#### Distribution Choice

**Gamma distribution** (recommended default):

- Computationally efficient
- Stable parameter estimation
- Adequate for most precipitation distributions
- Recommended by McKee et al. (1993)

**Pearson Type III distribution** (alternative):

- Better fit for highly skewed data
- More parameters = more flexible
- Requires more data for stable estimation
- Recommended when gamma fit is poor (Kolmogorov-Smirnov test p < 0.05)

#### Calibration Period

The calibration period defines the climatological baseline for "normal" precipitation.

**Requirements:**

- **Minimum length**: 30 years (≥30 years strongly recommended)
- **Maximum missing data**: \<20% of calibration period values
- **Goodness-of-fit validation**: Kolmogorov-Smirnov test at p = 0.05 significance level

**Selection guidance:**

- Use WMO-recommended 30-year climate normals (e.g., 1991-2020)
- Ensure calibration period is representative of current climate
- Avoid periods with major data gaps or quality issues
- Consider updating calibration period periodically for climate change studies

## Standardized Precipitation Evapotranspiration Index (SPEI)

### Overview

The Standardized Precipitation Evapotranspiration Index (SPEI) is an extension of SPI that incorporates the effect of temperature on drought through potential evapotranspiration (PET). Developed by Vicente-Serrano, Begueria, and Lopez-Moreno (2010), SPEI provides a more complete picture of water balance by accounting for both water supply (precipitation) and atmospheric water demand (PET).

SPEI is particularly valuable in warming climates where increased temperatures amplify drought severity through enhanced evaporative demand.

### When to Use SPEI

**Ideal applications:**

- Drought monitoring in warming climates
- Agricultural drought assessment (crop water stress)
- Climate change impact studies
- Regions with significant temperature variability
- Studies requiring both water supply and demand

**Advantages over SPI:**

- Captures temperature-driven drought intensification
- More realistic representation of agricultural water stress
- Better correlation with soil moisture and crop yields
- Suitable for climate change studies

**Limitations:**

- Requires quality temperature data in addition to precipitation
- PET estimation adds uncertainty (method-dependent)
- More complex computation than SPI
- PET methods may not be equally valid in all climates

### SPI vs SPEI Decision Guide

**Choose SPI when:**

- Only precipitation data is available
- Focus is on meteorological drought
- Comparing historical periods with consistent climate
- Operational monitoring with limited data sources
- Temperature data quality is questionable

**Choose SPEI when:**

- Both precipitation and temperature data are available
- Focus is on agricultural or ecological drought
- Studying climate change impacts
- Temperature trends are significant
- Crop water stress is the primary concern

**Performance comparison:**

| Characteristic     | SPI                | SPEI                      |
| ------------------ | ------------------ | ------------------------- |
| Data requirements  | Precipitation only | Precipitation + temp      |
| Drought type       | Meteorological     | Agricultural/ecological   |
| Climate change     | Less sensitive     | More sensitive            |
| Computation        | Simpler            | More complex              |
| Agricultural use   | Good               | Better                    |
| Historical studies | Good               | Better for recent decades |

### Parameter Selection

SPEI uses the same parameter selection guidelines as SPI:

- **Scale**: Same timescale considerations (1-72 months for monthly data, 1-2196 days for daily data)
- **Distribution**: Gamma (default) or Pearson Type III
- **Calibration period**: ≥30 years, \<20% missing data, K-S test p=0.05

**Additional consideration - PET method:**

See the Potential Evapotranspiration section below for guidance on choosing between Thornthwaite and Hargreaves methods.

See {doc}`algorithm-reference` for the SPEI computation steps.

## Potential Evapotranspiration (PET)

### Overview

Potential evapotranspiration (PET) represents the atmospheric water demand—the amount of water that would evaporate and transpire from a reference surface if sufficient water were available. PET is a key input for SPEI and water balance calculations.

This library implements two temperature-based PET methods:

1. **Thornthwaite (1948)**: Monthly timestep, temperature-only
2. **Hargreaves (1985)**: Daily timestep, temperature and radiation

The equations and implementation details are in {doc}`algorithm-reference`.

### Thornthwaite strengths and limitations

**Strengths:**

- Simple, requires only temperature and latitude
- Well-tested for temperate climates
- Computationally efficient
- Suitable when only temperature data is available

**Limitations:**

- Overestimates PET in arid/windy regions
- Underestimates PET in humid regions
- Does not account for wind speed, humidity, or solar radiation
- Monthly timestep only (not suitable for daily analysis)

### Hargreaves strengths and limitations

**Strengths:**

- More accurate than Thornthwaite in many climates
- Daily timestep suitable for high-resolution analysis
- Temperature range is a reasonable proxy for radiation
- Recommended by FAO for data-limited situations

**Limitations:**

- Requires daily min/max temperature (more data than Thornthwaite)
- Assumes temperature range correlates with radiation (not always true)
- May be less accurate in cloudy/humid climates
- Not suitable when only mean temperature is available

### Thornthwaite vs Hargreaves Decision Guide

**Choose Thornthwaite when:**

- Only monthly mean temperature is available
- Working with historical datasets (common format)
- Computational efficiency is critical
- Consistency with legacy analyses is required
- Operating in temperate climates

**Choose Hargreaves when:**

- Daily temperature data (min/max) is available
- Higher temporal resolution is needed
- Operating in semi-arid to arid climates
- Following FAO guidelines for crop water requirements
- More accurate PET estimates are required

**Performance comparison:**

| Characteristic      | Thornthwaite          | Hargreaves        |
| ------------------- | --------------------- | ----------------- |
| Temporal resolution | Monthly               | Daily             |
| Temperature inputs  | Mean only             | Min, max, mean    |
| Additional inputs   | Latitude              | Latitude          |
| Accuracy (general)  | Moderate              | Good              |
| Accuracy (arid)     | Poor (overestimates)  | Good              |
| Accuracy (humid)    | Poor (underestimates) | Moderate          |
| Data requirements   | Minimal               | Moderate          |
| Computational cost  | Low                   | Low               |
| Climate suitability | Temperate             | Semi-arid to arid |

:::{note}
For the most accurate PET estimates, consider using more sophisticated methods
like Penman-Monteith (FAO-56) if wind speed, humidity, and solar radiation data
are available. However, temperature-based methods remain the most practical choice
for large-scale or historical analyses.
:::

## Palmer Drought Indices

### Overview

The Palmer Drought Severity Index (PDSI) and related indices were developed by Wayne Palmer (1965) for the U.S. National Weather Service. The Palmer system includes five related indices:

- **PDSI**: Palmer Drought Severity Index (long-term moisture conditions)
- **PHDI**: Palmer Hydrological Drought Index (groundwater, reservoir levels)
- **PMDI**: Palmer Modified Drought Index (short-term agricultural drought)
- **Z-Index**: Palmer Z-Index (monthly moisture anomaly)
- **scPDSI**: Self-Calibrated PDSI (Wells et al., 2004), available through
  {func}`climate_indices.palmer.scpdsi`, improves spatial comparability by
  calibrating to each location's climate.

These indices are based on a water balance model that accounts for precipitation, evapotranspiration, soil moisture, and runoff.

### When to Use Palmer Indices

**Ideal applications:**

- Operational drought monitoring (NOAA, USDA)
- Agricultural drought assessment
- Long-term water resource planning (PHDI)
- Comparing drought across different U.S. regions (scPDSI)

**Strengths:**

- Comprehensive water balance approach
- Well-established in operational use
- Multiple indices for different applications

**Limitations:**

- Complex calibration requiring soil data
- Fixed temporal scale (not multi-scalar like SPI/SPEI)
- Slow response to emerging drought
- Original PDSI not directly comparable across climates
- Designed for U.S. climates (may need adjustment elsewhere)

:::{note}
The Palmer indices require significant calibration data including soil water holding
capacity, which may not be available for all locations. SPI and SPEI are often
preferred for global applications due to their simpler data requirements and
multi-scalar nature.
:::

The soil-moisture accounting model those indices build on is described in
{doc}`algorithm-reference`.

## Additional Indices

### Precipitation Concentration Index (PCI)

The Precipitation Concentration Index (PCI) quantifies the temporal distribution of precipitation throughout the year (Oliver, 1980).

**Interpretation:**

- PCI < 10: Uniform precipitation distribution
- 10 ≤ PCI < 15: Moderate precipitation concentration
- 15 ≤ PCI < 20: Irregular precipitation distribution
- PCI ≥ 20: Strong precipitation concentration

**Use cases:**

- Characterizing precipitation seasonality
- Climate classification
- Assessing temporal variability of water supply
- Agricultural planning (growing season water availability)

The PCI equation and implementation details are in {doc}`algorithm-reference`.

### Percentage of Normal Precipitation (PNP)

The Percentage of Normal Precipitation (PNP) expresses precipitation as a ratio of the long-term average for a given location and time period.
The output has units `"1"`: `1.0` means 100% of normal, and `0.5` means 50%.
Multiply by 100 to display a percentage; the library does not apply this conversion.

**Interpretation:**

- PNP > 1.0: Above-normal precipitation
- PNP = 1.0: Normal precipitation
- PNP < 1.0: Below-normal precipitation
- PNP < 0.7: Drought conditions (rule of thumb)

**Use cases:**

- Simple, intuitive drought indicator
- Public communication (easily understood)
- Operational monitoring
- Agricultural extension services

**Limitations:**

- Not standardized across different climates
- Skewed distribution (not suitable for statistical analysis)
- Less sophisticated than SPI/SPEI
- Sensitive to calibration period selection

The PNP equation and implementation details are in {doc}`algorithm-reference`.

## Validation Datasets and Methods

### Reference Datasets

The climate_indices library has been validated against several authoritative sources:

**NCAR Climate Data Guide**:

- SPI reference implementations
- Multi-scale validation datasets
- URL: <https://climatedataguide.ucar.edu/climate-data/standardized-precipitation-index-spi>

**NOAA/NCEI Climate Monitoring**:

- Operational SPI, SPEI, PDSI values for U.S. climate divisions
- Monthly updates for verification
- URL: <https://www.ncdc.noaa.gov/temp-and-precip/drought/>

**Global SPEI Database**:

- Vicente-Serrano et al. global SPEI dataset
- 0.5° resolution, 1901-present
- URL: <https://spei.csic.es/database.html>

### Validation Approach

**Test methodology**:

1. **Unit tests**: Test individual functions with known inputs/outputs
2. **Integration tests**: Test complete index calculations
3. **Regression tests**: Compare against reference implementations
4. **Property-based tests**: Validate mathematical properties using Hypothesis
5. **Benchmark tests**: Verify computational efficiency and memory usage

**Validation metrics**:

- Pearson correlation coefficient (r > 0.99 for reference datasets)
- Root mean square error (RMSE)
- Maximum absolute error
- Computation time and memory usage

### External Validation Resources

**NCAR Command Language (NCL)**:

- Mature SPI/SPEI implementations
- URL: <https://www.ncl.ucar.edu/>

**Climate Indices in Python (other libraries)**:

- `climate-indices` (NCAR/UCAR)
- `spei` (CRAN R package, Python wrapper)

:::{tip}
When validating this library's output against other implementations, ensure
consistent parameter choices (distribution, calibration period, scale) and
input data preprocessing (missing value handling, time alignment).
:::

## docs/algorithm-reference.md

# Algorithm Reference

Formulas, computation steps, and implementation details for the indices this
library implements. For what each index is, when to use it, and how to choose
parameters, see {doc}`algorithms`.

```{contents} On this page
:backlinks: none
:local: true
```

## SPI computation

The SPI computation follows these steps:

1. **Input validation and preprocessing**

   - Accept 1-D or 2-D arrays of precipitation values (any units)
   - Clip negative precipitation values to zero
   - Handle missing data appropriately

2. **Temporal aggregation**

   - Compute sliding sums over the specified timescale (e.g., 1, 3, 6, 12, 24 months)
   - For monthly data: reshape to (years, 12)
   - For daily data: reshape to (years, 366) assuming leap year format

3. **Distribution fitting**

   - Fit scaled precipitation to a probability distribution (gamma or Pearson Type III)
   - Compute distribution parameters separately for each calendar month/day
   - Use only calibration period data for parameter estimation

   **Gamma distribution** (default):

   - Two-parameter gamma distribution with shape (α) and scale (β) parameters
   - Parameters estimated using maximum likelihood estimation
   - Suitable for most precipitation distributions

   **Pearson Type III distribution** (alternative):

   - Three-parameter distribution with location, scale, and skewness
   - Parameters estimated using L-moments method (Hosking, 1997)
   - Better for highly skewed precipitation distributions
   - Automatic fallback to gamma if fitting fails

4. **Cumulative probability transformation**

   - Transform fitted values to cumulative probabilities using the distribution's CDF
   - Account for zero precipitation values in probability calculations. A zero
     accumulation cannot be transformed by the fitted continuous distribution, so it
     carries the probability mass `p0` and the fitted value is `p = p0 + (1 − p0)·F(x)`;
     where within that mass a zero is placed is set by `zero_handling` (see below).

5. **Inverse normal transformation**

   - Apply inverse normal (Gaussian) CDF to obtain standardized values
   - This transforms probabilities to standard normal deviates
   - `output_scale="probability"` stops after step 4 and returns the fitted
     cumulative probability in [0, 1] (the probability integral transform, or
     PIT, value); `output_scale="bounded"` returns its `2p - 1` mapping in
     [-1, 1]. Neither is passed through the inverse normal transform

6. **Clipping to valid range**

   - Constrain final SPI values to [-3.09, 3.09] on the default
     `output_scale="normal"` only
   - This range represents probabilities from 0.001 to 0.999

### Zero handling

`spi()` and `standardized_index()` take a `zero_handling` argument that decides where
within the zero mass a zero accumulation scores (ADR-0015). Normal scale:

| Value | A zero scores | Property |
| --- | --- | --- |
| `"classic"` (default) | `Φ⁻¹(p0)` | The top of the zero mass; NOAA/NCEI and SPEIbase convention |
| `"center_of_mass"` | `Φ⁻¹(p0 / 2)` | Ideal probability-scale mean of 1/2 (Stagge et al., 2015) |
| `"mean_zero"` | `−φ(Φ⁻¹(p0)) / p0` | Ideal normal-scale mean of 0; the mean of `Z` conditional on `Z < Φ⁻¹(p0)` (Allen & Otero, 2024) |

`p0` is the zero fraction of the calibration period's non-missing values. Only the
zero positions move: a value with a positive accumulation keeps its classic
transform, and for Pearson Type III a trace value below 0.0005 shares a zero's
position where `0 < p0 < 1`. On `output_scale="probability"` and `"bounded"` both
non-classic modes place a zero at the centre of the mass, `p0 / 2`, mapped to the
requested scale. A mode applies only where `0 < p0 < 1`: where the calibration
period holds no zero or is entirely zero, every mode keeps the classic result.
SPEI is excluded, because its `P − PET` series is offset before fitting and has no
physical zero mass. `standardized_index()` with `Distribution.beta` applies the mode
to both of its boundary masses, placing the mass at 1 at the mirror of the mass at
0 ([ADR-0024](docs/adr/0024-beta-distribution-for-bounded-series.md)).

## SPEI computation

The SPEI computation follows these steps:

1. **Water balance computation**

   - Compute climatic water balance: D = P - PET
   - Where P is precipitation and PET is potential evapotranspiration
   - Add constant offset (+1000 mm) to ensure all values are positive

2. **Apply SPI methodology**

   - Apply the same steps as SPI (aggregation, fitting, transformation)
   - Use water balance (D) instead of precipitation (P) as input
   - Support both gamma and Pearson Type III distributions

3. **Same constraints as SPI**

   - Clip to valid range [-3.09, 3.09] on the default `output_scale="normal"`
     only; `"probability"` and `"bounded"` are not z-clipped
   - Return dimensionless standardized values

**Key differences from SPI:**

- Incorporates temperature effects through PET
- Requires both precipitation and temperature data
- More sensitive to warming trends
- Better represents agricultural drought (crop water stress)

## PET computation

The two temperature-based methods this library implements, and their equations;
for when to choose each, see {doc}`algorithms`.

### Thornthwaite Method

The Thornthwaite method estimates monthly PET using mean temperature and latitude.

**Equation:**

$$
PET = 1.6 \left( \frac{L}{12} \right) \left( \frac{N}{30} \right) \left( \frac{10 T_a}{I} \right)^a
$$

Where:

- *PET* = potential evapotranspiration (mm/month)
- *L* = mean day length for the month (hours)
- *N* = number of days in the month
- *T*{sub}`a` = mean daily air temperature (°C, clipped to ≥0)
- *I* = annual heat index = Σ(*T*{sub}`ai`/5){sup}`1.514` for all 12 months
- *a* = (6.75×10{sup}`-7`)\*I\*{sup}`3` - (7.71×10{sup}`-5`)\*I\*{sup}`2` + (1.792×10{sup}`-2`)\*I\* + 0.49239

**Implementation details:**

- Automatically adjusts for leap years
- Negative temperatures are clipped to 0°C (no evaporation below freezing)
- Day length computed from latitude using solar geometry
- Returns monthly PET values in mm/month

### Hargreaves Method

The Hargreaves method estimates daily PET using temperature range as a proxy for solar radiation.

**Equation (based on FAO-56 equation 52):**

$$
ET_o = 0.0023 (T_{mean} + 17.8) (T_{max} - T_{min})^{0.5} \times 0.408 \times R_a
$$

Where:

- *ET*{sub}`o` = reference evapotranspiration over grass (mm/day)
- *T*{sub}`mean` = mean daily temperature (°C)
- *T*{sub}`max` = maximum daily temperature (°C)
- *T*{sub}`min` = minimum daily temperature (°C)
- *R*{sub}`a` = extraterrestrial radiation (MJ m{sup}`-2` day{sup}`-1`)

**Implementation details:**

- Computes extraterrestrial radiation from day of year and latitude
- Accounts for Earth-Sun distance variation
- Returns daily PET values in mm/day
- Validates temperature relationships (T{sub}`min` ≤ T{sub}`mean` ≤ T{sub}`max`)

## Palmer index computation

The Palmer drought indices use a two-layer soil moisture accounting model:

1. **Water balance model**

   - Compute potential recharge, runoff, and loss for each time step
   - Track moisture in surface and subsurface soil layers
   - Calculate moisture departure from climatically appropriate conditions

2. **CAFEC procedure**

   - Derive Climatically Appropriate For Existing Conditions (CAFEC) coefficients
   - Calibrate moisture departure to local climate

3. **Index calculation**

   - PDSI: Cumulative moisture anomaly with persistence
   - Z-Index: Monthly moisture departure (no persistence)
   - PHDI: Emphasizes long-term moisture deficits
   - PMDI: Responds quickly to short-term changes

## PCI computation

**Equation** (Oliver, 1980):

$$
PCI = \frac{\sum_{i=1}^{12} P_i^2}{\left(\sum_{i=1}^{12} P_i\right)^2} \times 100
$$

Where *P*{sub}`i` is the precipitation in month *i*.

**Implementation details:**

- Requires complete annual cycle (365 or 366 daily values)
- Rejects incomplete years or years with missing data
- Returns single PCI value per year

## PNP computation

**Equation:**

$$
PNP = \frac{P_{observed}}{P_{normal}}
$$

Where:

- *P*{sub}`observed` = observed precipitation for the period
- *P*{sub}`normal` = long-term average precipitation for the same calendar period

The output is a dimensionless ratio with units `"1"`: `1.0` is 100% of normal,
and `0.5` is 50% of normal. Multiply by 100 to display a percentage.
This convention preserves the values returned by earlier releases.
For nonnegative precipitation, the ratio is nonnegative and has no finite upper
bound. A zero or missing normal produces NaN.

**Implementation details:**

- Supports multi-step scales (e.g., 3-month, 6-month, or daily equivalents such as 90-day)
- Computes normals separately for each calendar time step
- Standard calibration period: 1981-2010 (U.S. climate normals)
- Can use any calibration period ≥30 years
- CLI PNP variables carry `valid_min: 0.0` and no `valid_max` because the ratio is unbounded.
- With `--pack`, the CLI selects storage per variable. PNP ratios above 3.2767
  use float32 instead of packed int16, so different timescales can use different storage types.

## Statistical Methods

### Distribution Fitting

**Gamma distribution**:

The two-parameter gamma distribution is the default for SPI and SPEI:

$$
f(x; \alpha, \beta) = \frac{1}{\beta^\alpha \Gamma(\alpha)} x^{\alpha-1} e^{-x/\beta}
$$

Where:

- α = shape parameter
- β = scale parameter
- Γ(α) = gamma function

Parameters are estimated using maximum likelihood estimation (MLE) on the calibration period data for each calendar month independently.

**Pearson Type III distribution**:

The three-parameter Pearson Type III distribution is an alternative for skewed data:

$$
f(x; \mu, \sigma, \gamma) = \frac{1}{\sigma \Gamma(\alpha) \beta^\alpha} (x - \xi)^{\alpha-1} e^{-(x-\xi)/\beta}
$$

Where:

- μ = location parameter
- σ = scale parameter
- γ = skewness parameter

Parameters are estimated using L-moments method (see below) on the calibration period data.

### L-moments Estimation

L-moments (Linear moments) are used to estimate Pearson Type III distribution parameters. L-moments have advantages over conventional moments:

- More robust to outliers
- Better for small sample sizes
- Unbiased estimators
- Exist even when conventional moments do not

**Implementation**:

This library implements the L-moments estimation procedures from Hosking (1997):

1. Compute sample L-moments (λ₁, λ₂, τ₃) from calibration data
2. Estimate Pearson Type III parameters from L-moments using analytical relationships
3. Validate parameter estimates

**Reference**: Hosking, J. R. M. (1997). *FORTRAN Routines for Use with the Method of L-Moments, Version 3*. IBM Research Report RC20525. IBM Research Division, T. J. Watson Research Center, Yorktown Heights, NY.

**Minimum data requirement**: At least 4 non-zero values per calendar month for stable L-moments estimation.

### Goodness-of-Fit Validation

The library performs goodness-of-fit testing to validate that the fitted distribution adequately represents the observed data.

**Kolmogorov-Smirnov (K-S) test**:

- Tests the null hypothesis that data follows the fitted distribution
- Applied separately for each calendar month
- Significance level: α = 0.05
- Warning issued if p < 0.05 (poor fit)

**Quality checks**:

1. **Calibration period length**: Warning if < 30 years

   - Constant: `MIN_CALIBRATION_YEARS = 30`

2. **Missing data threshold**: Warning if > 20% missing

   - Constant: `MISSING_DATA_THRESHOLD = 0.20`

3. **Goodness-of-fit threshold**: Warning if p-value < 0.05

   - Constant: `GOODNESS_OF_FIT_P_VALUE_THRESHOLD = 0.05`

**Fallback strategy**:

A calendar period with fewer than four non-zero values, or whose L-moment fit fails, gets zeroed Pearson Type III parameters and a logged warning; that alone does not trigger a fallback. SPI and `standardized_index()` refit the whole input with gamma, logging `distribution_fallback`, only when the Pearson fit raises a fitting error or its transform loses more than half of the input's valid values. SPEI does not fall back. See {doc}`data_requirements` for what a failed period outputs.

### Run Theory (Event Identification)

Run theory characterizes drought and wet events as *runs*: maximal contiguous
periods on one side of a threshold (Yevjevich, 1967). `climate_indices.runs`
reports detection, event metrics, and period summaries. Severity classes and
spatial aggregation stay downstream.

For a run of elapsed duration $d$ with values $x_i$, of which the set $A$ lies on
the chosen side of a threshold $T$ (all $d$ values when `max_gap=0`):

$$
\text{magnitude} = \sum_{i \in A} |x_i - T|
\qquad
\text{intensity} = \frac{\text{magnitude}}{d}
$$

The **peak** is the minimum value in the run when below the threshold, the
maximum when above. **Interarrival** is the number of time steps from one run's
start to the next run's start, so the final run has none (NaN).

Conventions:

- Direction is explicit (`direction="below"` or `"above"`) and never inferred
  from the threshold's sign, so a threshold such as "above -0.5" is expressible
  and wet runs take their peak on the maximum side.
- A NaN value is not in any run and terminates the run it would otherwise
  continue, so gaps split runs.
- A value exactly equal to the threshold is not in a run.
- `max_gap=0` disables pooling. A positive `max_gap` joins successive runs when
  their intervening gap is no longer than that many input steps.
- Pooling never crosses NaN or masked values. A threshold-equal value can form
  part of a pooled gap, but it contributes no magnitude.
- Duration includes pooled gaps. Magnitude sums only the deviations on the
  selected side of the threshold. Intensity divides that magnitude by elapsed
  duration, including the gaps.
- `min_duration` filters runs after pooling. Interarrival is computed between
  surviving runs only.

The primitive is index-agnostic: a dry SPI spell, a wet SPEI spell, and a KBDI
drought spell (``direction="above"``, since KBDI rises in drought) are the same
run concept on a different series.

#### Period summaries

`runs.summarize_runs()` returns a `RunStatistics` containing a `RunSet` of
retained events and summary arrays. Supply one `periods` label per input step,
such as its year or calendar month. Repeated labels form one group, and the
summary arrays follow sorted unique labels. Without labels, one summary covers
the whole input, including an empty input.

Each complete event belongs to its onset period. The summary reports event
count, total elapsed duration, maximum event magnitude, and valid input count.
An event that crosses a period boundary contributes all its duration and
magnitude to its onset period, so total duration can exceed that period's length.
Frequency is event count divided by the period's non-NaN, unmasked input count.
It is events per valid input step, not the fraction of steps in an event, an
annual rate, or a return period. Infinities retain the original run comparisons
and count as observations; callers must mask them if they represent missing data.

Periods with no events have zero count, duration, and maximum magnitude.
Frequency is NaN when a period has no valid observations. Missing, masked, or
infinite period labels are rejected rather than silently dropped. Missing means
NaN or infinity in a float array, NaT in a datetime64 or timedelta64 array, or
a masked element. NumPy converts mixed string and NaN labels to a string array
that holds the text `"nan"`, so such input is an ordinary label, not missing.
These conventions also apply to
{func}`climate_indices.energy.energy_drought_events`.

```{testcode}
import numpy as np
from climate_indices.runs import summarize_runs

result = summarize_runs(
    [-2, 0, -3, np.nan, -4],
    periods=[2020, 2020, 2021, 2021, 2021],
    max_gap=1,
)
assert result.events.start_index.tolist() == [0, 4]
assert result.count.tolist() == [1, 1]
assert result.total_duration.tolist() == [3, 1]
assert result.max_magnitude.tolist() == [3.0, 3.0]
assert result.frequency.tolist() == [0.5, 0.5]
```

A probability-scale series (`output_scale="probability"`) reports a lower-tail
probability, so its drought runs also use `direction="below"`. The 10%, 5%, and
2.5% one-sided levels (0.1, 0.05, 0.025) correspond to z-scores of about -1.28,
-1.64, and -1.96. Note this is a different ladder from the McKee/WMO drought
classes of -1.0, -1.5, and -2.0, whose probability-scale equivalents are about
0.159, 0.067, and 0.023. Probability-scale outputs are not clipped, so a run can
be counted all the way into the tail.

**Reference**: Yevjevich, V. (1967). *An objective approach to definitions and
investigations of continental hydrologic droughts*. Hydrology Papers 23,
Colorado State University. Drought classes: McKee et al. (1993) and WMO (2012),
both in the bibliography below.

### Compound Event Matching (Concurrent and Sequential)

`climate_indices.compound.compound_events()` matches two run sets that were
identified on the same time axis, for example a drought run set and a heat run
set. Both sources are given as [`RunSet`](#run-theory-event-identification)
objects from `identify_runs()`, plus the shared time axis, its declared cadence,
and one Boolean observation-validity mask. A row names one source run from each
set and never merges or scores them: the result holds indices, durations, lags,
and censoring flags only, so no drought magnitude, fire severity, or heat
intensity is combined into an unrequested severity score.

The time axis must be a complete cadence grid. An absent date is an unknown
observation, not a step to close, so a gap is rejected instead of being
bridged. Runs are validated before matching: each run must hold positive
durations equal to `end_index - start_index + 1`, lie inside the record, be
sorted and nonoverlapping, and contain no unknown observation. A run asserts a
contiguous observed spell, so a run that crosses an unknown step is rejected
rather than silently truncated.

**Concurrent matching** keeps every pair whose inclusive intervals intersect
and reports the intersection, with no one-to-one assignment:

$$
\text{start} = \max(a_{\text{start}}, b_{\text{start}}), \qquad
\text{end} = \min(a_{\text{end}}, b_{\text{end}}), \qquad
\text{overlap\_duration} = \text{end} - \text{start} + 1
$$

Touching endpoints count, so a pair sharing one time step has overlap duration
one. Adjacent pair records are never merged, and a pair is never discarded
because another run already matched, so one long run can match several short
ones and the reverse.

**Sequential matching** keeps every pair with $b_{\text{start}} >
a_{\text{end}}$, so an overlapping pair, or one starting exactly on the first
run's end, is not a follow-on. The lag is the separating distance:

$$
\text{lag} = b_{\text{start}} - a_{\text{end}}
$$

Adjacent runs therefore have lag one. A caller-supplied inclusive window
`lag_window=(low, high)` keeps every pair with `low <= lag <= high`; `None`
matches every lag at or above one. Every step from $a_{\text{end}}$ through
$b_{\text{start}}$ must be observed. A pair separated by an unknown step is not
reported, because a cascade is never inferred across a gap, and the rejected
count is recorded as `n_unknown_bridges` in the result metadata. The reported
interval is the first run's onset through the second run's end, and the lag is
reported separately: the enclosing interval is not an invented joint duration,
and `overlap_duration` is zero because sequential runs are disjoint by
definition.

**Censoring.** A run boundary that touches the record boundary or an unknown
neighbouring step is censored: the physical event may extend beyond the
observed record or across the unknown step. Every row records four flags —
`first_left_censored`, `first_right_censored`, `second_left_censored`, and
`second_right_censored` — so an observed duration is never presented as a
complete physical event duration. Censoring does not remove a run.

Rows appear in first source-run order and then in second source-run order, and
`first_run_id` and `second_run_id` are the caller's own positions, so
`first.start_index[first_run_id]` recovers the source run and its metrics. Every
result array is read-only, and an empty catalog keeps its integer, float, and
Boolean dtypes. The manifest records the calendar, the matching, censoring, and
missingness rules, the effective lag window, the observed and unknown step
counts, the run and pair counts, the rejected unknown-bridge count, and each
source's declared threshold, direction, and `min_duration`. Those three declarations are not retained by `RunSet`, so
any that the caller leaves out is recorded as `"unknown"` and is never
reconstructed from run metrics.

**Limits.** Pair counts are not independent observations and must not feed a
binomial significance test. Matching inherits the run definition of the source
indices and adds no onset probability or causal claim: a lagged association is
a descriptive co-occurrence, not a validated cascade. The catalog has no primary
index CF key, because it is an event catalog and not an index field.

**Reference**: the event contract and its censoring, lag, and all-pairs rules
are recorded in
[ADR-0022](https://github.com/monocongo/climate_indices/blob/main/docs/adr/0022-compound-analysis-contracts.md).

## Complete Bibliography

**SPI and Drought Indices:**

- McKee, T. B., Doesken, N. J., & Kleist, J. (1993). The relationship of drought frequency and duration to time scales. *Proceedings of the 8th Conference on Applied Climatology*, 17-22 January, Anaheim, CA. American Meteorological Society, Boston, MA, 179-184.
- Vicente-Serrano, S. M., Begueria, S., & Lopez-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>

**Potential Evapotranspiration:**

- Thornthwaite, C. W. (1948). An approach toward a rational classification of climate. *Geographical Review*, 38(1), 55-94. <https://doi.org/10.2307/210739>
- Hargreaves, G. H., & Samani, Z. A. (1985). Reference crop evapotranspiration from temperature. *Applied Engineering in Agriculture*, 1(2), 96-99. <https://doi.org/10.13031/2013.26773>
- Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). *Crop evapotranspiration: Guidelines for computing crop water requirements*. FAO Irrigation and Drainage Paper 56. Food and Agriculture Organization of the United Nations, Rome. ISBN 92-5-104219-5. Available at: <http://www.fao.org/3/x0490e/x0490e00.htm>

**Palmer Drought Indices:**

- Palmer, W. C. (1965). *Meteorological Drought*. U.S. Weather Bureau Research Paper No. 45. Washington, D.C.
- Wells, N., Goddard, S., & Hayes, M. J. (2004). A Self-Calibrating Palmer Drought Severity Index. *Journal of Climate*, 17(12), 2335-2351. <https://doi.org/10.1175/1520-0442(2004)017%3C2335:ASPDSI%3E2.0.CO;2>

**Precipitation Concentration Index:**

- Oliver, J. E. (1980). Monthly precipitation distribution: A comparative index. *The Professional Geographer*, 32(3), 300-309. <https://doi.org/10.1111/j.0033-0124.1980.00300.x>

**Statistical Methods:**

- Hosking, J. R. M. (1997). *FORTRAN Routines for Use with the Method of L-Moments, Version 3*. IBM Research Report RC20525. IBM Research Division, T. J. Watson Research Center, Yorktown Heights, NY. (Note: Available through various online archives and the R package `lmomco`)

**Additional Resources:**

- World Meteorological Organization (WMO). (2012). *Standardized Precipitation Index User Guide* (WMO-No. 1090). Geneva, Switzerland. Available at: <https://library.wmo.int/doc_num.php?explnum_id=7768>
- American Meteorological Society. (2020). Drought. *Glossary of Meteorology*. <https://glossary.ametsoc.org/wiki/Drought>
