Add flexibility to handling probability of zero-precipitation events in Standardized Indices (SPI)

[Hem-W]

Addressing a Problem?

The calculation of the Standardized Precipitation Index (SPI) involves transforming a fitted probability distribution (often Gamma) into a standard normal distribution. A critical step in this process is handling "zero-inflated" data (months with no precipitation).

Different methodologies exist for assigning a probability to these zero events before the normal transformation:

• Traditional/McKee Method ("Upper"): Assigns the full cumulative probability of zero, $P(x=0) = p_0$. This is known to introduce a positive bias in the mean SPI, especially in arid regions where $p_0$ is large, as it effectively treats all zero events as "wetter" than they should be.

• Center of Mass Method ("Center"): Assigns the probability $p_0 / 2$. This is recommended by Stagge et al. (2015) to correct the bias and ensure the SPI distribution remains centered at 0.

Reference:
Stagge et al. (2015) Candidate Distributions for Climatological Drought Indices (SPI and SPEI)

Potential Solution

Adding a prob_zero_method parameter to standardized_index and standardized_precipitation_index, which allow users to specify whether use "upper", "center" or customized callable as the probability corresponding to zero-precipitation.

Additional context

# %%
import xarray as xr
import numpy as np
import scipy.stats
from xclim.indices import standardized_precipitation_index
from xclim.indices.stats import standardized_index_fit_params

# Use cftime to create a NOLEAP calendar (365 days every year)
# This ensures a strict 365-day periodicity.
time = xr.cftime_range("2000-01-01", "2010-12-31", calendar="noleap", freq="D")
pr = xr.DataArray(
    np.random.rand(len(time)) * 10,
    coords={"time": time},
    dims="time",
    name="pr",
    attrs={"units": "mm/day"}
)

# Force DOY 180 to be ZERO for ALL 11 years.
pr[{"time": slice(179, None, 365)}] = 0

# Parameters for fitting
fit_params = dict(
    freq=None,
    window=1,
    dist="gamma",
    method="ML",
    doy_bounds=(180, 180),
)

# 1. Get fitting parameters
params = standardized_index_fit_params(pr, zero_inflated=True, **fit_params)
prob_of_zero = params.prob_of_zero.sel(dayofyear=180).values
print(f"Probability of Zero (p0): {prob_of_zero:.4f}")

# 2. Compute SPI
# Original zero-precipitation probability: "upper" (Biased) -> P = p0
spi_upper = standardized_precipitation_index(
    pr, params=params, **fit_params
)

# Verification: Check the SPI value for the first occurrence (index 179)
t_idx = 179
val_upper = spi_upper.isel(time=t_idx).values

print(f"\nResults for Zero-Precipitation Day (p0=1.0):")
print(f"Upper (p0):      SPI={val_upper:.4f}")

This will leads to Probability of Zero (p0): 1.0000 and SPI=8.2100, which may leads to misleading "extreme wet".

Contribution

• I would be willing/able to open a Pull Request to contribute this feature.

Code of Conduct

• I agree to follow this project's Code of Conduct