Run Theory for Climate Extremes

Author

Benny Istanto

Published

January 31, 2026

Overview

This tutorial demonstrates run theory (Yevjevich, 1967) for identifying and analyzing climate extreme events — both droughts and wet extremes — using real TerraClimate data from Bali, Indonesia (1958-2024).

What is Run Theory?

Run theory provides a systematic approach to analyzing climate extremes:

  • Events are identified when SPI/SPEI crosses a threshold
  • Duration (D): How long the event lasts (months)
  • Magnitude (M): Cumulative deficit/surplus below/above threshold
  • Intensity (I): Average severity = M / D
  • Peak (P): Most extreme value during event
  • Inter-arrival (T): Time between successive events

This goes beyond simple threshold exceedance by capturing the full evolution of events.

What you’ll learn:

  1. Load previously calculated SPI data
  2. Visualize SPI with threshold lines
  3. Identify drought events using run theory
  4. Analyze event characteristics (duration, magnitude, intensity, peak)
  5. Identify and compare wet events
  6. Time-series monitoring (month-by-month evolution)
  7. Understand dual magnitude (cumulative vs instantaneous)
  8. Calculate gridded period statistics for spatial analysis
  9. Export events to CSV and NetCDF

Three Modes of Analysis:

Mode Function Use Case Output
Event-Based identify_events() Historical analysis DataFrame of events
Time-Series calculate_timeseries() Real-time monitoring DataFrame by month
Period Stats calculate_period_statistics() Decision support Gridded statistics

1. Setup and Imports

Code
import sys
import os
from pathlib import Path

# Get the notebook directory and find the repository root
notebook_dir = Path.cwd()
repo_root = notebook_dir
while not (repo_root / 'src').exists() and repo_root != repo_root.parent:
    repo_root = repo_root.parent

# Add src to path
src_path = str(repo_root / 'src')
if src_path not in sys.path:
    sys.path.insert(0, src_path)

import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime

# Import event analysis functions
from runtheory import (
    identify_events,
    calculate_timeseries,
    calculate_period_statistics,
    summarize_events
)

# Import visualization functions
from visualization import (
    plot_index,
    plot_events,
    plot_event_timeline
)

# Plotting settings
plt.rcParams['figure.figsize'] = (14, 6)
plt.rcParams['font.size'] = 10

# Create output directories
output_plots = repo_root / 'output' / 'plots'
output_netcdf = repo_root / 'output' / 'netcdf'
output_csv = repo_root / 'output' / 'csv'
output_plots.mkdir(parents=True, exist_ok=True)
output_netcdf.mkdir(parents=True, exist_ok=True)
output_csv.mkdir(parents=True, exist_ok=True)

print("All imports successful!")
print(f"Repository root: {repo_root}")
All imports successful!
Repository root: /home/runner/work/precip-index/precip-index

2. Load SPI-12 Data

Load previously calculated SPI-12 data from Tutorial 01.

Code
# Load SPI-12
spi_file = repo_root / 'output' / 'netcdf' / 'spi_12_bali.nc'

ds = xr.open_dataset(spi_file)
spi_12 = ds['spi_gamma_12_month']

print("SPI-12 loaded successfully!")
print(f"\nDataset Information:")
print(f"  Shape: {spi_12.shape}")
print(f"  Dimensions: {spi_12.dims}")
print(f"  Time range: {spi_12.time[0].values} to {spi_12.time[-1].values}")
print(f"  Spatial extent: {len(spi_12.lat)} x {len(spi_12.lon)} grid (Bali, Indonesia)")
print(f"  Lat range: {float(spi_12.lat.min()):.2f} to {float(spi_12.lat.max()):.2f}")
print(f"  Lon range: {float(spi_12.lon.min()):.2f} to {float(spi_12.lon.max()):.2f}")
SPI-12 loaded successfully!

Dataset Information:
  Shape: (804, 24, 35)
  Dimensions: ('time', 'lat', 'lon')
  Time range: 1958-01-01T00:00:00.000000000 to 2024-12-01T00:00:00.000000000
  Spatial extent: 24 x 35 grid (Bali, Indonesia)
  Lat range: -8.94 to -7.98
  Lon range: 114.35 to 115.77

