Mengetahui seberapa tidak stabil kolom atmosfer adalah langkah awal yang kritis dalam prakiraan konveksi dan guntur. Salah satu indeks yang sudah lama digunakan untuk keperluan ini adalah Showalter Stability Index (SSI). Dengan data ERA5 dari ECMWF dan library MetPy, kita bisa menghitung SSI untuk seluruh domain Indonesia secara otomatis—tidak hanya untuk satu titik sounding, tapi untuk ribuan titik grid sekaligus.
Apa itu Indeks Showalter?
SSI pertama kali diperkenalkan oleh A. K. Showalter dalam makalahnya tahun 1953. Formula dasarnya:
$$\text{SI} = T_{500} - T'_{500}$$
di mana \(T_{500}\) adalah suhu lingkungan aktual pada 500 hPa, dan \(T'_{500}\) adalah suhu parsel udara pada 500 hPa setelah diangkat secara adiabatik.
Proses kalkulasinya terdiri dari tiga langkah: (1) angkat parsel dari 850 hPa secara dry-adiabatic hingga LCL (Lifting Condensation Level), yaitu ketinggian di mana suhu parsel sama dengan dewpoint 850 hPa; (2) angkat parsel dari LCL secara moist-adiabatic ke 500 hPa; (3) kurangi suhu parsel tersebut dari suhu lingkungan di 500 hPa. Nilai SI negatif berarti parsel lebih hangat dari lingkungannya—kolom tidak stabil dan konveksi mudah berkembang.
Threshold yang digunakan dalam tutorial ini mengikuti panduan EUMeTrain:
| Rentang SI | Interpretasi |
|---|---|
| \(\text{SI} > +3\) | Stabil, guntur tidak mungkin |
| \(0 < \text{SI} \leq +3\) | Marginal, guntur mungkin dengan forcing kuat |
| \(-3{,}0 \leq \text{SI} \leq 0\) | Tidak stabil, potensi tinggi guntur lebat |
| \(-6 \leq \text{SI} < -3\) | Sangat tidak stabil, guntur parah sangat mungkin |
| \(\text{SI} < -6\) | Ekstrem tidak stabil, potensi cuaca buruk hebat |
SSI berbeda dari Lifted Index (LI) hanya pada titik awal parsel: SSI menggunakan parsel fixed dari 850 hPa, sementara LI menggunakan parsel representatif dari lapisan terendah 50 hPa. Ini membuat SSI lebih tepat dipakai saat terdapat lapisan dingin dangkal di bawah 850 hPa yang menutupi potensi instabilitas sesungguhnya—misalnya saat cold surge memasuki wilayah tropis.
Satu keterbatasan penting: SSI tidak memperhitungkan pemanasan diurnal atau sumber kelembaban di bawah 850 hPa. Jika lapisan dingin meluas di atas 850 hPa, SSI cenderung meremehkan potensi konvektif sesungguhnya.
Menyiapkan dan Membaca Data ERA5
ERA5 menyediakan suhu atmosfer (\(T\)) dan specific humidity (\(q\)) pada 37 pressure level standar, termasuk 500 dan 850 hPa, dengan resolusi spasial ~0,25° (~31 km) dan temporal hourly sejak Januari 1940. Kita butuh dua file: temperature dan specific humidity di level 500 dan 850 hPa.
Download keduanya menggunakan cdsapi—daftar akun gratis di cds.climate.copernicus.eu jika belum punya. Guard if not os.path.exists(OUT) memastikan kita hanya download sekali; jalankan ulang snippet dan cdsapi langsung skip ke pembacaan file lokal.
import os
import cdsapi
import xarray as xr
import numpy as np
OUT_T = "era5_t_pl500-850_indonesia_2024_d.nc"
OUT_Q = "era5_q_pl500-850_indonesia_2024_d.nc"
if not os.path.exists(OUT_T):
c = cdsapi.Client(quiet=True)
c.retrieve(
"reanalysis-era5-pressure-levels",
{
"product_type": "reanalysis",
"variable": "temperature",
"pressure_level": ["500", "850"],
"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",
"area": [6, 95, -11, 141],
"format": "netcdf",
},
OUT_T,
)
if not os.path.exists(OUT_Q):
c = cdsapi.Client(quiet=True)
c.retrieve(
"reanalysis-era5-pressure-levels",
{
"product_type": "reanalysis",
"variable": "specific_humidity",
"pressure_level": ["500", "850"],
"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",
"area": [6, 95, -11, 141],
"format": "netcdf",
},
OUT_Q,
)
ds_t = xr.open_dataset(OUT_T)
ds_q = xr.open_dataset(OUT_Q)
print("=== Temperature dataset ===")
print(ds_t)
print("\n=== Specific humidity dataset ===")
print(ds_q)
=== Temperature dataset ===
<xarray.Dataset> Size: 37MB
Dimensions: (valid_time: 366, pressure_level: 2, latitude: 69,
longitude: 185)
Coordinates:
* valid_time (valid_time) datetime64[ns] 3kB 2024-01-01 ... 2024-12-31
expver (valid_time) <U4 6kB ...
* pressure_level (pressure_level) float64 16B 850.0 500.0
* latitude (latitude) float64 552B 6.0 5.75 5.5 ... -10.5 -10.75 -11.0
* longitude (longitude) float64 1kB 95.0 95.25 95.5 ... 140.8 141.0
number int64 8B ...
Data variables:
t (valid_time, pressure_level, latitude, longitude) float32 37MB ...
Attributes:
GRIB_centre: ecmf
GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts
GRIB_subCentre: 0
Conventions: CF-1.7
institution: European Centre for Medium-Range Weather Forecasts
history: 2026-05-10T04:10 GRIB to CDM+CF via cfgrib-0.9.1...
=== Specific humidity dataset ===
<xarray.Dataset> Size: 37MB
Dimensions: (valid_time: 366, pressure_level: 2, latitude: 69,
longitude: 185)
Coordinates:
* valid_time (valid_time) datetime64[ns] 3kB 2024-01-01 ... 2024-12-31
expver (valid_time) <U4 6kB ...
* pressure_level (pressure_level) float64 16B 850.0 500.0
* latitude (latitude) float64 552B 6.0 5.75 5.5 ... -10.5 -10.75 -11.0
* longitude (longitude) float64 1kB 95.0 95.25 95.5 ... 140.8 141.0
number int64 8B ...
Data variables:
q (valid_time, pressure_level, latitude, longitude) float32 37MB ...
Attributes:
GRIB_centre: ecmf
GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts
GRIB_subCentre: 0
Conventions: CF-1.7
institution: European Centre for Medium-Range Weather Forecasts
history: 2026-05-10T04:16 GRIB to CDM+CF via cfgrib-0.9.1...
Output di atas memperlihatkan dimensi kedua dataset: waktu (satu snapshot per hari selama 2024), pressure level (500 dan 850 hPa), latitude, dan longitude. Koordinat inilah yang akan kita gunakan untuk mengekstrak nilai pada masing-masing level.
Mengkonversi Specific Humidity menjadi Dewpoint
ERA5 menyimpan kandungan uap air sebagai specific humidity \(q\) (kg/kg), bukan dewpoint. Sebelum kita bisa menghitung SSI, \(q\) perlu dikonversi ke dewpoint terlebih dahulu. Konversi ini melewati tiga langkah: mixing ratio \(w = q/(1-q)\), tekanan uap parsial \(e = (w/(\varepsilon + w)) \cdot p\), lalu dewpoint dari \(e\) menggunakan formula Magnus—semua langkah ini sudah ditangani oleh metpy.calc.dewpoint_from_specific_humidity.
Kita pilih snapshot pertama di dataset (1 Januari 2024, 00Z) sebagai contoh representatif. Seluruh array 2D akan digunakan di snippet-snippet berikutnya.
from metpy.calc import dewpoint_from_specific_humidity
from metpy.units import units
# Deteksi nama koordinat waktu dan pressure level secara otomatis
time_coord = "valid_time" if "valid_time" in ds_t.dims else "time"
p_coord = "pressure_level" if "pressure_level" in ds_t.dims else "level"
# Ambil snapshot pertama (1 Januari 2024 00Z)
t_step = ds_t.isel({time_coord: 0})
q_step = ds_q.isel({time_coord: 0})
# Ekstrak array numpy per level
t850 = t_step["t"].sel({p_coord: 850}).values # K, shape (lat, lon)
t500 = t_step["t"].sel({p_coord: 500}).values # K
q850 = q_step["q"].sel({p_coord: 850}).values # kg/kg
q500 = q_step["q"].sel({p_coord: 500}).values # kg/kg
lats = t_step.latitude.values
lons = t_step.longitude.values
# Tampilkan beberapa sampel dewpoint 850 hPa
print("Contoh dewpoint 850 hPa (°C):")
for ii, jj in [(0, 0), (10, 20), (20, 40), (30, 60)]:
td = dewpoint_from_specific_humidity(
850 * units.hPa,
t850[ii, jj] * units.K,
q850[ii, jj] * units("kg/kg"),
)
print(
f" lat={lats[ii]:.2f}°, lon={lons[jj]:.2f}°"
f" T850={t850[ii, jj] - 273.15:.1f}°C"
f" Td850={td.to('degC').magnitude:.1f}°C"
)
Contoh dewpoint 850 hPa (°C):
lat=6.00°, lon=95.00° T850=18.3°C Td850=14.0°C
lat=3.50°, lon=100.00° T850=17.9°C Td850=17.3°C
lat=1.00°, lon=105.00° T850=18.3°C Td850=16.4°C
lat=-1.50°, lon=110.00° T850=18.5°C Td850=17.2°C
Selisih antara T850 dan Td850 di atas mencerminkan dewpoint depression—gap yang lebih kecil berarti udara lebih lembab di 850 hPa, yang penting untuk menentukan LCL dalam perhitungan SSI.
Menghitung Indeks Showalter untuk Seluruh Grid
metpy.calc.showalter_index menerima profil 1D dengan dua level (850 dan 500 hPa, urutan dari tekanan tinggi ke rendah). Untuk menerapkannya ke seluruh grid ERA5, kita loop atas setiap titik. Karena resolusi 0,25° menghasilkan grid ~69 × 185 titik, kita subsample setiap 4 titik ([::4, ::4]) untuk menjaga runtime di bawah batas sandbox—hasilnya tetap merepresentasikan pola regional dengan baik.
import warnings
warnings.filterwarnings("ignore")
from metpy.calc import showalter_index, dewpoint_from_specific_humidity
from metpy.units import units
import numpy as np
# Subsample setiap 4 titik untuk efisiensi
lats_sub = lats[::4]
lons_sub = lons[::4]
t850_sub = t850[::4, ::4]
t500_sub = t500[::4, ::4]
q850_sub = q850[::4, ::4]
q500_sub = q500[::4, ::4]
nlat, nlon = t850_sub.shape
si = np.full((nlat, nlon), np.nan)
for i in range(nlat):
for j in range(nlon):
try:
p = np.array([850.0, 500.0]) * units.hPa
T = np.array([t850_sub[i, j] - 273.15,
t500_sub[i, j] - 273.15]) * units.degC
td850 = dewpoint_from_specific_humidity(
850 * units.hPa,
t850_sub[i, j] * units.K,
q850_sub[i, j] * units("kg/kg"),
)
td500 = dewpoint_from_specific_humidity(
500 * units.hPa,
t500_sub[i, j] * units.K,
q500_sub[i, j] * units("kg/kg"),
)
td = np.array([td850.to("degC").magnitude,
td500.to("degC").magnitude]) * units.degC
ssi = showalter_index(p, T, td)
si[i, j] = ssi[0].magnitude
except Exception:
pass
print(f"Grid subsampled : {nlat} x {nlon} = {nlat * nlon} titik")
print(f"SI min : {np.nanmin(si):.1f}°C")
print(f"SI max : {np.nanmax(si):.1f}°C")
print(f"SI mean : {np.nanmean(si):.1f}°C")
Grid subsampled : 18 x 47 = 846 titik
SI min : -2.6°C
SI max : 9.9°C
SI mean : 0.8°C
Nilai min dan max di atas menggambarkan rentang instabilitas pada tanggal tersebut. Nilai SI sangat negatif biasanya terkonsentrasi di daerah konveksi aktif—Maluku, Kalimantan tengah, atau wilayah ITCZ yang sedang aktif.
Memvisualisasikan Peta Instabilitas Atmosfer
Peta 2D SSI jauh lebih informatif dibanding angka statistik saja. Kita plot menggunakan matplotlib dan cartopy di atas proyeksi PlateCarree, dengan colormap diverging yang memperlihatkan batas threshold secara langsung. Warna merah menandai daerah yang paling tidak stabil (\(\text{SI} < -3\)), dan biru menandai daerah stabil.
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
lon2d, lat2d = np.meshgrid(lons_sub, lats_sub)
# Level dan warna per band threshold
levels = [-9, -6, -3, 0, 3, 6]
cmap = plt.cm.RdBu_r
norm = mcolors.BoundaryNorm(levels, cmap.N)
fig, ax = plt.subplots(
figsize=(13, 6),
subplot_kw={"projection": ccrs.PlateCarree()},
)
cf = ax.pcolormesh(
lon2d, lat2d, si,
cmap=cmap, norm=norm,
transform=ccrs.PlateCarree(),
)
ax.add_feature(cfeature.COASTLINE, linewidth=0.9, color="black")
ax.add_feature(cfeature.BORDERS, linewidth=0.4, linestyle=":", color="gray")
ax.add_feature(cfeature.LAND, facecolor="lightgray", alpha=0.15)
ax.set_extent([95, 141, -11, 6], crs=ccrs.PlateCarree())
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color="gray", alpha=0.5, linestyle="--")
gl.top_labels = False
gl.right_labels = False
cbar = fig.colorbar(cf, ax=ax, orientation="vertical", pad=0.02, shrink=0.85, extend="both")
cbar.set_label("Showalter Index (°C)", fontsize=11)
cbar.set_ticks(levels)
ax.set_title(
"Peta Showalter Stability Index — Indonesia\nERA5 data tahun 2024, 1 Jan 00Z (grid subsampled 1°)",
fontsize=12, pad=10,
)
# Annotasi threshold
for val, label, ypos in [
(-7.5, "Ekstrem tidak stabil", 0.15),
(-4.5, "Sangat tidak stabil", 0.30),
(-1.5, "Tidak stabil", 0.48),
( 1.5, "Marginal", 0.65),
( 4.5, "Stabil", 0.82),
]:
cbar.ax.text(
2.5, ypos, label,
transform=cbar.ax.transAxes,
fontsize=7.5, va="center", color="black",
)
plt.tight_layout()
plt.savefig("si_indonesia_map.png", dpi=120, bbox_inches="tight")
print("Peta disimpan: si_indonesia_map.png")
Peta di atas menunjukkan distribusi spasial instabilitas atmosfer. Pada colormap ini, daerah merah tua (\(\text{SI} < -6\)) akan menandai zona konveksi potensial paling kuat bila ada, sementara warna biru menunjukkan kolom atmosfer yang stabil. Perhatikan bahwa untuk snapshot 1 Januari 2024 00Z ini SI minimum hanya mencapai \(-2{,}6\)°C, sehingga peta didominasi warna biru-hijau (marginal hingga tidak stabil) tanpa zona merah—kondisi yang lazim di awal tahun. Coba tanggal pada Februari–Maret atau Oktober–November untuk melihat kontras instabilitas yang lebih tajam.
Interpretasi Hasil dan Penerapan Praktis
Setelah kita punya array si dua dimensi, langkah praktis berikutnya adalah menghitung berapa banyak grid point yang jatuh ke tiap band threshold. Ini memberi gambaran cepat tentang proporsi wilayah Indonesia yang berada dalam kondisi stabil vs tidak stabil pada tanggal tersebut.
import numpy as np
valid = si[~np.isnan(si)]
total = len(valid)
bands = [
(">+3", valid > 3, "Stabil"),
("0 s.d. +3", (valid > 0) & (valid <= 3), "Marginal"),
("-3 s.d. 0", (valid >= -3) & (valid <= 0), "Tidak stabil"),
("<-3", valid < -3, "Sangat/Ekstrem tdk stabil"),
]
print(f"Total grid point valid: {total}")
print("-" * 55)
print(f"{'Band SI':<14} {'Label':<30} {'Titik':>6} {'%':>6}")
print("-" * 55)
for rng, mask, label in bands:
n = int(np.sum(mask))
print(f"{rng:<14} {label:<30} {n:>6} {100*n/total:>5.1f}%")
print("-" * 55)
Total grid point valid: 846
-------------------------------------------------------
Band SI Label Titik %
-------------------------------------------------------
>+3 Stabil 71 8.4%
0 s.d. +3 Marginal 476 56.3%
-3 s.d. 0 Tidak stabil 299 35.3%
<-3 Sangat/Ekstrem tdk stabil 0 0.0%
-------------------------------------------------------
Persentase ini berguna untuk membandingkan kondisi antar tanggal atau antar musim. Misalnya, kita bisa looping snippet-3 dan snippet-5 atas seluruh 366 hari di ERA5 data tahun 2024 dan membuat time series proporsi daerah tidak stabil—cara sederhana untuk melihat pola seasonal instabilitas atmosfer di Indonesia.
Dua batasan penting perlu diingat saat menginterpretasi hasil SSI. Pertama, SSI tidak memperhitungkan kelembaban dan pemanasan di bawah 850 hPa, sehingga bisa meremehkan instabilitas saat boundary layer sangat lembab atau saat lapisan dingin dangkal berada di bawah 850 hPa. Kedua, SSI adalah indeks termodinamik statis—ia tidak menangkap forcing dinamis (wind shear, divergensi, vortisitas) yang menentukan apakah konveksi benar-benar terjadi atau tidak. Gunakan SSI bersama indeks lain seperti CAPE, CIN, dan K-Index untuk gambaran yang lebih lengkap.
Kesimpulan dan Langkah Berikutnya
Kita telah menjalani workflow lengkap: download data ERA5 pressure-level via cdsapi, buka file dengan xarray, konversi specific humidity ke dewpoint menggunakan MetPy, hitung Showalter Index per grid point dengan loop efisien, lalu visualisasikan distribusi spasialnya di atas Indonesia menggunakan Cartopy. Seluruh workflow berjalan dari satu snapshot tunggal (1 Januari 2024 00Z) tapi mudah dikembangkan ke analisis multitemporal.
Langkah berikutnya yang layak dicoba: bandingkan SSI dengan Lifted Index untuk domain yang sama menggunakan metpy.calc.lifted_index, buat time series harian proporsi grid tidak stabil sepanjang 2024, atau integrasikan SSI ke dalam decision support sederhana yang mempertemukan threshold SSI dengan data curah hujan IMERG.
Eksplorasi artikel meteorologi lainnya di meteo.my.id — kunjungi meteo.my.id.
Referensi
- Skew-T Parameters and Indices — NOAA National Weather Service — Penjelasan lengkap NWS tentang Showalter Index, prosedur tiga langkah, tabel threshold, dan perbandingannya dengan Lifted Index.
- Glossary: Showalter Index — NOAA National Weather Service — Definisi resmi SSI sebagai indeks stabilitas untuk prakiraan guntur beserta keterbatasannya.
- Stability Indices — EUMeTrain / EUMETSAT — Tabel threshold SSI yang digunakan dalam tutorial ini: \(\text{SI} > +3\) stabil, \(-3\) s.d. \(0\) tidak stabil, \(< -3\) sangat tidak stabil.
- showalter_index — MetPy 1.7 (Unidata/UCAR) — Dokumentasi API resmi MetPy untuk
showalter_indexdandewpoint_from_specific_humidity, termasuk persyaratan pint.Quantity dan langkah konversi dari specific humidity ERA5. - ERA5: data documentation — ECMWF / Copernicus Climate Change Service — Dokumentasi resmi ERA5: 37 pressure level standar, variabel yang tersedia (termasuk \(T\) dan \(q\)), resolusi spasial ~0,25°, dan cakupan temporal.
Tidak ada komentar:
Posting Komentar