Skip to content

Sub-Monthly SPI/SPEI Support #3

Open
Open
Sub-Monthly SPI/SPEI Support#3
Labels
enhancementNew feature or request

Description

@bennyistanto
bennyistanto
opened on Apr 3, 2026

Summary

Add support for computing SPI and SPEI from dekadal, pentadal, weekly, and daily input data, where the accumulation scale remains in months (e.g., SPI-12 = 12-month accumulation) but the output is produced at the input's native temporal resolution.

Primary focus: SPI from CHIRPS data (available as daily, pentad, and dekad).
Secondary: SPEI from daily ERA5-Land, with temporal aggregation helpers to match sub-monthly precipitation inputs.

Example: SPI-12 from CHIRPS dekadal → an SPI value every dekad, each representing a 12-month drought signal. "SPI-12 as of Dekad 2 March 2026" covers 21 Mar 2025 – 20 Mar 2026.

Motivation

Currently, users computing sub-monthly SPI must:

  1. Split dekadal data into 3 separate streams (dekad 1, 2, 3)
  2. Aggregate each to monthly
  3. Compute SPI separately for each stream
  4. Merge results back (cdo mergetime)

This is error-prone, slow, and wastes disk space. Native sub-monthly support would:

  • Accept raw dekadal/pentadal/daily data directly
  • Produce a continuous sub-monthly SPI series in one call
  • Detect drought onset/recovery between monthly boundaries
  • Enable weekly/dekadal monitoring products from daily inputs

Data Source Context

Dataset Available Frequencies Primary Index Notes
CHIRPS Daily, Pentad (5-day), Dekad (10-day), Monthly SPI Most common use case
IMERG Half-hourly, Daily, Monthly SPI Can aggregate to any frequency
ERA5-Land Hourly, Daily, Monthly SPEI (with PET) Provides T, PET for SPEI
CMORPH 30-min, Daily, Monthly SPI High-resolution precip
MSWEP 3-hourly, Daily, Monthly SPI Multi-source ensemble

SPI is the primary target - precipitation-only datasets are widely available at sub-monthly frequency. SPEI requires matching PET/temperature data, which is mostly available as daily (ERA5-Land) or monthly. For dekadal/pentadal SPEI, users would need to aggregate daily ERA5-Land to match the precipitation frequency - this requires temporal aggregation helpers (see Phase 3).

Design Principle

Scale is always in months. The temporal resolution of the input determines:

  • How many time steps the accumulation window spans
  • How many distributions are fitted (one per period-of-year)
  • How frequently an output value is produced
Input Frequency Periods/Year Scale=12 → Steps Distributions Fitted
Monthly 12 12 12
Dekadal 36 36 36
Pentadal 73 73 73
Weekly 52 ~52 52
Daily 366 ~365 366

Why pentadal = 73? Each month has 6 pentads (days 1–5, 6–10, 11–15, 16–20, 21–25, 26–end). 12 months × 6 = 72, but the last pentad of each month varies in length (3–6 days), and some definitions use 73 to cover leap years consistently. CHIRPS uses 72 pentads/year (fixed 6 per month). We adopt 73 to safely accommodate all conventions, with NaN padding for the unused slot if needed.

Update: After further review, CHIRPS pentad uses exactly 72 pentads/year (6 per month, last pentad absorbs remaining days). We should use pentadal = 72.

Architecture (Current State)

The core algorithm is already generic. The data flow is:

spi(scale=12, periodicity=monthly)
  → sum_to_scale(data, scale=12)          # Rolling sum of 12 time steps
  → reshape to (years, periods_per_year)  # (years, 12)
  → fit distribution per period           # 12 distributions
  → transform per period                  # 12 transformations

sum_to_scale(), _rolling_sum_3d(), distribution fitting, and run theory are all parameterized by integers - no month-specific logic. The changes are mostly about:

  1. Converting "scale in months" → "scale in time steps"
  2. Expanding the Periodicity enum
  3. Relaxing a few hardcoded validations
  4. Adding date arithmetic for sub-monthly periods
  5. Temporal aggregation helpers for SPEI data alignment

Implementation Plan

Phase 1: Core Infrastructure

1.1 Expand Periodicity Enum

File: src/config.py