3. Select Sample Location

For event-based and time-series analysis, extract a single location in central Bali.

Code
# Select center of Bali
lat_idx = len(spi_12.lat) // 2
lon_idx = len(spi_12.lon) // 2

# Extract location
spi_loc = spi_12.isel(lat=lat_idx, lon=lon_idx)
lat_val = float(spi_12.lat.values[lat_idx])
lon_val = float(spi_12.lon.values[lon_idx])

print(f"Selected location: {lat_val:.2f}, {lon_val:.2f} (Central Bali)")
print(f"SPI time series length: {len(spi_loc)} months ({len(spi_loc)/12:.1f} years)")
print(f"Mean SPI: {float(spi_loc.mean()):.3f}")
print(f"Std SPI: {float(spi_loc.std()):.3f}")
Selected location: -8.48, 115.06 (Central Bali)
SPI time series length: 804 months (67.0 years)
Mean SPI: -0.049
Std SPI: 1.006

Visualize SPI-12 Time Series

Code
fig, ax = plt.subplots(figsize=(14, 5))
spi_loc.plot(ax=ax, linewidth=0.8, color='steelblue')
ax.axhline(y=0, color='k', linestyle='-', linewidth=0.8, alpha=0.3)
ax.axhline(y=-1.2, color='red', linestyle='--', linewidth=0.8, alpha=0.5, label='Dry threshold -1.2')
ax.axhline(y=1.2, color='blue', linestyle='--', linewidth=0.8, alpha=0.5, label='Wet threshold +1.2')
ax.fill_between(spi_loc.time, -5, 0, alpha=0.1, color='red', label='Dry')
ax.fill_between(spi_loc.time, 0, 5, alpha=0.1, color='blue', label='Wet')
ax.set_ylim(-3, 3)
ax.set_title(f'SPI-12 at {lat_val:.2f}, {lon_val:.2f} (Central Bali, 1958-2024)')
ax.set_ylabel('SPI-12')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

The dashed lines at +-1.2 mark the thresholds for identifying moderate drought and wet events. Periods below -1.2 are drought events; periods above +1.2 are wet events.

4. Mode 1: Event-Based Analysis (Dry Events)

Identify complete dry events using run theory. Each event is characterized by duration, magnitude, intensity, peak, and inter-arrival time.

Code
# Identify dry events (drought)
threshold_dry = -1.2  # Moderate dry threshold
min_duration = 3  # Minimum 3 months to be considered an event

print(f"Identifying dry events with:")
print(f"  Threshold: {threshold_dry}")
print(f"  Minimum duration: {min_duration} months")
print()

events_dry = identify_events(spi_loc, threshold=threshold_dry, min_duration=min_duration)

print(f"Found {len(events_dry)} dry events")
print()
print("Dry Event Summary:")
print(events_dry[['start_date', 'end_date', 'duration', 'magnitude', 'intensity', 'peak']])
Identifying dry events with:
  Threshold: -1.2
  Minimum duration: 3 months

Found 14 dry events

Dry Event Summary:
   start_date   end_date  duration  magnitude  intensity      peak
0  1961-09-01 1962-07-01        11   4.158311   0.378028 -1.976741
1  1965-10-01 1966-01-01         4   1.563179   0.390795 -1.883761
2  1972-10-01 1972-12-01         3   1.060205   0.353402 -1.761785
3  1976-11-01 1977-12-01        14   3.015663   0.215405 -1.726251
4  1980-06-01 1980-11-01         6   0.842345   0.140391 -1.538919
5  1982-11-01 1983-09-01        11   4.622329   0.420212 -2.224815
6  1987-07-01 1987-11-01         5   2.381117   0.476223 -1.795504
7  1988-05-01 1988-09-01         5   0.601202   0.120240 -1.448271
8  1991-08-01 1991-10-01         3   0.781494   0.260498 -1.601405
9  2002-06-01 2002-12-01         7   2.373209   0.339030 -1.967303
10 2005-05-01 2005-07-01         3   0.212718   0.070906 -1.389078
11 2014-06-01 2014-12-01         7   3.313735   0.473391 -1.905663
12 2023-12-01 2024-02-01         3   0.221492   0.073831 -1.351519
13 2024-04-01 2024-09-01         6   1.324795   0.220799 -1.572638
NoteThreshold Selection

