Memantau Kekeringan dengan Indeks SPI dan SPEI

Musim kering 2023 mencatat salah satu episode kekeringan terburuk Indonesia dalam beberapa tahun terakhir. El Niño mendorong anomali SST Niño 3.4 ke +1,42°C pada Agustus 2023, bertepatan dengan IOD positif +0,98°C. Akibatnya 23.451 hektare sawah mengalami kekeringan dan 6.964 hektare gagal panen total. Di Kalimantan, lahan gambut yang mengering memicu kebakaran yang membakar sekitar 267.900 hektare pada awal Oktober 2023 — melampaui total area yang terbakar sepanjang 2022.

Dua indeks terstandarisasi — Standardized Precipitation Index (SPI) dan ekstensinya SPEI — menjadi alat standar WMO untuk mendeteksi kondisi seperti itu jauh sebelum sawah layu. Tutorial ini membangun keduanya dari awal: download ERA5 via cdsapi, hitung SPI dan SPEI dengan scipy, visualisasikan hasilnya. BMKG sendiri mengoperasikan pemantauan SPI nasional dalam Climate Early Warning System (CEWS) — fondasi yang sama yang kita pasang di sini.

Sumber: NASA Earth Observatory — Asap tebal dari kebakaran di Kalimantan Selatan dan Tengah, 2 Oktober 2023, akibat kekeringan El Niño 2023. (artikel asli)

Implementasi SPI dengan Python dan ERA5

SPI dikembangkan McKee, Doesken, dan Kleist (1993). Curah hujan mentah mengikuti distribusi gamma — miring ke kanan, tidak bisa langsung dibandingkan antar wilayah. Solusinya: fit gamma per kalender bulan, lalu transformasi ke standar normal via inverse CDF:

$$\text{SPI} = \Phi^{-1}(F_\gamma(x))$$

Hasilnya z-score yang diinterpretasikan via skala McKee:

SPI	Kategori
\(\leq -2{,}0\)	Kekeringan Ekstrem
\(-1{,}99\) s.d. \(-1{,}50\)	Kekeringan Parah
\(-1{,}49\) s.d. \(-1{,}00\)	Kekeringan Sedang
\(-0{,}99\) s.d. \(+0{,}99\)	Normal
\(\geq +1{,}00\)	Basah

Karena skala simetris dan distribusi standar normal universal, satu nilai SPI-3 = −1,7 di Jakarta dapat dibandingkan langsung dengan SPI-3 = −1,7 di Australia — sesuatu yang mustahil dilakukan dengan defisit milimeter mentah. Pilihan timescale menentukan jenis kekeringan yang terdeteksi: SPI-1 untuk stres tanaman jangka pendek, SPI-3 untuk kelembaban tanah musiman, SPI-6 untuk debit sungai, SPI-12 untuk pengisian air tanah. WMO merekomendasikan menghitung beberapa timescale dari dataset yang sama agar berbagai jenis kekeringan tertangkap secara paralel.

Setup dan fungsi pembantu

import os, warnings
import numpy as np
import xarray as xr
from scipy import stats
import pandas as pd

warnings.filterwarnings("ignore")

def fit_and_transform_spi(series: np.ndarray) -> np.ndarray:
    """Fit gamma, return SPI array. Zeros get p=0."""
    spi = np.full(len(series), np.nan)
    valid = series[series > 0]
    if len(valid) < 4:
        return spi
    alpha, loc, beta = stats.gamma.fit(valid, floc=0)
    p = np.where(series <= 0, 0.0,
                 stats.gamma.cdf(series, alpha, loc=loc, scale=beta))
    p = np.clip(p, 1e-6, 1 - 1e-6)
    spi[:] = stats.norm.ppf(p)
    return spi

def classify_spi(val):
    if np.isnan(val): return "N/A"
    if val <= -2.0:   return "Ekstrem"
    if val <= -1.5:   return "Parah"
    if val <= -1.0:   return "Sedang"
    if val >= 1.0:    return "Basah"
    return "Normal"

# Titik grid Jakarta
LAT, LON = -6.2, 106.8
print("Library dimuat. Fungsi helper dan konstanta titik Jakarta siap.")

Library dimuat. Fungsi helper dan konstanta titik Jakarta siap.

Download ERA5 presipitasi dan agregasi bulanan