class Periodicity(Enum):
    monthly  = 12
    dekadal  = 36
    pentadal = 72
    weekly   = 52
    daily    = 366

    def unit(self) -> str:
        return {
            Periodicity.monthly:  "month",
            Periodicity.dekadal:  "dekad",
            Periodicity.pentadal: "pentad",
            Periodicity.weekly:   "week",
            Periodicity.daily:    "day",
        }[self]

    def steps_per_month(self) -> float:
        """Number of time steps per month for this periodicity."""
        return {
            Periodicity.monthly:  1.0,
            Periodicity.dekadal:  3.0,
            Periodicity.pentadal: 6.0,
            Periodicity.weekly:   52.0 / 12.0,  # ~4.333
            Periodicity.daily:    366.0 / 12.0,  # ~30.5
        }[self]

    def scale_to_steps(self, scale_months: int) -> int:
        """Convert scale in months to number of time steps."""
        return round(scale_months * self.steps_per_month())

    @staticmethod
    def from_string(s: str) -> 'Periodicity':
        try:
            return Periodicity[s.lower()]
        except KeyError:
            valid = ', '.join(f"'{p.name}'" for p in Periodicity)
            raise ValueError(
                f"Invalid periodicity: '{s}'. Must be one of: {valid}."
            )

Scale conversion table (scale_to_steps):

Periodicity scale=1 scale=3 scale=6 scale=12 scale=24
monthly 1 3 6 12 24
dekadal 3 9 18 36 72
pentadal 6 18 36 72 144
weekly 4 13 26 52 104
daily 31 91 183 366 731

Note on rounding: Dekadal and pentadal conversions are exact integers. Weekly and daily use round(). The rounding error is ≤3.5 days, acceptable for drought monitoring. A calendar-aware approach is discussed in Risks & Edge Cases.

1.2 Add Scale Conversion Layer

File: src/compute.py

Add conversion at the entry point of compute_index_parallel() and compute_spi_1d():

def compute_index_parallel(
    data, scale, ..., periodicity, ...
):
    # Convert scale from months to time steps for sub-monthly data
    if periodicity != Periodicity.monthly:
        actual_scale = periodicity.scale_to_steps(scale)
    else:
        actual_scale = scale

    # Rest of function uses actual_scale instead of scale
    scaled_data = _rolling_sum_3d(data, actual_scale, dtype)
    ...

Decision: The scale parameter in spi() always means months. The conversion happens internally. Users think "SPI-12" or "SPI-3", never "sum 36 dekads".

Files affected: compute.py (3 functions: compute_index_parallel, compute_spi_1d, compute_spei_1d), chunked.py (2 functions: compute_spi_chunked, compute_spei_chunked).

1.3 Fix Hardcoded Validation

File: src/utils.py, validate_array() line 146

# Before
if values.shape[1] not in (12, 366):

# After
valid_periods = tuple(p.value for p in Periodicity)
if values.shape[1] not in valid_periods:

1.4 Variable Naming

The existing VAR_NAME_PATTERN = "{index}_{distribution}_{scale}_{periodicity}" auto-generates:

spi_gamma_12_dekad     # SPI-12 from dekadal data
spi_gamma_3_pentad     # SPI-3 from pentadal data
spi_gamma_6_week       # SPI-6 from weekly data
spi_gamma_1_day        # SPI-1 from daily data

However, long_name currently generates "..., 12-dekad" which could confuse users. Update get_long_name() to clarify:

# Before
f"{index_full} ({dist_name}), {scale}-{periodicity.unit()}"
# "Standardized Precipitation Index (Gamma), 12-dekad"

# After - for sub-monthly periodicities, show both
if periodicity == Periodicity.monthly:
    return f"{index_full} ({dist_name}), {scale}-month"
else:
    return f"{index_full} ({dist_name}), {scale}-month scale, {periodicity.unit()}al resolution"
# "Standardized Precipitation Index (Gamma), 12-month scale, dekadal resolution"

Phase 2: Date Handling

2.1 Dekadal Date Arithmetic

File: src/utils.py, compute_time_values()

Dekad definitions (WMO standard):

Dekad Day range
1 1–10
2 11–20
3 21–end of month 
elif periodicity == Periodicity.dekadal:
    dekad_starts = [1, 11, 21]
    current_year, current_month = initial_year, initial_month
    for i in range(total_periods):
        dekad_in_month = i % 3  # 0, 1, 2
        day = dekad_starts[dekad_in_month]
        current_date = datetime(current_year, current_month, day)
        days[i] = (current_date - start_date).days
        if dekad_in_month == 2:  # advance to next month after dekad 3
            current_month += 1
            if current_month > 12:
                current_month = 1
                current_year += 1