Common thresholds for SPI/SPEI event analysis:

  • -0.7: Mild drought (more events, less severe)
  • -1.0: Moderate drought (WMO recommendation)
  • -1.2: Significant drought (used in this example)
  • -1.5: Severe drought (fewer but more intense events)

Dry Event Statistics

Code
# Summarize dry events
summary_dry = summarize_events(events_dry)

print("Dry Event Statistics:")
print("=" * 50)
print(f"Total events: {summary_dry['num_events']}")
print(f"Mean duration: {summary_dry['mean_duration']:.1f} months")
print(f"Max duration: {summary_dry['max_duration']} months")
print(f"Mean magnitude: {summary_dry['mean_magnitude']:.2f}")
print(f"Max magnitude: {summary_dry['max_magnitude']:.2f}")
print(f"Most severe peak: {summary_dry['most_severe_peak']:.2f}")
print(f"Mean inter-arrival: {summary_dry['mean_interarrival']:.1f} months")
Dry Event Statistics:
==================================================
Total events: 14.0
Mean duration: 6.3 months
Max duration: 14.0 months
Mean magnitude: 1.89
Max magnitude: 4.62
Most severe peak: -2.22
Mean inter-arrival: 57.8 months

Visualize Dry Events

Code
fig = plot_events(spi_loc, events_dry, threshold=threshold_dry,
                  title=f'Dry Events at {lat_val:.2f}, {lon_val:.2f} (Bali, 1958-2024)')
plt.tight_layout()
plt.show()

Drought events are highlighted in the time series. Each shaded region represents a complete event that exceeded the threshold for at least 3 consecutive months.

5. Wet Events Analysis

The same functions work for wet events — just use a positive threshold.

Code
# Identify wet events
threshold_wet = 1.2  # Moderate wet threshold (positive!)

print(f"Identifying wet events with:")
print(f"  Threshold: {threshold_wet}")
print(f"  Minimum duration: {min_duration} months")
print()

events_wet = identify_events(spi_loc, threshold=threshold_wet, min_duration=min_duration)

print(f"Found {len(events_wet)} wet events")
print()
print("Wet Event Summary:")
print(events_wet[['start_date', 'end_date', 'duration', 'magnitude', 'intensity', 'peak']])
Identifying wet events with:
  Threshold: 1.2
  Minimum duration: 3 months

Found 8 wet events

Wet Event Summary:
  start_date   end_date  duration  magnitude  intensity      peak
0 1968-09-01 1969-04-01         8   2.090736   0.261342  1.618894
1 1973-10-01 1974-06-01         9   2.115463   0.235051  1.619856
2 1975-04-01 1976-04-01        13   7.253528   0.557964  1.994779
3 1981-11-01 1982-06-01         8   1.791800   0.223975  1.631948
4 1998-09-01 1999-09-01        13  11.949839   0.919218  2.633997
5 2010-09-01 2011-08-01        12  12.105968   1.008831  2.445220
6 2016-10-01 2017-09-01        12   9.510695   0.792558  2.291671
7 2022-10-01 2023-02-01         5   0.573008   0.114602  1.535313

Wet Event Statistics

Code
summary_wet = summarize_events(events_wet)

print("Wet Event Statistics:")
print("=" * 50)
print(f"Total events: {summary_wet['num_events']}")
print(f"Mean duration: {summary_wet['mean_duration']:.1f} months")
print(f"Max duration: {summary_wet['max_duration']} months")
print(f"Mean magnitude: {summary_wet['mean_magnitude']:.2f}")
print(f"Max magnitude: {summary_wet['max_magnitude']:.2f}")
print(f"Most extreme peak: {summary_wet['most_severe_peak']:.2f}")
print(f"Mean inter-arrival: {summary_wet['mean_interarrival']:.1f} months")
Wet Event Statistics:
==================================================
Total events: 8.0
Mean duration: 10.0 months
Max duration: 13.0 months
Mean magnitude: 5.92
Max magnitude: 12.11
Most extreme peak: 1.54
Mean inter-arrival: 92.7 months