ERA5 total_precipitation adalah akumulasi 1-jam pada setiap timestamp. Sampling 6-jam (00, 06, 12, 18 UTC) menangkap satu jam dari setiap blok enam-jam, sehingga total bulanan yang kita peroleh hanya proporsional — bagus untuk mempelajari pipeline SPI tetapi bukan angka klimatologis yang lengkap. Untuk produksi, gunakan data hourly (atau produk monthly-means ERA5).

import os, cdsapi, xarray as xr

OUT_TP = "era5_tp_indonesia_2024_6h.nc"
if not os.path.exists(OUT_TP):
    c = cdsapi.Client(quiet=True)
    c.retrieve(
        "reanalysis-era5-single-levels",
        {
            "product_type": "reanalysis",
            "variable": ["total_precipitation"],
            "year": "2024",
            "month": [f"{m:02d}" for m in range(1, 13)],
            "day":   [f"{d:02d}" for d in range(1, 32)],
            "time":  ["00:00", "06:00", "12:00", "18:00"],
            "area":  [6, 95, -11, 141],   # N, W, S, E — Indonesia
            "format": "netcdf",
        },
        OUT_TP,
    )

ds_tp = xr.open_dataset(OUT_TP)
if "valid_time" in ds_tp.dims:
    ds_tp = ds_tp.rename({"valid_time": "time"})

tp_var = "tp" if "tp" in ds_tp else list(ds_tp.data_vars)[0]
tp_monthly = ds_tp[tp_var].resample(time="1ME").sum() * 1000
tp_monthly.attrs["units"] = "mm/month"

print(f"Dimensi  : {dict(tp_monthly.sizes)}")
print(f"Periode  : {str(tp_monthly.time.values[0])[:7]} → {str(tp_monthly.time.values[-1])[:7]}")
print(f"Rentang  : {float(tp_monthly.min()):.1f} – {float(tp_monthly.max()):.1f} mm/bulan")

Dimensi  : {'time': 12, 'latitude': 69, 'longitude': 185}
Periode  : 2024-01 → 2024-12
Rentang  : 0.0 – 484.7 mm/bulan

Kalkulasi SPI-3 untuk wilayah Jakarta

Kita extract deret waktu bulanan di titik Jakarta, terapkan rolling 3-bulan, lalu fit gamma. Karena data hanya mencakup 2024 (12 bulan), kita gunakan pooled fit atas semua bulan — cara yang valid secara pedagogis. Produksi SPI yang benar membutuhkan baseline minimal 30–50 tahun sesuai panduan WMO.

import numpy as np, pandas as pd
from scipy import stats

tp_jkt = tp_monthly.sel(latitude=LAT, longitude=LON, method="nearest")
tp_arr  = tp_jkt.values.astype(float)

# Rolling 3-bulan (leading 2 bulan = NaN)
tp_roll3 = np.array([
    np.sum(tp_arr[max(0, i-2):i+1]) if i >= 2 else np.nan
    for i in range(len(tp_arr))
])

valid_vals = tp_roll3[~np.isnan(tp_roll3)]
spi3_valid = fit_and_transform_spi(valid_vals)

spi3_full = np.full(len(tp_roll3), np.nan)
spi3_full[2:] = spi3_valid

print("=== SPI-3 Jakarta 2024 (DEMO — baseline hanya 12 bln, bukan 30 thn) ===")
print(f"{'Bulan':<10}{'P-3bln (mm)':<16}{'SPI-3':<10}{'Kategori'}")
print("-" * 48)
months = pd.date_range("2024-01", periods=12, freq="ME")
for m, p, s in zip(months, tp_roll3, spi3_full):
    s_str = f"{s:+.3f}" if not np.isnan(s) else "   — "
    p_str = f"{p:.1f}"  if not np.isnan(p) else " — "
    print(f"{str(m)[:7]:<10}{p_str:<16}{s_str:<10}{classify_spi(s)}")

=== SPI-3 Jakarta 2024 (DEMO — baseline hanya 12 bln, bukan 30 thn) ===
Bulan     P-3bln (mm)     SPI-3     Kategori
------------------------------------------------
2024-01    —                 —      N/A
2024-02    —                 —      N/A
2024-03   130.4           +1.517    Basah
2024-04   112.9           +1.210    Basah
2024-05   80.8            +0.559    Normal
2024-06   53.6            -0.147    Normal
2024-07   42.2            -0.516    Normal
2024-08   29.0            -1.044    Sedang
2024-09   27.0            -1.138    Sedang
2024-10   20.2            -1.498    Sedang
2024-11   60.3            +0.048    Normal
2024-12   101.6           +0.997    Normal