2.2 Pentadal Date Arithmetic

Pentad definitions (CHIRPS convention):

Pentad Day range
1 1–5
2 6–10
3 11–15
4 16–20
5 21–25
6 26–end of month 
elif periodicity == Periodicity.pentadal:
    pentad_starts = [1, 6, 11, 16, 21, 26]
    current_year, current_month = initial_year, initial_month
    for i in range(total_periods):
        pentad_in_month = i % 6  # 0-5
        day = pentad_starts[pentad_in_month]
        current_date = datetime(current_year, current_month, day)
        days[i] = (current_date - start_date).days
        if pentad_in_month == 5:  # advance to next month after pentad 6
            current_month += 1
            if current_month > 12:
                current_month = 1
                current_year += 1

2.3 Weekly Date Arithmetic

elif periodicity == Periodicity.weekly:
    current_date = datetime(initial_year, initial_month, 1)
    for i in range(total_periods):
        days[i] = (current_date - start_date).days
        current_date += timedelta(days=7)

Edge case: 52 vs 53 weeks/year. Pad 53rd week with NaN like daily pads Feb 29.

2.4 Daily Date Arithmetic

Already implemented. No changes needed.

Phase 3: Temporal Aggregation Helpers

These are essential for two use cases:

  1. SPI from daily data - user may want dekadal or weekly SPI, needs to aggregate daily precip first
  2. SPEI from sub-monthly precip - ERA5-Land provides daily PET/temperature, must aggregate to match dekadal/pentadal precip frequency

File: src/utils.py (new section)

3.1 Daily to Dekadal

def aggregate_to_dekadal(
    da: xr.DataArray,
    method: str = 'sum'
) -> xr.DataArray:
    """
    Aggregate daily data to dekadal (WMO 10-day periods).

    Dekad 1: days 1-10, Dekad 2: days 11-20, Dekad 3: days 21-end.

    :param da: daily DataArray with time coordinate
    :param method: 'sum' for precipitation, 'mean' for temperature/PET
    :return: dekadal DataArray (36 values per year)
    """
    # Assign dekad labels: day 1-10 → dekad 1, 11-20 → 2, 21+ → 3
    days = da.time.dt.day.values
    dekad_label = np.where(days <= 10, 1, np.where(days <= 20, 2, 3))

    # Create grouping key: year-month-dekad
    ym = da.time.dt.strftime('%Y-%m').values
    group_key = [f"{y}-{d}" for y, d in zip(ym, dekad_label)]

    da_grouped = da.assign_coords(dekad_group=('time', group_key))
    if method == 'sum':
        return da_grouped.groupby('dekad_group').sum('time')
    else:
        return da_grouped.groupby('dekad_group').mean('time')

3.2 Daily to Pentadal

def aggregate_to_pentadal(
    da: xr.DataArray,
    method: str = 'sum'
) -> xr.DataArray:
    """
    Aggregate daily data to pentadal (5-day periods, 6 per month).

    Pentads: 1-5, 6-10, 11-15, 16-20, 21-25, 26-end.
    Last pentad absorbs remaining days (3-6 days depending on month).

    :param da: daily DataArray with time coordinate
    :param method: 'sum' for precipitation, 'mean' for temperature/PET
    :return: pentadal DataArray (72 values per year)
    """
    days = da.time.dt.day.values
    pentad_label = np.where(days <= 5, 1,
                   np.where(days <= 10, 2,
                   np.where(days <= 15, 3,
                   np.where(days <= 20, 4,
                   np.where(days <= 25, 5, 6)))))

    ym = da.time.dt.strftime('%Y-%m').values
    group_key = [f"{y}-{p}" for y, p in zip(ym, pentad_label)]

    da_grouped = da.assign_coords(pentad_group=('time', group_key))
    if method == 'sum':
        return da_grouped.groupby('pentad_group').sum('time')
    else:
        return da_grouped.groupby('pentad_group').mean('time')

3.3 Daily to Weekly

