# climate_indices

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

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

## README.md

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

# climate_indices

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

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

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

### Indices

**Drought and moisture balance**

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

**Evapotranspiration**

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

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

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

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

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

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

Planned, none shipped yet in the package:

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

**Supporting tools**

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

### Highlights

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

### Project goals

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

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


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

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

## Installation

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

Install the released 2.x package from PyPI:

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

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

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

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

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

## Developer Workflow

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

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

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

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

## Release Recipe

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

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

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

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

### Maintainer Quick Commands

Read-only preflight:

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

Safe PR branch setup:

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

Approval-required release tag commands:

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

## Supported Python Versions

Supports Python **3.10–3.14**.

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

### Version Support Policy

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

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

### API Stability

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

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

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

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

## Validation

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

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

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

## Upgrading to 3.0.0

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

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

## VALIDATION.md

# Validation Status

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

## Test Commands

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

## Per-Index Evidence

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

## Flood Evidence Classification

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

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

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

## EDDI Fixture Policy

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

Required fixture layout:

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

Each directory must contain:

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

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

## SPEI Plausibility Classification

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

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

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

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

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

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

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

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

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

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

## SPEI Log-Logistic Cross-Implementation

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

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

## Palmer Validation Classification

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

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

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

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

## scPDSI Calibration-Anchor Measurement

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

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

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

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

## References

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

## docs/xarray_compatibility.md

# xarray Compatibility Matrix

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

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

## Stateful fire indices

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

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

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

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

## Operational Guidance

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

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

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

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

## Chunking guidance for gridded indices

Three rules fix the shape of a chunked input:

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

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

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

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

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

Recommendations for monthly grids:

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

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

## docs/algorithm_refs/eddi.md

# Evaporative Demand Drought Index (EDDI)

## Scope

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

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

## Public API

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

## Validation

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

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

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

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

## Tolerance

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

## References

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

## docs/algorithm_refs/pet.md

# Potential Evapotranspiration (PET)

## Scope

`climate_indices` implements three PET methods:

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

## Public API

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

## FAO-56 Penman-Monteith

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

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

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

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

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

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

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

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

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

## Validation

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

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

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

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

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

## Tolerance

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

## References

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

## docs/algorithm_refs/palmer.md

# Palmer Drought Indices

## Scope

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

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

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

## Algorithm Notes

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

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

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

### Exact-zero comparisons

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

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

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

## Validation

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

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

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

## Release Decision

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

## References

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

## 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)
- **Current 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-core/            # pure Rust numerical kernels
crates/climate-py/              # PyO3/NumPy bindings -> climate_indices._native
src/climate_indices/_native.pyi # type stub for the extension
```

- **`climate-core`** is pure numerical Rust. It has no PyO3, NumPy bindings,
  Python exceptions, or CPython assumptions, receives already-validated
  numbers, and returns numbers. It could later serve other Rust consumers, but
  it is not published to crates.io.
- **`climate-py`** is the only crate that knows about Python. It converts
  arrays and errors at the boundary and holds no climate algorithm.
- **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-core`) |
|---|---|
| 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-core`) |
|---|---|
| `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-core`) |
|---|---|
| `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-core`) |
|---|---|
| `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`, 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-core/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 |

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-core/src/`, in the module that owns the index
family, and declare that module in `crates/climate-core/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-core` 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
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

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

:::{warning}
**Beta Feature**

The xarray DataArray API described in this guide is **beta** in 3.0.0 and is
promoted no earlier than 3.1.0
([ADR-0012](docs/adr/0012-xarray-api-stays-beta-through-3.0.0.md)). The interface
(parameter inference, metadata attributes, coordinate handling) may change with a
minor version, never in a patch release. The underlying computation results are
intended to be identical to the stable NumPy API.
:::

:::{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, beta)

**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 xarray API is beta in 3.0.0. The NumPy API remains the stable
integration surface; this page covers the Dask execution path for the xarray
adapters described in the [xarray compatibility matrix](docs/xarray_compatibility.md).
:::

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 percentage of the long-term average for a given location and time period.

**Interpretation:**

- PNP > 100%: Above-normal precipitation
- PNP = 100%: Normal precipitation
- PNP < 100%: Below-normal precipitation
- PNP < 70%: 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.

## 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}} \times 100
$$

Where:

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

**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

## 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**:

If Pearson Type III fitting fails (due to insufficient data or numerical issues), the library automatically falls back to gamma distribution with warning messages in the log.

### 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 and metrics only; severity classes, spatial aggregation, and
period slicing stay downstream.

For a run of $d$ values $x_i$ on the chosen side of a threshold $T$:

$$
\text{magnitude} = \sum_{i=1}^{d} |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.
- `min_duration` filters runs after identification, and 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.

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.

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