Visualize Wet Events

Code
fig = plot_events(spi_loc, events_wet, threshold=threshold_wet,
                  title=f'Wet Events at {lat_val:.2f}, {lon_val:.2f} (Bali, 1958-2024)')
plt.tight_layout()
plt.show()

Compare Dry vs Wet Events

Code
# Comparison summary
print("=" * 50)
print("Comparison: Dry vs Wet Events")
print("=" * 50)
print(f"Dry events: {summary_dry['num_events']:.0f}  |  Wet events: {summary_wet['num_events']:.0f}")
print(f"Dry avg duration: {summary_dry['mean_duration']:.1f} months  |  Wet avg duration: {summary_wet['mean_duration']:.1f} months")
print(f"Dry avg magnitude: {summary_dry['mean_magnitude']:.2f}  |  Wet avg magnitude: {summary_wet['mean_magnitude']:.2f}")

# Visualize comparison
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Duration comparison
ax1.hist(events_dry['duration'], bins=15, alpha=0.7, label='Drought', color='orange')
ax1.hist(events_wet['duration'], bins=15, alpha=0.7, label='Wet', color='blue')
ax1.set_xlabel('Duration (months)')
ax1.set_ylabel('Frequency')
ax1.set_title('Event Duration Distribution')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Magnitude comparison
ax2.hist(events_dry['magnitude'], bins=15, alpha=0.7, label='Drought', color='orange')
ax2.hist(events_wet['magnitude'], bins=15, alpha=0.7, label='Wet', color='blue')
ax2.set_xlabel('Magnitude')
ax2.set_ylabel('Frequency')
ax2.set_title('Event Magnitude Distribution')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
==================================================
Comparison: Dry vs Wet Events
==================================================
Dry events: 14  |  Wet events: 8
Dry avg duration: 6.3 months  |  Wet avg duration: 10.0 months
Dry avg magnitude: 1.89  |  Wet avg magnitude: 5.92

ImportantBidirectional Analysis

The same functions work for both extremes:

  • Droughts: negative thresholds (e.g., -1.2)
  • Wet events: positive thresholds (e.g., +1.2)

This unified framework makes analysis consistent and reproducible.

6. Mode 2: Time-Series Monitoring

Calculate month-by-month characteristics for real-time monitoring. This mode tracks the evolution of events as they unfold — useful for operational systems.

Code
# Calculate dry event time series
print("Calculating dry event time series...")
ts_dry = calculate_timeseries(spi_loc, threshold=threshold_dry)

print(f"Time series calculated")
print(f"  Length: {len(ts_dry)} months")
print(f"  Columns: {list(ts_dry.columns)}")
Calculating dry event time series...
Time series calculated
  Length: 804 months
  Columns: ['index_value', 'is_event', 'event_id', 'duration', 'magnitude_cumulative', 'magnitude_instantaneous', 'intensity', 'peak_so_far', 'deviation']

Current Status Check

Code
# Check current status (last month in data)
current = ts_dry.iloc[-1]

print("Current Status (Latest Data):")
print("=" * 50)
print(f"Date: {ts_dry.index[-1]}")
print(f"SPI-12: {current['index_value']:.2f}")

if current['is_event']:
    print(f"\nDRY EVENT IN PROGRESS")
    print(f"  Event ID: {current['event_id']}")
    print(f"  Duration: {current['duration']} months")
    print(f"  Cumulative magnitude: {current['magnitude_cumulative']:.2f}")
    print(f"  Current severity: {current['magnitude_instantaneous']:.2f}")
    print(f"  Intensity: {current['intensity']:.2f}")
else:
    print(f"\nNO DRY EVENT - Normal conditions")
Current Status (Latest Data):
==================================================
Date: 2024-12-01 00:00:00
SPI-12: -0.29

NO DRY EVENT - Normal conditions

Event Evolution Timeline (5-Panel Plot)

The 5-panel event timeline shows:

  1. Index value (SPI-12)
  2. Duration (months)
  3. Cumulative magnitude (total deficit, always increasing during event)
  4. Instantaneous magnitude (current severity, varies within event)
  5. Intensity (average severity)
Code
fig = plot_event_timeline(ts_dry)