def aggregate_to_weekly(
    da: xr.DataArray,
    method: str = 'sum'
) -> xr.DataArray:
    """
    Aggregate daily data to weekly (7-day periods).

    :param da: daily DataArray with time coordinate
    :param method: 'sum' for precipitation, 'mean' for temperature/PET
    :return: weekly DataArray (~52 values per year)
    """
    if method == 'sum':
        return da.resample(time='7D').sum()
    else:
        return da.resample(time='7D').mean()

3.4 Use Case: Dekadal SPEI from CHIRPS Dekad + ERA5-Land Daily

import xarray as xr
from utils import aggregate_to_dekadal
from indices import spei

# Load data
precip_dekad = xr.open_dataset('chirps_dekad.nc')['precip']  # Already dekadal
pet_daily = xr.open_dataset('era5land_daily.nc')['pet']       # Daily

# Aggregate PET to match precipitation frequency
pet_dekad = aggregate_to_dekadal(pet_daily, method='sum')  # Sum PET over dekad

# Compute SPEI-12 at dekadal resolution
spei_12 = spei(precip_dekad, pet=pet_dekad, scale=12, periodicity='dekadal')

Phase 4: Auto-Detection

4.1 Auto-Detect Periodicity from xarray Input

File: src/utils.py

def detect_periodicity(time_coord) -> Periodicity:
    """Infer periodicity from time coordinate spacing."""
    if len(time_coord) < 2:
        return Periodicity.monthly  # fallback

    # Compute median time delta in days
    deltas = np.diff(time_coord.values).astype('timedelta64[D]').astype(int)
    median_delta = np.median(deltas)

    if median_delta >= 28:       # 28-31 days → monthly
        return Periodicity.monthly
    elif median_delta >= 9:      # 9-11 days → dekadal
        return Periodicity.dekadal
    elif median_delta >= 6:      # 6-8 days → weekly
        return Periodicity.weekly
    elif 4 <= median_delta <= 5: # 4-5 days → pentadal
        return Periodicity.pentadal
    elif median_delta <= 2:      # 1-2 days → daily
        return Periodicity.daily
    else:
        raise ValueError(
            f"Cannot auto-detect periodicity (median delta: {median_delta} days). "
            f"Please specify periodicity explicitly."
        )

Usage in spi():

def spi(precip, scale=12, periodicity='auto', ...):
    if periodicity == 'auto':
        periodicity = detect_periodicity(precip.time)
    elif isinstance(periodicity, str):
        periodicity = Periodicity.from_string(periodicity)

Note: Default periodicity changes from Periodicity.monthly to 'auto'. This is backward-compatible - monthly data auto-detects as monthly. But must test carefully.

Phase 5: Calibration Period Handling

5.1 Year-Based Calibration (No Changes Needed)

The current calibration logic slices the reshaped (years, periods_per_year) array by year indices:

cal_start_idx = calibration_start_year - data_start_year
cal_end_idx = calibration_end_year - data_start_year + 1
calib_data = scaled_4d[cal_start_idx:cal_end_idx, period_idx, :, :]

This works for any periodicity because:

  • Reshaping to (years, periods_per_year) handles the mapping
  • cal_start_idx and cal_end_idx are year indices
  • All periods within a calibration year are included

5.2 Input Length Validation

For sub-monthly data: n_time should ideally be divisible by periods_per_year. If not, reshape_to_2d() already pads with NaN.

Current reshape_to_2d() handles this generically - no changes needed.

Phase 6: Run Theory & Downstream

6.1 Run Theory

identify_runs(), identify_events(), calculate_events_spatial() are all index-based - duration is "number of time steps." This is already correct.

Addition: Add metadata attributes to run theory output so users know the time unit:

# In calculate_events_spatial():
ds['mean_duration'].attrs['units'] = periodicity.unit() + 's'  # "dekads", "pentads", etc.

Add a conversion helper for users who want months:

def duration_to_months(duration_steps: int, periodicity: Periodicity) -> float:
    """Convert duration in time steps to approximate months."""
    return duration_steps / periodicity.steps_per_month()
Periodicity duration=6 → months
monthly 6.0
dekadal 2.0
pentadal 1.0
weekly 1.38
daily 0.20

6.2 GeoTIFF Export

No changes needed - time-agnostic.

6.3 Visualization

No code changes needed. Sub-monthly data plots correctly with more data points on the time axis. Consider adding optional rolling mean smoothing for daily plots (noisy).

Phase 7: Documentation & Naming

7.1 Variable Naming Convention