Output menampilkan klasifikasi SPI-3 per bulan. Bulan-bulan musim kering (Juli–Oktober) umumnya menunjukkan SPI negatif — sinyal defisit presipitasi yang diperparah El Niño.

Visualisasi time series SPI-3

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.dates as mdates
import numpy as np, pandas as pd

months_plot = pd.date_range("2024-01", periods=12, freq="ME")

def bar_color(s):
    if np.isnan(s):  return "#aaaaaa"
    if s <= -2.0:    return "#b22222"
    if s <= -1.5:    return "#e87020"
    if s <= -1.0:    return "#f5c518"
    if s >= 1.0:     return "#2e7d32"
    return "#90caf9"

colors = [bar_color(s) for s in spi3_full]

fig, ax = plt.subplots(figsize=(11, 5))
ax.bar(months_plot, spi3_full, width=20, color=colors, edgecolor="white", linewidth=0.5)
ax.axhline(-1.0, color="#f5c518", lw=1.2, ls="--")
ax.axhline(-1.5, color="#e87020", lw=1.2, ls="--")
ax.axhline(-2.0, color="#b22222", lw=1.4, ls="--")
ax.axhline(0,    color="black",   lw=0.8)

patches = [
    mpatches.Patch(color="#b22222", label="Ekstrem ≤ −2,0"),
    mpatches.Patch(color="#e87020", label="Parah −1,5…−2,0"),
    mpatches.Patch(color="#f5c518", label="Sedang −1,0…−1,5"),
    mpatches.Patch(color="#90caf9", label="Normal"),
    mpatches.Patch(color="#2e7d32", label="Basah ≥ +1,0"),
]
ax.legend(handles=patches, loc="upper right", fontsize=8, framealpha=0.8)
ax.set_ylabel("SPI-3", fontsize=12)
ax.set_xlabel("Bulan", fontsize=12)
ax.set_title("SPI-3 Jakarta (lat ≈ −6,2°, lon ≈ 106,8°) — ERA5 2024", fontsize=11)
ax.set_ylim(-3.5, 3.5)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax.xaxis.set_major_locator(mdates.MonthLocator())
plt.tight_layout()
plt.savefig("snippet-4.png", dpi=120, bbox_inches="tight")
print("saved")

Kolom merah atau oranye yang muncul di pertengahan–akhir tahun mengonfirmasi defisit presipitasi di musim kering — persis periode saat El Niño menggeser rainfall belt dari Indonesia ke Samudra Pasifik tengah.

Konsep dan Implementasi SPEI

SPI hanya melihat supply air (curah hujan), tetapi mengabaikan demand atmosfer. Ketika suhu naik, evapotranspirasi potensial ikut naik — artinya bahkan dengan curah hujan yang sama, kekeringan bisa terasa lebih parah karena air menguap lebih cepat.

SPEI, yang diperkenalkan Vicente-Serrano et al. (2010), menggunakan water balance:

$$D = P - \text{PET}$$

di mana \(D\) negatif berarti demand lebih besar dari supply. Karena \(D\) bisa negatif, distribusi gamma tidak cocok. SPEI memakai distribusi log-logistic tiga parameter (scipy.stats.fisk), yang telah terbukti cocok untuk deret \(D\) di semua timescale dan iklim.

Tutorial ini menggunakan Thornthwaite PET — hanya butuh suhu bulanan dan lintang:

$$\text{PET} = 16 \left(\frac{10\,T}{I}\right)^a \cdot C$$

di mana \(I = \sum_{i=1}^{12} (T_i/5)^{1{,}514}\) adalah heat index tahunan, \(a = 6{,}75 \times 10^{-7} I^3 - 7{,}71 \times 10^{-5} I^2 + 1{,}792 \times 10^{-2} I + 0{,}49239\), dan \(C\) adalah faktor koreksi panjang hari.

Download ERA5 temperature dan hitung PET Thornthwaite

import os, cdsapi, xarray as xr, numpy as np

OUT_T2M = "era5_t2m_indonesia_2024_6h.nc"
if not os.path.exists(OUT_T2M):
    c = cdsapi.Client(quiet=True)
    c.retrieve(
        "reanalysis-era5-single-levels",
        {
            "product_type": "reanalysis",
            "variable": ["2m_temperature"],
            "year": "2024",
            "month": [f"{m:02d}" for m in range(1, 13)],
            "day":   [f"{d:02d}" for d in range(1, 32)],
            "time":  ["00:00", "06:00", "12:00", "18:00"],
            "area":  [6, 95, -11, 141],
            "format": "netcdf",
        },
        OUT_T2M,
    )

ds_t = xr.open_dataset(OUT_T2M)
if "valid_time" in ds_t.dims:
    ds_t = ds_t.rename({"valid_time": "time"})
t_var = "t2m" if "t2m" in ds_t else list(ds_t.data_vars)[0]
t_monthly = (ds_t[t_var].resample(time="1ME").mean() - 273.15)
t_jkt = t_monthly.sel(latitude=LAT, longitude=LON, method="nearest").values

def thornthwaite_pet(T_monthly, latitude_deg):
    T = np.maximum(T_monthly, 0.0)
    I = np.sum((T / 5) ** 1.514)
    a = 6.75e-7 * I**3 - 7.71e-5 * I**2 + 1.792e-2 * I + 0.49239
    lat_rad = np.deg2rad(abs(latitude_deg))
    decl_deg = np.array([-20.0, -10.8, 0.0, 11.6, 20.0, 23.45,
                          23.45, 20.0, 11.6, 0.0, -10.8, -20.0])
    days_pm  = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
    photoperiod = np.array([
        (24 / np.pi) * np.arccos(
            np.clip(-np.tan(lat_rad) * np.tan(np.deg2rad(d)), -1, 1)
        )
        for d in decl_deg
    ])
    C = (photoperiod / 12.0) * (days_pm / 30.0)
    return np.where(T > 0, 16.0 * ((10.0 * T / I) ** a) * C, 0.0)

pet_jkt = thornthwaite_pet(t_jkt, LAT)

import pandas as pd
print("=== Suhu bulanan dan PET Thornthwaite — Jakarta 2024 ===")
print(f"{'Bulan':<8}{'T (°C)':<10}{'PET (mm/bln)'}")
print("-" * 32)
months_s = pd.date_range("2024-01", periods=12, freq="ME")
for m, T, P in zip(months_s, t_jkt, pet_jkt):
    print(f"{str(m)[:7]:<8}{T:<10.1f}{P:.1f}")
print(f"\nPET tahunan: {pet_jkt.sum():.0f} mm")
print(f"P tahunan  : {float(tp_monthly.sel(latitude=LAT, longitude=LON, method='nearest').sum()):.0f} mm")

=== Suhu bulanan dan PET Thornthwaite — Jakarta 2024 ===
Bulan   T (°C)    PET (mm/bln)
--------------------------------
2024-01 27.2      146.2
2024-02 27.3      134.4
2024-03 27.2      150.0
2024-04 27.8      161.6
2024-05 28.3      182.0
2024-06 27.5      155.3
2024-07 27.0      149.9
2024-08 27.8      166.9
2024-09 28.2      171.4
2024-10 28.8      189.5
2024-11 27.7      153.2
2024-12 27.2      146.6

PET tahunan: 1907 mm
P tahunan  : 313 mm

Suhu Jakarta konsisten ~26–28°C sepanjang tahun (khas tropis), sehingga PET Thornthwaite relatif stabil. Perbandingan PET tahunan vs. presipitasi tahunan menunjukkan posisi Jakarta pada neraca air tahunan.

Kalkulasi SPEI-3 dan perbandingan dengan SPI-3

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import numpy as np, pandas as pd
from scipy import stats

# Water balance dan rolling 3-bulan
tp_jkt_arr = tp_monthly.sel(latitude=LAT, longitude=LON, method="nearest").values.astype(float)
D_monthly  = tp_jkt_arr - pet_jkt
D_roll3    = np.array([
    np.sum(D_monthly[max(0, i-2):i+1]) if i >= 2 else np.nan
    for i in range(len(D_monthly))
])

def compute_spei(series):
    spei_out = np.full(len(series), np.nan)
    valid_idx = ~np.isnan(series)
    v = series[valid_idx]
    if len(v) < 4: return spei_out
    c, loc, scale = stats.fisk.fit(v)
    p = np.clip(stats.fisk.cdf(v, c, loc=loc, scale=scale), 1e-6, 1 - 1e-6)
    spei_out[valid_idx] = stats.norm.ppf(p)
    return spei_out