# Customize title
axes = fig.get_axes()
axes[0].set_title('')
fig.suptitle(f'Dry Event Evolution at {lat_val:.2f}, {lon_val:.2f} (Bali)',
             fontsize=14, fontweight='bold', y=0.995)
plt.subplots_adjust(top=0.97)
plt.tight_layout()
plt.show()

Magnitude Comparison: Cumulative vs Instantaneous

Two types of magnitude provide different insights:

  • Cumulative (blue): Total deficit, monotonically increases during event — like debt accumulation
  • Instantaneous (red): Current severity, varies with SPI pattern — like NDVI crop phenology
Code
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 8), sharex=True)

# Get event mask
event_mask = ts_dry['is_event']

# Cumulative magnitude
cumulative_event = ts_dry['magnitude_cumulative'].where(event_mask)
ax1.fill_between(ts_dry.index, 0, cumulative_event,
                 alpha=0.4, color='steelblue', label='Cumulative')
ax1.plot(ts_dry.index, cumulative_event,
         color='darkblue', linewidth=1.5)
ax1.set_ylabel('Cumulative Magnitude', fontsize=11)
ax1.set_title('Cumulative Magnitude (Total Deficit - Always Increasing)', fontsize=12)
ax1.set_ylim(bottom=0)
ax1.legend()
ax1.grid(True, alpha=0.3)

# Instantaneous magnitude
instantaneous_event = ts_dry['magnitude_instantaneous'].where(event_mask)
ax2.fill_between(ts_dry.index, 0, instantaneous_event,
                 alpha=0.4, color='coral', label='Instantaneous')
ax2.plot(ts_dry.index, instantaneous_event,
         color='darkred', linewidth=1.5)
ax2.set_ylabel('Instantaneous Magnitude', fontsize=11)
ax2.set_xlabel('Time', fontsize=11)
ax2.set_title('Instantaneous Magnitude (Current Severity - Varies Within Event)', fontsize=12)
ax2.set_ylim(bottom=0)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.suptitle('Dual Magnitude Comparison - Dry Events (Bali)', fontsize=14, y=0.995)
plt.tight_layout()
plt.show()

The cumulative magnitude grows throughout each event, showing total accumulated deficit. The instantaneous magnitude shows current severity at each timestep — it can increase or decrease as the drought intensifies or eases within a single event.

7. Mode 3: Period Statistics (Gridded)

Calculate spatial statistics for specific time periods across all grid cells. This answers decision-maker questions like:

  • “How many drought events occurred in 2020?”
  • “Where was the worst drought in the last 5 years?”
  • “What percentage of time was each area in drought?”

Single-Year Statistics (2020)

Code
print("Calculating dry event statistics for 2020...")
stats_2020 = calculate_period_statistics(spi_12, threshold=threshold_dry,
                                          start_year=2020, end_year=2020,
                                          min_duration=min_duration)

print(f"Statistics calculated")
print(f"  Variables: {list(stats_2020.data_vars)}")
Calculating dry event statistics for 2020...
Statistics calculated
  Variables: ['num_events', 'total_event_months', 'total_magnitude', 'mean_magnitude', 'max_magnitude', 'worst_peak', 'mean_intensity', 'max_intensity', 'pct_time_in_event']
Code
# Plot number of events and worst severity
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

stats_2020['num_events'].plot(ax=ax1, cmap='YlOrRd',
                               cbar_kwargs={'label': 'Event Count', 'format': '%d'})
ax1.set_title('Number of Dry Events in 2020 - Bali', fontsize=12)
ax1.set_xlabel('Longitude')
ax1.set_ylabel('Latitude')

stats_2020['worst_peak'].plot(ax=ax2, cmap='Reds_r',
                               cbar_kwargs={'label': 'Worst SPI-12 Value'})
ax2.set_title('Worst Dry Event Severity in 2020 - Bali', fontsize=12)
ax2.set_xlabel('Longitude')
ax2.set_ylabel('Latitude')

plt.tight_layout()
plt.show()

Multi-Year Statistics (2020–2024)

Code
print("Calculating dry event statistics for 2020-2024...")
stats_5yr = calculate_period_statistics(spi_12, threshold=threshold_dry,
                                         start_year=2020, end_year=2024,
                                         min_duration=min_duration)

print(f"5-year statistics calculated")
print(f"  Variables: {list(stats_5yr.data_vars)}")
Calculating dry event statistics for 2020-2024...
5-year statistics calculated
  Variables: ['num_events', 'total_event_months', 'total_magnitude', 'mean_magnitude', 'max_magnitude', 'worst_peak', 'mean_intensity', 'max_intensity', 'pct_time_in_event']
Code
# 2x2 panel of key statistics
fig, axes = plt.subplots(2, 2, figsize=(16, 12))

variables = ['num_events', 'worst_peak', 'total_magnitude', 'pct_time_in_event']
titles = ['Event Count', 'Worst Severity (Most Negative SPI)', 'Total Magnitude', '% Time in Dry Event']
cmaps = ['YlOrRd', 'Reds_r', 'YlOrRd', 'Reds']

for ax, var, title, cmap in zip(axes.flat, variables, titles, cmaps):
    if var == 'num_events':
        cbar_kwargs = {'label': title, 'format': '%d'}
    elif var == 'pct_time_in_event':
        cbar_kwargs = {'label': title, 'format': '%.0f'}
    else:
        cbar_kwargs = {'label': title}

    stats_5yr[var].plot(ax=ax, cmap=cmap, cbar_kwargs=cbar_kwargs)
    ax.set_title(title, fontsize=12)
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')

plt.suptitle('Dry Event Statistics Summary - Bali (2020-2024)', fontsize=14, y=0.995)
plt.tight_layout()
plt.show()

The four-panel summary shows: (1) how many drought events occurred, (2) the worst severity reached, (3) the total cumulative magnitude across all events, and (4) the percentage of time each grid cell spent in drought. Areas with higher values in all panels are the most drought-prone.

8. Export Results

Save events to CSV and statistics to NetCDF for further analysis.

Code
# Save drought events
events_dry.to_csv(output_csv / 'drought_events_bali.csv', index=False)
print(f"Drought events saved to: {output_csv / 'drought_events_bali.csv'}")

# Save wet events
events_wet.to_csv(output_csv / 'wet_events_bali.csv', index=False)
print(f"Wet events saved to: {output_csv / 'wet_events_bali.csv'}")

# Save time series
ts_dry.to_csv(output_csv / 'drought_timeseries_bali.csv')
print(f"Time series saved to: {output_csv / 'drought_timeseries_bali.csv'}")

# Save gridded statistics
stats_5yr.to_netcdf(str(output_netcdf / 'dry_stats_2020-2024_bali.nc'))
print(f"Gridded statistics saved to: {output_netcdf / 'dry_stats_2020-2024_bali.nc'}")
Drought events saved to: /home/runner/work/precip-index/precip-index/output/csv/drought_events_bali.csv
Wet events saved to: /home/runner/work/precip-index/precip-index/output/csv/wet_events_bali.csv
Time series saved to: /home/runner/work/precip-index/precip-index/output/csv/drought_timeseries_bali.csv
Gridded statistics saved to: /home/runner/work/precip-index/precip-index/output/netcdf/dry_stats_2020-2024_bali.nc

9. Summary

Key Takeaways

  1. Run theory provides a systematic framework for identifying complete extreme events with duration, magnitude, intensity, and peak characteristics
  2. Bidirectional analysis: The same functions work for both drought (negative threshold) and wet events (positive threshold)
  3. Three analysis modes serve different purposes: event-based (historical), time-series (monitoring), and period statistics (decision support)
  4. Dual magnitude: Cumulative tracks total accumulated deficit; instantaneous tracks current severity
  5. Gridded statistics enable spatial analysis of drought-prone areas

Best Practices

Parameter Recommendation
Threshold +-1.0 or +-1.2 for operational monitoring
Min duration 3 months for SPI-12 (captures sustained events)
Period statistics Use specific time windows for targeted questions
Both extremes Analyze both dry and wet for a complete picture

Next Steps

References

  • Yevjevich, V. (1967): An objective approach to definitions and investigations of continental hydrologic droughts
  • McKee et al. (1993): SPI methodology
Back to top