# Monthly (current, unchanged)
spi_gamma_12_month        → SPI-12 monthly

# Dekadal
spi_gamma_12_dekad        → SPI-12 at dekadal resolution

# Pentadal
spi_gamma_3_pentad        → SPI-3 at pentadal resolution

# Weekly
spi_gamma_6_week          → SPI-6 at weekly resolution

# Daily
spi_gamma_1_day           → SPI-1 at daily resolution

7.2 Long Name / Metadata

long_name = "Standardized Precipitation Index (Gamma), 12-month scale, dekadal resolution"

# Additional variable attributes
scale = 12             # Always in months
scale_unit = "month"
output_resolution = "dekadal"
accumulation_steps = 36  # Actual number of time steps in accumulation window

7.3 User-Facing API Examples

# ── SPI from CHIRPS products ──────────────────────────────────────

# Dekadal SPI-12 from CHIRPS dekadal
spi_12 = spi(chirps_dekad, scale=12, periodicity='dekadal')

# Pentadal SPI-3 from CHIRPS pentad
spi_3 = spi(chirps_pentad, scale=3, periodicity='pentadal')

# Daily SPI-1 from CHIRPS daily
spi_1 = spi(chirps_daily, scale=1, periodicity='daily')

# Auto-detect periodicity (works for any input)
spi_12 = spi(chirps_dekad, scale=12)  # auto-detects dekadal


# ── Dekadal SPI from CHIRPS daily (aggregate first) ──────────────

from utils import aggregate_to_dekadal

precip_dekad = aggregate_to_dekadal(chirps_daily, method='sum')
spi_12 = spi(precip_dekad, scale=12, periodicity='dekadal')


# ── Dekadal SPEI from CHIRPS dekad + ERA5-Land daily ─────────────

from utils import aggregate_to_dekadal
from indices import spei

precip_dekad = xr.open_dataset('chirps_dekad.nc')['precip']
pet_daily = xr.open_dataset('era5land_pet_daily.nc')['pev']

# Aggregate daily PET to dekadal to match precip
pet_dekad = aggregate_to_dekadal(pet_daily, method='sum')

spei_12 = spei(precip_dekad, pet=pet_dekad, scale=12, periodicity='dekadal')

7.4 Documentation Updates

File What to add
docs/user-guide/spi.qmd Sub-monthly SPI examples with CHIRPS
docs/user-guide/spei.qmd Note on SPEI feasibility, ERA5-Land aggregation
docs/technical/methodology.qmd Scale conversion explanation, sub-monthly fitting approach
docs/get-started/configuration.qmd New periodicity options
docs/technical/cf-compliance.qmd Sub-monthly variable names
docs/tutorials/ (new) Tutorial: dekadal SPI from CHIRPS

Phase 8: Testing

8.1 Unit Tests

Test Description
test_scale_to_steps Verify: scale=12 dekadal → 36, pentadal → 72, weekly → 52
test_dekadal_reshape 36 dekads/year reshape and padding
test_pentadal_reshape 72 pentads/year reshape and padding
test_weekly_reshape 52 weeks/year reshape, 53-week year edge case
test_monthly_unchanged Regression: existing monthly behavior preserved
test_daily_unchanged Regression: existing daily behavior preserved
test_detect_periodicity Auto-detection from time coordinate for all frequencies
test_dekadal_spi_1d SPI-12 from 1D dekadal series
test_dekadal_spi_3d SPI-12 from 3D dekadal gridded data
test_pentadal_spi SPI-3 from pentadal data
test_calibration_dekadal Calibration period works with 36 periods/year
test_params_save_load_dekadal Fitting parameters roundtrip for dekadal
test_aggregate_daily_to_dekadal Daily → dekadal aggregation (sum and mean)
test_aggregate_daily_to_pentadal Daily → pentadal aggregation
test_aggregate_daily_to_weekly Daily → weekly aggregation
test_run_theory_dekadal Event duration in dekads, conversion to months

8.2 Validation Against Monthly

Using CHIRPS monthly as ground truth:

  1. Take CHIRPS dekadal, compute SPI-12 at dekadal resolution
  2. Extract dekad 3 of each month (end-of-month value)
  3. Compare against SPI-12 computed from CHIRPS monthly
  4. Should be close (not identical - slight window alignment differences)

Same validation for pentadal (extract pentad 6 of each month).

8.3 Sample Data