spei3_full = compute_spei(D_roll3)
months_plot = pd.date_range("2024-01", periods=12, freq="ME")

fig, axes = plt.subplots(2, 1, figsize=(11, 8), sharex=True)
configs = [
    (axes[0], spi3_full,  "SPI-3",  "#1976d2", "#e53935",
     "SPI-3 Jakarta 2024 — hanya presipitasi"),
    (axes[1], spei3_full, "SPEI-3", "#388e3c", "#bf360c",
     "SPEI-3 Jakarta 2024 — P minus PET (Thornthwaite)"),
]
for ax, values, label, cpos, cneg, title in configs:
    bc = [cpos if (not np.isnan(v) and v >= 0) else cneg for v in values]
    ax.bar(months_plot, values, width=20, color=bc, edgecolor="white", lw=0.4)
    for yval, col in [(-1.0, "#f5c518"), (-1.5, "#e87020"), (-2.0, "#b22222")]:
        ax.axhline(yval, color=col, lw=1.0, ls="--")
    ax.axhline(0, color="black", lw=0.8)
    ax.set_ylabel(label, fontsize=11)
    ax.set_ylim(-3.5, 3.5)
    ax.set_title(title, fontsize=10)

axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%b"))
axes[1].xaxis.set_major_locator(mdates.MonthLocator())
axes[1].set_xlabel("Bulan 2024", fontsize=11)
fig.suptitle(
    "SPI-3 vs SPEI-3 — Jakarta (lat ≈ −6,2°, lon ≈ 106,8°)\n"
    "Garis putus: −1,0 Sedang / −1,5 Parah / −2,0 Ekstrem  |  DEMO: baseline 12 bln, bukan 30 thn",
    fontsize=9, y=1.01,
)
plt.tight_layout()
plt.savefig("snippet-6.png", dpi=120, bbox_inches="tight")
print("saved")

SPEI-3 cenderung menunjukkan nilai lebih negatif dari SPI-3 di musim kering. Suhu Jakarta yang tetap tinggi sepanjang tahun mempertahankan PET yang juga tinggi, sehingga defisit air terasa lebih dalam dari yang ditangkap SPI saja. Ini inti argumen Vicente-Serrano: "SPEI captures the main impact of increased temperatures on water demand."

Catatan penting: fitting pada 10–12 titik saja tidak cukup stabil untuk operasional. Dalam produksi, gunakan baseline ≥30 tahun per bulan kalender sesuai panduan WMO.

SPI versus SPEI — Kapan Menggunakan Setiap Indeks

Aspek	SPI	SPEI
Input	Presipitasi	\(P - \text{PET}\)
Distribusi	Gamma	Log-logistic
Sensitivitas suhu	Tidak ada	Ada
Cocok untuk	Data sparse, monitoring rutin	Studi perubahan iklim
Digunakan BMKG	Ya (CEWS operasional)	Terbatas

Gunakan SPI saat data suhu tidak tersedia atau saat membutuhkan baseline historis panjang di stasiun sederhana. Gunakan SPEI saat menganalisis dampak perubahan iklim atau membandingkan skenario future projection GCM — karena SPEI secara eksplisit menangkap efek "warming intensifies droughts."

Langkah Selanjutnya

Dua langkah paling langsung setelah tutorial ini:

• Perluas baseline: unduh ERA5 untuk 1991–2020 (periode referensi klimatologis WMO). Dengan 30 tahun, fitting gamma dan log-logistic jauh lebih stabil dan bisa dibandingkan langsung dengan produk BMKG.
• Gunakan library siap pakai: climate-indices di PyPI mengimplementasikan SPI dan SPEI dengan scipy backend yang sudah divalidasi dan menangani edge case lebih robust dari implementasi manual.

Eksplorasi artikel meteorologi lainnya di meteo.my.id (https://meteo.my.id) — dari analisis ENSO, pemrosesan ERA5, hingga visualisasi cuaca dengan cartopy.

Referensi

• Standardized Precipitation Index User Guide (WMO-No. 1090) — WMO, 2012.
• Standardized Precipitation Index (SPI) — UCAR/NCAR Climate Data Guide, 2022.
• Standardized Precipitation Evapotranspiration Index (SPEI) — UCAR/NCAR Climate Data Guide, 2022.
• Indeks Presipitasi Terstandarisasi — BMKG, 2026.
• Indonesian Fires Return in 2023 — NASA Earth Observatory, 2023.