Sumber: NOAA Climate Prediction Center, Stratospheric Research Group ( Isentropic Analysis ) Kenapa Suhu Potensial Penting Bayangkan sebu...
Rabu, 24 Juni 2026
Sumber: NOAA Climate Prediction Center, Stratospheric Research Group (Isentropic Analysis)
Kenapa Suhu Potensial Penting
Bayangkan sebuah parsel udara di atas Jawa yang perlahan naik ke ketinggian 5 km. Suhu udara itu turun — bukan karena kehilangan panas, melainkan karena tekanan atmosfer berkurang sehingga parsel berekspansi. Kalau kita hanya melihat suhu biasa (temperature), kita tidak bisa membedakan apakah perubahan itu akibat gerakan adiabatik atau pertukaran panas sejati dengan lingkungan sekitar.
Di sinilah potential temperature (suhu potensial), biasa dilambangkan \(\theta\), menjadi alat yang sangat berguna. Jika udara naik atau turun tanpa pertukaran panas, bagaimana kita masih bisa mengetahui dari mana asalnya? Jawabannya: \(\theta\) tetap konstan dalam proses adiabatik kering. Setiap parsel udara "membawa" nilai \(\theta\)-nya seperti sidik jari. Dengan melacak bidang isentropis (permukaan \(\theta\) = konstan), kita bisa menelusuri jalur pergerakan massa udara meskipun ketinggiannya berubah drastis.
NOAA Climate Prediction Center menggunakan pendekatan ini untuk memantau pergerakan udara di stratosfer, termasuk melacak deplesi ozon saat musim dingin di kutub. Analisis isentropis juga menjadi fondasi untuk menghitung CAPE dan CIN, dua parameter ketidakstabilan atmosfer yang kritis untuk prakiraan konveksi dan cuaca ekstrem.
Dalam tutorial ini, kita download data temperature ERA5 pada dua level tekanan (500 dan 850 hPa) untuk wilayah Indonesia, hitung \(\theta\) secara vectorized menggunakan xarray dan NumPy, lalu buat penampang vertikal meridional untuk melihat struktur isentropisnya secara langsung.
Rumus Suhu Potensial dan Setup Data
Formula untuk potential temperature sudah dikenal sejak akhir abad ke-19, diturunkan dari persamaan Poisson untuk proses adiabatik:
$$\theta = T \left(\frac{P_0}{P}\right)^{\kappa}$$
Di sini: - \(T\) adalah suhu dalam Kelvin - \(P\) adalah tekanan dalam hPa - \(P_0 = 1000\ \text{hPa}\) adalah tekanan referensi standar - \(\kappa = R_d / c_p\) adalah konstanta Poisson
Untuk udara kering, \(R_d \approx 287\ \text{J/(kg·K)}\) dan \(c_p \approx 1004\ \text{J/(kg·K)}\), sehingga \(\kappa \approx 2/7 \approx 0{,}2857\). Dalam praktik, nilai \(0{,}286\) sering digunakan sebagai aproksimasi; ECMWF sendiri menggunakan \(0{,}28571\) dalam model IFS mereka. Kita pakai nilai eksak \(2/7\) yang sesuai konvensi MetPy dan ECMWF.
Sebagai sanity check yang berguna: parsel dengan \(T = 273\ \text{K}\) di level 850 hPa akan memiliki \(\theta \approx 286\ \text{K}\) — artinya parsel itu akan lebih panas sekitar 13 K jika kita bawa adiabatik ke permukaan referensi 1000 hPa. Intuisi ini berguna untuk memverifikasi hasil komputasi kita.
ERA5 menyediakan variabel temperature pada pressure levels dalam satuan Kelvin, tersedia di 37 level tekanan mulai dari 1000 hingga 1 hPa. Untuk tutorial ini, kita download data temperature di 500 dan 850 hPa untuk Indonesia sepanjang tahun 2024. Snippet berikut melakukan download dari CDS (hanya dijalankan sekali jika file belum ada) dan langsung membuka dataset-nya:
import os
import cdsapi
import xarray as xr
OUT = "era5_t_pl500-850_indonesia_2024_d.nc"
if not os.path.exists(OUT):
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,
)
ds = xr.open_dataset(OUT)
print(ds)
<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...
Output di atas menampilkan struktur dataset ERA5: dimensi valid_time (seluruh hari di 2024), pressure_level (dua nilai: 500 dan 850 hPa), latitude, dan longitude dengan bounding box Indonesia. Variabel utamanya adalah t — temperature dalam Kelvin.
Menghitung Theta Adiabatik untuk Semua Grid
Keunggulan xarray adalah kemampuannya untuk melakukan broadcasting secara otomatis antar dimensi. Kita tidak perlu loop manual: operasi \(\theta = T \times (P_0 / P)^\kappa\) bisa langsung diterapkan ke seluruh DataArray sekaligus, dan xarray secara otomatis menyelaraskan dimensi pressure_level dengan koordinat yang benar.
Nilai \(\kappa = 2/7\) kita pakai sebagai bilangan rasional eksak daripada float desimal, menghindari akumulasi kesalahan pembulatan saat dikalikan ke ribuan titik grid:
import numpy as np
P0 = 1000.0 # referensi tekanan dalam hPa
kappa = 2 / 7 # Rd/cp untuk udara kering (konvensi ECMWF: 0.28571)
# Ambil variabel temperature (dalam Kelvin) dan koordinat pressure_level
T = ds["t"] # shape: (time, pressure_level, latitude, longitude)
P = ds["pressure_level"] # shape: (pressure_level,) dalam hPa
# xarray broadcast P ke semua dimensi T secara otomatis
theta = T * (P0 / P) ** kappa
theta.name = "theta"
theta.attrs["units"] = "K"
theta.attrs["long_name"] = "Potential Temperature"
# Verifikasi: ambil sampel di waktu pertama, satu titik grid
print("Suhu potensial pada waktu pertama (K):")
for lev in theta.pressure_level.values:
val = float(theta.isel(valid_time=0, latitude=10, longitude=10).sel(pressure_level=lev))
T_val = float(T.isel(valid_time=0, latitude=10, longitude=10).sel(pressure_level=lev))
print(f" {int(lev):4d} hPa T = {T_val:.2f} K → θ = {val:.2f} K")
Suhu potensial pada waktu pertama (K):
850 hPa T = 291.16 K → θ = 305.00 K
500 hPa T = 268.89 K → θ = 327.78 K
Output di atas memperlihatkan nilai \(\theta\) pada dua level tekanan untuk satu titik grid. Perhatikan bahwa \(\theta\) di 850 hPa lebih tinggi dari \(T\)-nya (karena parsel itu "harus dipanaskan" bila dibawa ke 1000 hPa), sementara di 500 hPa perbedaannya lebih besar lagi karena tekanannya jauh lebih rendah. Inilah property konservasi yang kita cari: dua parsel dengan \(\theta\) yang sama berasal dari "akar adiabatik" yang sama, meski sekarang berada di ketinggian berbeda.
Visualisasi Theta dalam Penampang Vertikal
Dengan hanya dua level tekanan (500 dan 850 hPa), kita bisa membuat penampang meridional yang menampilkan distribusi \(\theta\) sepanjang garis bujur 115°E — memotong Borneo, Jawa, hingga selatan Indonesia. Irisan ini memungkinkan kita melihat bagaimana nilai \(\theta\) bervariasi terhadap lintang dan ketinggian.
Karena hanya ada dua level, kita gunakan pcolormesh untuk visualisasi filled background dan scatter untuk menandai nilai aktual di setiap level — lebih jujur daripada menginterpolasi kontur di antara dua titik saja:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
# Rata-rata theta sepanjang waktu, ambil irisan meridional di sekitar 115°E
theta_mean = theta.mean("valid_time")
lon_idx = int(np.abs(theta_mean.longitude.values - 115.0).argmin())
cross = theta_mean.isel(longitude=lon_idx) # shape: (pressure_level, latitude)
lats = cross.latitude.values
levels = cross.pressure_level.values # [500, 850] hPa
fig, ax = plt.subplots(figsize=(10, 5))
# Warna latar untuk setiap level tekanan
colors = ["#4575b4", "#d73027"] # biru = 500 hPa (lebih tinggi θ), merah = 850 hPa
for i, (lev, color) in enumerate(zip(levels, colors)):
vals = cross.sel(pressure_level=lev).values
ax.plot(lats, vals, color=color, linewidth=2.0,
label=f"{int(lev)} hPa (θ rerata {vals.mean():.1f} K)")
ax.fill_between(lats, vals, alpha=0.15, color=color)
ax.set_xlabel("Lintang (°)", fontsize=12)
ax.set_ylabel("θ (K)", fontsize=12)
ax.set_title("Potential Temperature θ — Penampang Meridional 115°E\nRerata Harian 2024 (ERA5, Indonesia)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, linestyle="--", alpha=0.5)
ax.set_xlim(lats.min(), lats.max())
plt.tight_layout()
plt.savefig("snippet-3-1.webp", dpi=120)
print("snippet-3-1.webp")
Plot di atas menampilkan profil \(\theta\) rata-rata sepanjang 2024 pada dua level tekanan, sepanjang garis bujur 115°E. Beberapa pola menarik yang kita harapkan terlihat: nilai \(\theta\) di 500 hPa secara konsisten lebih tinggi daripada di 850 hPa (artinya lapisan atas lebih stabil secara potensial terhadap lapisan bawah). Gradien \(\theta\) terhadap lintang yang relatif kecil di wilayah tropis mencerminkan karakteristik atmosfer tropis yang hampir barotropik — berbeda dengan mid-latitudes di mana gradien \(\theta\) jauh lebih tajam.
Analisis Isentropic dan Interpretasi
Setelah kita punya field \(\theta\) lengkap, pertanyaan yang lebih menarik muncul: pada level tekanan berapa isentrope \(\theta = 300\ \text{K}\) berada di atas Indonesia? Isentrope ini penting karena mewakili batas bawah yang khas antara lapisan campuran (mixed layer) dan lapisan stabil di atas Jawa dan Borneo.
Perlu diingat prinsip dasarnya: dalam proses adiabatik kering, entropy berbanding lurus dengan \(\theta\), sehingga "entropy is directly proportional to potential temperature" — parsel udara yang bergerak di sepanjang permukaan \(\theta = \text{konstan}\) tidak mengalami pertukaran panas dengan lingkungannya. Prinsip inilah yang membuat isentropis menjadi "jalur alami" pergerakan massa udara dalam atmosfer bebas.
Karena dataset kita hanya memiliki dua level (500 dan 850 hPa), kita identifikasi di level mana \(\theta\) paling dekat ke 300 K, lalu ambil statistik distribusinya di seluruh domain:
import numpy as np
TARGET_THETA = 300.0 # K — isentrope khas untuk atmosfer tropis bawah
# Hitung jarak absolut ke target di setiap grid dan waktu
diff = np.abs(theta - TARGET_THETA) # shape: (time, pressure_level, lat, lon)
# Temukan index level yang paling dekat ke 300 K, per titik grid dan waktu
nearest_idx = diff.argmin(dim="pressure_level") # shape: (time, lat, lon), nilai 0 atau 1
# Presisi level: 0 = 500 hPa, 1 = 850 hPa (tergantung urutan pressure_level)
plevels = theta.pressure_level.values # [500., 850.] atau [850., 500.]
print(f"Urutan level dalam dataset: {plevels} hPa")
# Rata-rata nearest_idx terhadap waktu
avg_idx = nearest_idx.mean("valid_time") # nilai 0.0 berarti selalu 500 hPa, 1.0 = selalu 850 hPa
frac_850 = float((nearest_idx == 1).mean("valid_time").mean())
frac_500 = float((nearest_idx == 0).mean("valid_time").mean())
print(f"\nIsentrope θ = {TARGET_THETA} K ditemukan lebih dekat ke:")
print(f" 500 hPa: {frac_500*100:.1f}% dari semua titik grid × waktu")
print(f" 850 hPa: {frac_850*100:.1f}% dari semua titik grid × waktu")
# Hitung θ aktual di 850 hPa dan 500 hPa untuk konteks
theta_850 = float(theta.sel(pressure_level=850).mean())
theta_500 = float(theta.sel(pressure_level=500).mean())
print(f"\nRerata θ domain sepanjang 2024:")
print(f" 850 hPa: {theta_850:.1f} K")
print(f" 500 hPa: {theta_500:.1f} K")
print(f"\nKesimpulan: isentrope θ = {TARGET_THETA} K berada {'di antara' if theta_850 < TARGET_THETA < theta_500 else 'di luar rentang'} 850–500 hPa.")
Urutan level dalam dataset: [850. 500.] hPa
Isentrope θ = 300.0 K ditemukan lebih dekat ke:
500 hPa: 100.0% dari semua titik grid × waktu
850 hPa: 0.0% dari semua titik grid × waktu
Rerata θ domain sepanjang 2024:
850 hPa: 305.4 K
500 hPa: 327.8 K
Kesimpulan: isentrope θ = 300.0 K berada di luar rentang 850–500 hPa.
Output snippet-4 memberikan gambaran kuantitatif posisi isentrope \(\theta = 300\ \text{K}\) di wilayah Indonesia. Rerata domain menunjukkan \(\theta\) di 850 hPa sudah mencapai 305{,}4 K, artinya isentrope 300 K berada di bawah level 850 hPa — di lapisan troposfer paling bawah, sekitar 900–950 hPa di atas Indonesia. Ini masuk akal: atmosfer tropis Indonesia hangat dan lembap, sehingga udara dekat permukaan sudah memiliki \(\theta\) yang cukup tinggi meskipun secara teknik berada di lapisan batas (boundary layer).
Parsel udara yang bergerak sepanjang isentrope ini — misalnya udara maritim yang mengalir dari Samudra Hindia menuju pantai barat Sumatra — dalam kondisi ideal tidak bertukar panas dengan sekitarnya. Analisis isentropis memungkinkan kita "mengikuti" massa udara tersebut bahkan ketika jalurnya melintasi batas-batas level tekanan yang kita gunakan.
Langkah Selanjutnya
Tutorial ini memperkenalkan komputasi \(\theta\) dari ERA5 sebagai langkah pertama dalam analisis termodinamika atmosfer. Ada beberapa arah yang menarik untuk dieksplorasi lebih lanjut:
Equivalent potential temperature (\(\theta_e\)) memperluas konsep yang sama ke udara lembap. Alih-alih mengabaikan kandungan air, \(\theta_e\) memperhitungkan panas laten kondensasi — sangat relevan untuk atmosfer tropis Indonesia yang kaya uap air. Perbedaan antara \(\theta\) dan \(\theta_e\) langsung mencerminkan kandungan moisture massa udara.
Diagnosa stabilitas menjadi jauh lebih tepat saat menggunakan \(\theta\) daripada suhu biasa. Kalau \(\theta\) meningkat terhadap ketinggian (\(\partial\theta/\partial z > 0\)), lapisan itu stabil secara adiabatik kering. Profil \(\theta\) dari ERA5 bisa langsung dibandingkan dengan profil radiosonde untuk validasi model.
Output GFS real-time menggunakan format dan konvensi yang sangat mirip ERA5 pada pressure levels. Pola yang sama — buka NetCDF, ambil variabel temperature, hitung \(\theta\) dengan formula Poisson — berlaku langsung untuk data GFS yang kita download lewat NOMADS atau Open-Meteo. Ini memungkinkan perbandingan antara analisis ERA5 dengan prakiraan numerik cuaca terbaru.
Potential temperature adalah fondasi yang solid. Dari sini, CAPE dan CIN, analisis frontal, dan diagnosa konveksi semuanya menjadi lebih mudah dipahami secara fisis.
Eksplorasi artikel meteorologi lainnya di meteo.my.id — buka di sini.
Referensi
- NOAA National Weather Service Glossary — Potential Temperature — definisi resmi NWS: temperature a parcel of dry air would have if brought adiabatically to 1000 mb.
- NOAA Climate Prediction Center — Isentropic Analysis Information — penjelasan CPC tentang prinsip konservasi entropy dan penggunaan isentropis untuk memantau ozon stratosfer.
- EUMeTrain / EUMETSAT — Potential Temperature — derivasi formula Poisson (\(\theta = T \times (1000/p)^{0{,}286}\)) dan aplikasi analisis mesoskal dan sinoptik.
- Unidata MetPy — metpy.calc.potential_temperature — implementasi Python menggunakan \(\kappa = R_d/c_p \approx 0{,}2857\), kompatibel dengan xarray ERA5.
- ECMWF — ERA5 Data Documentation — dokumentasi lengkap ERA5: 37 pressure levels, resolusi 0,28125°, coverage 1940–sekarang.
Selasa, 23 Juni 2026
Sumber: GOES Project Science Office / NOAA (The Intertropical Convergence Zone)
Indonesia terletak persis di jalur ITCZ — Intertropical Convergence Zone, tempat angin pasat dari belahan bumi utara dan selatan saling bertemu. Pertemuan itu tidak berhenti di permukaan: udara yang berkonvergensi dipaksa naik, membentuk awan cumulonimbus raksasa, dan akhirnya menghasilkan hujan yang mendominasi musim basah kita. Tapi seberapa kuat konvergensi itu, dan di mana tepatnya lokasinya pada tanggal tertentu?
Jawabannya ada di data angin permukaan ERA5. Dengan menghitung divergensi horizontal dari komponen u10 dan v10, kita bisa mengukur langsung seberapa besar angin "terkumpul" (konvergensi, nilai negatif) atau "menyebar" (divergensi, nilai positif) di setiap titik grid. Tutorial ini memandu kita dari download data via cdsapi, koreksi koordinat sferis, visualisasi peta divergensi dengan cartopy, hingga perbandingan musiman DJF versus JJA.
Konsep Divergensi Horizontal
Divergensi horizontal mengukur seberapa cepat angin menyebar atau mengumpul di suatu area permukaan. Secara matematis, divergensi adalah dot product dari operator nabla dengan vektor angin horizontal:
$$\nabla \cdot \mathbf{V} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}$$
Dalam notasi sederhana ini, \(u\) adalah komponen angin zonal (barat-timur) dan \(v\) adalah komponen meridional (utara-selatan). Satuannya adalah \(\text{s}^{-1}\).
- Positif (\(\nabla \cdot \mathbf{V} > 0\)): divergensi — angin menyebar keluar dari suatu titik. Di permukaan, ini terkait dengan subsidensi dan cuaca cerah.
- Negatif (\(\nabla \cdot \mathbf{V} < 0\)): konvergensi — angin mengalir masuk ke suatu titik. Di permukaan, ini memicu gerakan naik dan pembentukan awan hujan.
Formula Cartesius di atas hanya valid untuk grid datar. Di grid ERA5 yang menggunakan koordinat lintang-bujur sferis, kita perlu koreksi:
$$\nabla \cdot \mathbf{V} = \frac{1}{a \cos\varphi} \frac{\partial u}{\partial \lambda} + \frac{1}{a \cos\varphi} \frac{\partial (v \cos\varphi)}{\partial \varphi}$$
Di sini \(a\) adalah jari-jari Bumi (\(6{,}371 \times 10^6\) m), \(\varphi\) adalah lintang, dan \(\lambda\) adalah bujur. Faktor \(\cos\varphi\) muncul karena jarak antar-bujur semakin menyempit menuju kutub. Tanpa koreksi ini, nilai divergensi di lintang rendah (seperti Indonesia) akan meleset karena grid ERA5 bukan proyeksi Cartesius.
Keterkaitan ITCZ sangat langsung: ketika ITCZ berada di atas Indonesia (DJF — musim hujan), divergensi permukaan bernilai sangat negatif. Ketika ITCZ bergeser ke belahan bumi utara (JJA), nilai divergensi di Indonesia cenderung mendekati nol atau bahkan positif.
Mengunduh Data ERA5 u10 dan v10
ERA5 menyediakan komponen angin permukaan 10 m — u10 (paramId 165) dan v10 (paramId 166) — pada resolusi \(0{,}25°\) setiap 6 jam. Kita download keduanya untuk seluruh tahun 2024 dengan bounding box Indonesia (\(6°\)N, \(95°\)E, \(-11°\)S, \(141°\)E). Daftar akun CDS di cds.climate.copernicus.eu lalu install cdsapi sebelum menjalankan snippet berikut.
import os
import cdsapi
import xarray as xr
OUT_U = "era5_u10_indonesia_2024_6h.nc"
OUT_V = "era5_v10_indonesia_2024_6h.nc"
if not os.path.exists(OUT_U):
c = cdsapi.Client(quiet=True)
c.retrieve("reanalysis-era5-single-levels", {
"product_type": "reanalysis",
"variable": ["10m_u_component_of_wind"],
"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_U)
if not os.path.exists(OUT_V):
c = cdsapi.Client(quiet=True)
c.retrieve("reanalysis-era5-single-levels", {
"product_type": "reanalysis",
"variable": ["10m_v_component_of_wind"],
"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_V)
ds_u = xr.open_dataset(OUT_U)
ds_v = xr.open_dataset(OUT_V)
# ERA5 dari CDS kadang memakai koordinat 'valid_time'; rename ke 'time' agar konsisten
if "valid_time" in ds_u.dims:
ds_u = ds_u.rename({"valid_time": "time"})
if "valid_time" in ds_v.dims:
ds_v = ds_v.rename({"valid_time": "time"})
print("Dataset u10:")
print(ds_u)
print("\nDataset v10:")
print(ds_v)
print(f"\nWaktu tersedia: {ds_u.time.values[0]} s.d. {ds_u.time.values[-1]}")
print(f"Resolusi grid: {abs(float(ds_u.latitude.values[1] - ds_u.latitude.values[0])):.2f}°")
Dataset u10:
<xarray.Dataset> Size: 75MB
Dimensions: (time: 1464, latitude: 69, longitude: 185)
Coordinates:
* time (time) datetime64[ns] 12kB 2024-01-01 ... 2024-12-31T18:00:00
expver (time) <U4 23kB ...
* latitude (latitude) float64 552B 6.0 5.75 5.5 5.25 ... -10.5 -10.75 -11.0
* longitude (longitude) float64 1kB 95.0 95.25 95.5 ... 140.5 140.8 141.0
number int64 8B ...
Data variables:
u10 (time, latitude, longitude) float32 75MB ...
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-10T03:59 GRIB to CDM+CF via cfgrib-0.9.1...
Dataset v10:
<xarray.Dataset> Size: 75MB
Dimensions: (time: 1464, latitude: 69, longitude: 185)
Coordinates:
* time (time) datetime64[ns] 12kB 2024-01-01 ... 2024-12-31T18:00:00
expver (time) <U4 23kB ...
* latitude (latitude) float64 552B 6.0 5.75 5.5 5.25 ... -10.5 -10.75 -11.0
* longitude (longitude) float64 1kB 95.0 95.25 95.5 ... 140.5 140.8 141.0
number int64 8B ...
Data variables:
v10 (time, latitude, longitude) float32 75MB ...
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:02 GRIB to CDM+CF via cfgrib-0.9.1...
Waktu tersedia: 2024-01-01T00:00:00.000000000 s.d. 2024-12-31T18:00:00.000000000
Resolusi grid: 0.25°
Output menunjukkan dimensi (time, latitude, longitude) dengan 1464 time step untuk u10 dan v10, masing-masing mencakup seluruh 2024. Kita juga bisa konfirmasi resolusi \(0{,}25°\) dari perbedaan antar-lintang.
Menghitung Divergensi Permukaan
Dengan data di tangan, kita pilih satu time slice representatif dari musim DJF — 15 Januari 2024 pukul 00Z — untuk menghitung divergensi menggunakan formula sferis. Implementasinya memakai numpy.gradient dengan koreksi faktor \(\cos\varphi\) secara eksplisit.
import numpy as np
# Select time: Jan 15, 2024 00Z (DJF — ITCZ aktif di atas Indonesia)
TIME = "2024-01-15T00:00:00"
u = ds_u["u10"].sel(time=TIME, method="nearest")
v = ds_v["v10"].sel(time=TIME, method="nearest")
lat = u.latitude.values # degrees
lon = u.longitude.values # degrees
R = 6371000.0 # jari-jari Bumi dalam meter
lat_rad = np.deg2rad(lat)
dlon = np.deg2rad(np.gradient(lon))
dlat = np.deg2rad(np.gradient(lat))
# Koreksi sferis: div = 1/(R cosφ) * ∂u/∂λ + 1/(R cosφ) * ∂(v cosφ)/∂φ
cos_lat = np.cos(lat_rad)[:, np.newaxis] # shape (nlat, 1)
v_cos = v.values * np.cos(lat_rad)[:, np.newaxis]
du_dlambda = np.gradient(u.values, axis=1) / dlon[np.newaxis, :]
dvcoslat_dphi = np.gradient(v_cos, axis=0) / dlat[:, np.newaxis]
div = (du_dlambda + dvcoslat_dphi) / (R * cos_lat)
div_scaled = div * 1e5 # konversi ke 10⁻⁵ s⁻¹
print(f"Waktu: {TIME}")
print(f"Domain mean divergensi : {div_scaled.mean():.4f} × 10⁻⁵ s⁻¹")
print(f"Minimum divergensi (konvergensi terkuat): {div_scaled.min():.4f} × 10⁻⁵ s⁻¹")
print(f"Maximum divergensi : {div_scaled.max():.4f} × 10⁻⁵ s⁻¹")
# Lokasi konvergensi terkuat
lat_idx, lon_idx = np.unravel_index(div_scaled.argmin(), div_scaled.shape)
print(f"Konvergensi terkuat di: lat={lat[lat_idx]:.2f}°, lon={lon[lon_idx]:.2f}°")
Waktu: 2024-01-15T00:00:00
Domain mean divergensi : -0.2492 × 10⁻⁵ s⁻¹
Minimum divergensi (konvergensi terkuat): -11.6014 × 10⁻⁵ s⁻¹
Maximum divergensi : 12.2946 × 10⁻⁵ s⁻¹
Konvergensi terkuat di: lat=-5.75°, lon=114.25°
Output berikut menampilkan nilai domain-mean divergensi dan lokasi konvergensi terkuat pada tanggal tersebut. Nilai minimum yang sangat negatif (biasanya di kisaran \(-5\) hingga \(-10 \times 10^{-5}\ \text{s}^{-1}\) pada DJF) mengindikasikan zona ITCZ yang aktif, sementara nilai positif di area lain menandai region dengan angin yang menyebar.
Visualisasi Peta Divergensi dan Angin
Peta divergensi lebih mudah dibaca ketika kita overlay dengan vektor angin, sehingga terlihat langsung ke mana angin bergerak dan di mana terjadi penumpukan massa udara. Kita pakai pcolormesh dengan colormap RdBu_r — merah untuk divergensi, biru untuk konvergensi — dan quiver untuk tanda panah angin.
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig, ax = plt.subplots(figsize=(12, 7), subplot_kw={"projection": ccrs.PlateCarree()})
ax.set_extent([95, 141, -11, 6], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.4, linestyle=":")
gl = ax.gridlines(draw_labels=True, linewidth=0.3, color="gray", alpha=0.5)
gl.top_labels = False
gl.right_labels = False
pcm = ax.pcolormesh(
lon, lat, div_scaled,
cmap="RdBu_r", vmin=-8, vmax=8,
transform=ccrs.PlateCarree()
)
cb = plt.colorbar(pcm, ax=ax, label="Divergensi (×10⁻⁵ s⁻¹)", shrink=0.75, pad=0.05)
cb.ax.tick_params(labelsize=9)
# Subsample vektor angin agar tidak terlalu padat
step = 4
ax.quiver(
lon[::step], lat[::step],
u.values[::step, ::step], v.values[::step, ::step],
transform=ccrs.PlateCarree(),
scale=50, width=0.003, color="black", alpha=0.7,
label="Angin 10 m"
)
ax.set_title(
"Divergensi Angin Permukaan ERA5 — 15 Januari 2024 00Z\n"
"(Merah = Divergensi, Biru = Konvergensi)",
fontsize=12, pad=10
)
plt.tight_layout()
plt.savefig("div_map.png", dpi=120, bbox_inches="tight")
print("Peta divergensi tersimpan sebagai div_map.png")
Pada peta DJF, zona biru pekat seharusnya tampak memanjang di sekitar ekuator dan beberapa derajat selatan — itu signature ITCZ. Perhatikan bagaimana vektor angin dari arah barat laut dan barat daya "bertemu" di zona biru tersebut, mengonfirmasi konvergensi fisik yang mendorong naiknya massa udara.
Analisis Musiman Konvergensi
Satu snapshot DJF memberi gambaran sesaat. Untuk memahami pola musiman, kita hitung rata-rata divergensi selama DJF (Januari–Februari 2024) dan JJA (Juni–Agustus 2024), lalu bandingkan nilai domain-mean-nya.
Hipotesisnya: selama DJF, ITCZ berada di dekat atau selatan ekuator sehingga divergensi domain Indonesia lebih negatif (konvergensi lebih kuat). Selama JJA, ITCZ bermigrasi ke belahan bumi utara, meninggalkan Indonesia dengan angin yang lebih divergen dan musim kemarau.
u_all = ds_u["u10"]
v_all = ds_v["v10"]
def compute_div_mean(u_sel, v_sel):
lat_r = np.deg2rad(u_sel.latitude.values)
lon_ = u_sel.longitude.values
dlon_ = np.deg2rad(np.gradient(lon_))
dlat_ = np.deg2rad(np.gradient(lat_r))
cos_ = np.cos(lat_r)[:, np.newaxis]
v_cos_ = v_sel.values * np.cos(lat_r)[:, np.newaxis]
du_dl = np.gradient(u_sel.values, axis=1) / dlon_[np.newaxis, :]
dvc_dp = np.gradient(v_cos_, axis=0) / (dlat_[:, np.newaxis] * R)
div_ = (du_dl / (R * cos_)) + dvc_dp
return div_ * 1e5
# DJF: Januari dan Februari 2024 (Desember 2024 di luar dataset)
djf_mask = u_all.time.dt.month.isin([1, 2])
u_djf = u_all.sel(time=djf_mask).mean("time")
v_djf = v_all.sel(time=djf_mask).mean("time")
div_djf = compute_div_mean(u_djf, v_djf)
# JJA: Juni–Agustus 2024
jja_mask = u_all.time.dt.month.isin([6, 7, 8])
u_jja = u_all.sel(time=jja_mask).mean("time")
v_jja = v_all.sel(time=jja_mask).mean("time")
div_jja = compute_div_mean(u_jja, v_jja)
print("=== Divergensi Musiman (Domain Indonesia, 2024) ===")
print(f"DJF (Jan–Feb 2024): mean = {div_djf.mean():.4f} × 10⁻⁵ s⁻¹, min = {div_djf.min():.4f}")
print(f"JJA (Jun–Agu 2024): mean = {div_jja.mean():.4f} × 10⁻⁵ s⁻¹, min = {div_jja.min():.4f}")
print()
if div_djf.mean() < div_jja.mean():
print("DJF lebih konvergen (nilai lebih negatif) → ITCZ aktif di atas Indonesia saat musim hujan.")
else:
print("JJA lebih konvergen → analisis lebih lanjut diperlukan.")
=== Divergensi Musiman (Domain Indonesia, 2024) ===
DJF (Jan–Feb 2024): mean = -10.3682 × 10⁻⁵ s⁻¹, min = -371.9647
JJA (Jun–Agu 2024): mean = -2.0960 × 10⁻⁵ s⁻¹, min = -391.9321
DJF lebih konvergen (nilai lebih negatif) → ITCZ aktif di atas Indonesia saat musim hujan.
Output menunjukkan perbedaan yang jelas antara dua musim. Nilai mean DJF yang lebih negatif mengonfirmasi posisi ITCZ di atas atau dekat Indonesia, sementara nilai JJA yang mendekati nol (atau bahkan positif) mencerminkan migrasi ITCZ ke utara. Ini konsisten dengan pola monsun Asia-Australia yang mendominasi iklim kepulauan kita. MJO juga bisa memodulasi sinyal ini — fase aktif MJO memperkuat konvergensi permukaan dan divergensi atas, sementara fase teredam membalik polanya.
Langkah Berikutnya
Tutorial ini menunjukkan cara dasar menghitung dan memvisualisasikan divergensi permukaan dari ERA5. Ada beberapa arah eksplorasi yang menarik sebagai kelanjutan:
-
Gerak vertikal dari persamaan kontinuitas — divergensi permukaan yang kita hitung bisa diintegrasikan secara vertikal untuk memperkirakan kecepatan vertikal (\(\omega\)). Ini cara sederhana memverifikasi apakah area konvergensi memang bertepatan dengan kolom udara yang naik.
-
Korelasi dengan curah hujan IMERG — overlay data konvergensi bulanan dengan IMERG precipitation akan memperlihatkan seberapa erat korelasi spasial antara zona konvergensi dan curah hujan aktual di Indonesia.
-
Pengaruh MJO — filter data per fase MJO (Wheeler-Hendon index) lalu hitung komposit divergensi per fase. Artikel terpisah di blog ini akan membahas cara membuat diagram Wheeler-Hendon dari data ERA5.
-
Divergensi level atas (200 hPa) — di troposfer atas, pola sebaliknya berlaku: divergensi kuat di 200 hPa terkait kolom konveksi aktif. MetPy menyediakan
metpy.calc.divergence()yang otomatis menangani koreksi sferis jika input-nya adalah xarray DataArray dengan metadata proyeksi.
Eksplorasi artikel meteorologi lainnya di meteo.my.id — dari analisis deret waktu cuaca hingga pemetaan angin dengan cartopy, kunjungi meteo.my.id.
Referensi
- NOAA JetStream — Convergence Zone — penjelasan ITCZ sebagai zona pertemuan angin pasat dari kedua belahan bumi beserta mekanisme pembentukan awan hujannya.
- NASA Earth Observatory — The Intertropical Convergence Zone — gambaran visual dan saintifik ITCZ termasuk variabilitasnya musiman dan dampak pada pola curah hujan tropis.
- NOAA Climate.gov — What is the MJO and why do we care? — penjelasan MJO dan bagaimana fase aktif/teredam memodulasi konvergensi permukaan serta curah hujan di kawasan tropis.
- MetPy API — metpy.calc.divergence — dokumentasi fungsi divergensi MetPy yang menangani koreksi sferis secara otomatis menggunakan metadata xarray.
- ECMWF — ERA5 Data Documentation — dokumentasi resmi ERA5 termasuk daftar variabel, resolusi, dan cara akses via CDS API.
Sumber: GOES Project Science Office / NOAA ( The Intertropical Convergence Zone ) Indonesia terletak persis di jalur ITCZ — Intertropical ...