Dataset Source Coverage Purpose
CHIRPS dekadal CHC UCSB Small region, 3+ years Primary SPI testing
CHIRPS pentad CHC UCSB Same region Pentadal SPI testing
CHIRPS daily CHC UCSB Same region Aggregation + daily SPI
CHIRPS monthly CHC UCSB Same region Validation baseline
ERA5-Land daily PET CDS Same region SPEI aggregation testing

Files to Modify

File Changes Phase
src/config.py Expand Periodicity enum: add dekadal, pentadal, weekly. Add steps_per_month(), scale_to_steps() methods 1
src/utils.py Fix validate_array() hardcoded check. Add dekadal/pentadal/weekly to compute_time_values(). Add detect_periodicity(). Add aggregate_to_dekadal(), aggregate_to_pentadal(), aggregate_to_weekly() 1, 2, 3, 4
src/compute.py Add scale conversion in compute_index_parallel(), compute_spi_1d(), compute_spei_1d() 1
src/indices.py Support periodicity='auto' in spi() and spei() 1, 4
src/chunked.py Pass scale conversion through (same pattern as compute.py) 1
src/runtheory.py Add duration unit metadata to output attributes 6
src/__init__.py Export new aggregation functions and updated Periodicity 1
docs/ Update methodology, user guide, examples, CF compliance, new tutorial 7

Implementation Order

Phase 1 (Core)              ← Must do first, enables everything else
  ├── 1.1 Periodicity enum (dekadal, pentadal, weekly)
  ├── 1.2 Scale conversion layer (compute.py, chunked.py)
  ├── 1.3 Fix validate_array() hardcoded check
  └── 1.4 Variable naming / long_name update

Phase 2 (Dates)              ← Needed for CF-compliant output
  ├── 2.1 Dekadal date arithmetic
  ├── 2.2 Pentadal date arithmetic
  └── 2.3 Weekly date arithmetic

Phase 3 (Aggregation)        ← Essential for SPEI + daily→sub-monthly
  ├── 3.1 aggregate_to_dekadal()
  ├── 3.2 aggregate_to_pentadal()
  └── 3.3 aggregate_to_weekly()

Phase 4 (Auto-detection)     ← Nice-to-have, improves usability
  └── 4.1 detect_periodicity()

Phase 5 (Calibration)        ← Verify only, likely no code changes

Phase 6 (Downstream)         ← Run theory labeling
  └── 6.1 Duration unit metadata

Phase 7 (Docs)               ← After implementation is stable

Phase 8 (Tests)              ← Throughout, formal suite at end

Risks & Edge Cases

Risk Impact Mitigation
Pentadal: 72 vs 73 per year Reshape padding CHIRPS uses 72 (6 per month). Use 72 as enum value. If other datasets use 73, handle via padding.
Weekly: 52 vs 53 weeks/year ~Every 5-6 years has 53 weeks Pad 53rd week with NaN, same approach as daily Feb 29
Daily: scale=1 → 31 steps round(1 × 30.5) = 31 Acceptable for SPI-1. Document that scale=1 daily ≈ 1 month
Rounding in weekly/daily conversion ±3 days difference Document. For exact accumulation, use calendar-aware approach (future)
Sample size for daily fitting 366 distributions × 30yr = 30 samples each Warn if calibration < 20 years. Recommend dekadal/pentadal over daily for short records
Memory for daily global 366 parameter grids (~30× monthly) Chunked processing handles this. Warn in docs
SPEI data alignment PET/precip at different frequencies Aggregation helpers solve this. Validate matching time coordinates
Backward compatibility Existing monthly calls must not break Default periodicity='auto' detects monthly. All existing tests must pass

Calendar-Aware Accumulation (Future Enhancement)

The round(scale × steps_per_month) approach is simple but approximate for weekly/daily. A more precise future approach:

def calendar_aware_steps(time_coord, idx, scale_months):
    """Count actual time steps in the preceding scale_months calendar months."""
    end_date = time_coord[idx]
    start_date = end_date - pd.DateOffset(months=scale_months)
    return np.sum(time_coord[:idx+1] >= start_date)

This gives a variable-length accumulation window - exactly 12 calendar months at every time step. More accurate but significantly more complex (each time step has a different window size, breaking the fixed-window rolling sum). Recommend deferring unless users report issues with the fixed-window approach.

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or request

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions