Relative humidity mudah dipahami, tapi punya satu kelemahan besar: nilainya berubah hanya karena suhu udara naik — bahkan tanpa ada penambahan atau pengurangan uap air sama sekali. Itu membuat RH kurang ideal sebagai tracer kelembaban ketika udara bergerak secara adiabatik. Specific humidity (\(q\), satuan kg/kg) tidak punya masalah itu. Nilai \(q\) konstan selama parsel udara naik atau turun secara adiabatik tanpa perubahan fase, sehingga \(q\) menjadi variabel pilihan untuk analisis moisture budget, trajectory, dan profil kelembaban vertikal.
ERA5 menyediakan \(q\) langsung di pressure level (500, 850 hPa, dst.), tapi di permukaan ceritanya berbeda. ERA5 tidak mengarsipkan specific humidity permukaan secara langsung — yang tersedia adalah 2m dewpoint temperature (d2m) dan surface pressure (sp). Dari dua variabel itu, kita bisa menurunkan \(q\) permukaan menggunakan persamaan termodinamika standar.
Tutorial ini berjalan dari dua arah: kita buka file \(q\) dari pressure level, lalu kita derivasi \(q\) permukaan dari d2m dan MSLP, kemudian bandingkan keduanya di titik Jakarta dan buat plot temporal.
Memuat Data ERA5
Kita mulai dengan mengunduh dua file ERA5 dari CDS. File pertama berisi \(q\) pada pressure level 500 dan 850 hPa — ini sudah dalam satuan kg/kg, tidak perlu derivasi. File kedua berisi d2m dan MSLP untuk derivasi \(q\) permukaan di bagian berikutnya.
Untuk mengakses CDS, kita perlu akun di cds.climate.copernicus.eu dan file konfigurasi ~/.cdsapirc yang berisi API key. Setelah itu, snippet di bawah akan mengunduh file sekali; jalankan ulang dan cache guard (if not os.path.exists(OUT)) akan melewati download yang sudah ada.
import os
import cdsapi
import xarray as xr
# --- pressure-level specific humidity (500 & 850 hPa, daily 00Z, 2024) ---
OUT_Q = "era5_q_pl500-850_indonesia_2024_d.nc"
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], # N, W, S, E — Indonesia
"format": "netcdf",
},
OUT_Q,
)
ds_q = xr.open_dataset(OUT_Q)
print("Pressure-level q dataset:")
print(ds_q)
print("\nDimensions:", dict(ds_q.dims))
print("Variables:", list(ds_q.data_vars))
print("q units:", ds_q["q"].attrs.get("units", "tidak tercantum"))
Pressure-level q 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...
Dimensions: {'valid_time': 366, 'pressure_level': 2, 'latitude': 69, 'longitude': 185}
Variables: ['q']
q units: kg kg**-1
ERA5 menyimpan \(q\) pada pressure level dalam satuan kg/kg, dengan 37 level tersedia dari 1 hPa hingga 1000 hPa. Kita hanya meminta 500 dan 850 hPa karena itu rentang paling relevan untuk analisis lapisan troposfer bawah-tengah di wilayah tropis.
Sekarang kita unduh file surface — 2m dewpoint temperature dan mean sea-level pressure (MSLP). Keduanya akan digunakan di bagian berikutnya untuk menurunkan \(q\) permukaan.
# --- surface variables: d2m + MSLP (6-hourly, 2024) ---
OUT_D2M = "era5_d2m_indonesia_2024_6h.nc"
if not os.path.exists(OUT_D2M):
c = cdsapi.Client(quiet=True)
c.retrieve(
"reanalysis-era5-single-levels",
{
"product_type": "reanalysis",
"variable": ["2m_dewpoint_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_D2M,
)
OUT_MSL = "era5_msl_indonesia_2024_6h.nc"
if not os.path.exists(OUT_MSL):
c = cdsapi.Client(quiet=True)
c.retrieve(
"reanalysis-era5-single-levels",
{
"product_type": "reanalysis",
"variable": ["mean_sea_level_pressure"],
"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_MSL,
)
ds_d2m = xr.open_dataset(OUT_D2M)
ds_msl = xr.open_dataset(OUT_MSL)
print("d2m dataset:")
print(ds_d2m)
print("\nmsl dataset:")
print(ds_msl)
print("\nd2m units:", ds_d2m["d2m"].attrs.get("units", "tidak tercantum"))
print("msl units:", ds_msl["msl"].attrs.get("units", "tidak tercantum"))
d2m dataset:
<xarray.Dataset> Size: 75MB
Dimensions: (valid_time: 1464, latitude: 69, longitude: 185)
Coordinates:
* valid_time (valid_time) datetime64[ns] 12kB 2024-01-01 ... 2024-12-31T18...
expver (valid_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:
d2m (valid_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:06 GRIB to CDM+CF via cfgrib-0.9.1...
msl dataset:
<xarray.Dataset> Size: 75MB
Dimensions: (valid_time: 1464, latitude: 69, longitude: 185)
Coordinates:
* valid_time (valid_time) datetime64[ns] 12kB 2024-01-01 ... 2024-12-31T18...
expver (valid_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:
msl (valid_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:54 GRIB to CDM+CF via cfgrib-0.9.1...
d2m units: K
msl units: Pa
ERA5 menyimpan d2m dalam Kelvin dan MSLP dalam Pascal. Kita akan mengonversi keduanya ke satuan yang sesuai sebelum memasukkannya ke formula. Catatan penting: kita menggunakan MSLP sebagai proxy tekanan permukaan. Untuk wilayah dataran rendah seperti sebagian besar Jakarta dan pesisir Indonesia, pendekatan ini cukup baik karena perbedaan antara MSLP dan surface pressure di wilayah datar sangat kecil. Untuk analisis di wilayah pegunungan, sebaiknya gunakan surface pressure (sp) yang tersedia terpisah di ERA5.
Menurunkan Kelembaban Spesifik di Permukaan
Karena ERA5 tidak mengarsipkan specific humidity permukaan secara langsung, kita harus menurunkannya dari d2m dan tekanan menggunakan dua persamaan.
Langkah pertama: hitung saturation vapor pressure \(e_s\) dari dewpoint temperature \(T_d\) menggunakan formula Bolton (1980):
$$e_s(T_d) = 6{,}112 \cdot \exp\!\left(\frac{17{,}67 \, T_{d,C}}{T_{d,C} + 243{,}5}\right)$$
di mana \(T_{d,C}\) adalah suhu dewpoint dalam °C dan \(e_s\) dalam hPa. Pada suhu dewpoint sekitar 25°C (khas wilayah pesisir tropis Indonesia), formula ini menghasilkan \(e_s \approx 31{,}7\) hPa.
Langkah kedua: gunakan nilai \(e_s\) ini sebagai estimasi tekanan parsial uap air aktual \(e \approx e_s(T_d)\) — karena di titik dewpoint, udara tepat jenuh sehingga tekanan uap aktual sama dengan \(e_s\). Kemudian masukkan ke formula specific humidity:
$$q = \frac{\varepsilon \, e}{p - (1 - \varepsilon)\, e}, \quad \varepsilon \approx 0{,}622$$
di mana \(p\) adalah tekanan total atmosfer (hPa) dan \(\varepsilon = R_d / R_v \approx 0{,}622\) adalah rasio konstanta gas udara kering terhadap uap air. Nilai eksak yang digunakan dalam IFS ERA5 adalah \(\varepsilon = 0{,}621981\). Formula ini menghasilkan \(q\) dalam kg/kg.
Persamaan ini merupakan formula standar AMS Glossary yang juga digunakan secara internal di model IFS ECMWF. Perhatikan bahwa formula ini sedikit berbeda dari pendekatan sederhana \(q \approx \varepsilon e / p\) — perbedaannya kecil (orde 1–2%) tapi lebih akurat.
import numpy as np
# --- konstanta ---
EPSILON = 0.622 # R_d / R_v
def bolton_es(t_k):
"""Saturation vapor pressure (hPa) menggunakan formula Bolton 1980.
Input: suhu dalam Kelvin."""
t_c = t_k - 273.15
return 6.112 * np.exp(17.67 * t_c / (t_c + 243.5))
def specific_humidity(es_hpa, p_hpa):
"""Specific humidity (kg/kg) dari vapor pressure dan tekanan total, keduanya dalam hPa."""
return EPSILON * es_hpa / (p_hpa - (1 - EPSILON) * es_hpa)
# --- terapkan ke seluruh grid ---
# Beberapa file ERA5 dari CDS menggunakan koordinat "valid_time" alih-alih "time".
# Standarkan namanya sehingga kode di bawah tetap konsisten.
if "valid_time" in ds_d2m.coords:
ds_d2m = ds_d2m.rename({"valid_time": "time"})
if "valid_time" in ds_msl.coords:
ds_msl = ds_msl.rename({"valid_time": "time"})
d2m_k = ds_d2m["d2m"] # Kelvin
msl_pa = ds_msl["msl"] # Pascal
msl_hpa = msl_pa / 100.0 # konversi ke hPa
es_2m = bolton_es(d2m_k) # hPa, shape: (time, lat, lon)
q_2m = specific_humidity(es_2m, msl_hpa) # kg/kg
# Beri atribut yang rapi
q_2m = q_2m.assign_attrs(units="kg/kg", long_name="2m specific humidity (derived)")
q_2m.name = "q_2m"
# --- cek satu timestep: nilai spasial rerata ---
q_sample = float(q_2m.isel(time=0).mean())
es_sample = float(es_2m.isel(time=0).mean())
msl_sample = float(msl_hpa.isel(time=0).mean())
d2m_c_sample = float((d2m_k.isel(time=0) - 273.15).mean())
print(f"Timestep 0 — rata-rata spasial Indonesia:")
print(f" d2m = {d2m_c_sample:.2f} °C")
print(f" e_s(d2m) = {es_sample:.2f} hPa")
print(f" MSLP = {msl_sample:.2f} hPa")
print(f" q_2m = {q_sample*1000:.3f} g/kg")
print(f"\nBentuk array q_2m: {q_2m.shape}")
print(f"Rentang nilai: min={float(q_2m.min())*1000:.2f} g/kg, max={float(q_2m.max())*1000:.2f} g/kg")
Timestep 0 — rata-rata spasial Indonesia:
d2m = 24.38 °C
e_s(d2m) = 30.57 hPa
MSLP = 1011.21 hPa
q_2m = 19.024 g/kg
Bentuk array q_2m: (1464, 69, 185)
Rentang nilai: min=5.19 g/kg, max=26.11 g/kg
Hasilnya menunjukkan \(q\) permukaan rata-rata spasial Indonesia. Nilai yang wajar untuk wilayah tropis basah berada di kisaran 16–20 g/kg di permukaan — lebih tinggi di atas laut dan lebih bervariasi di daratan.
Membandingkan Permukaan dan 850 hPa
Setelah punya \(q\) permukaan dari derivasi, mari kita bandingkan dengan \(q\) di 850 hPa dari dataset pressure level. Kita pilih titik Jakarta (\(-6{,}2°\) LS, \(106{,}85°\) BT) dan hitung rata-rata bulanan untuk melihat pola musiman.
import pandas as pd
# Titik Jakarta
LAT_JKT = -6.2
LON_JKT = 106.85
# Standarkan nama koordinat waktu jika file ERA5 menggunakan "valid_time".
if "valid_time" in ds_q.coords:
ds_q = ds_q.rename({"valid_time": "time"})
# --- q_850 dari pressure-level dataset ---
q_850 = ds_q["q"].sel(pressure_level=850, method="nearest")
# Pilih titik grid terdekat Jakarta
q_850_jkt = q_850.sel(latitude=LAT_JKT, longitude=LON_JKT, method="nearest")
q_2m_jkt = q_2m.sel(latitude=LAT_JKT, longitude=LON_JKT, method="nearest")
# Rerata bulanan
q_850_monthly = q_850_jkt.resample(time="1ME").mean()
q_2m_monthly = q_2m_jkt.resample(time="1ME").mean()
# Tampilkan tabel per bulan
months = ["Jan","Feb","Mar","Apr","Mei","Jun","Jul","Agu","Sep","Okt","Nov","Des"]
print(f"{'Bulan':<6} {'q_2m (g/kg)':>12} {'q_850 (g/kg)':>14} {'Selisih (g/kg)':>15}")
print("-" * 50)
for i in range(len(q_850_monthly.time)):
q2 = float(q_2m_monthly.isel(time=i)) * 1000
q8 = float(q_850_monthly.isel(time=i)) * 1000
print(f"{months[i]:<6} {q2:>12.2f} {q8:>14.2f} {(q2 - q8):>15.2f}")
print()
print(f"Rata-rata tahunan q_2m : {float(q_2m_monthly.mean())*1000:.2f} g/kg")
print(f"Rata-rata tahunan q_850 : {float(q_850_monthly.mean())*1000:.2f} g/kg")
print(f"Gradien vertikal rata-rata: {(float(q_2m_monthly.mean()) - float(q_850_monthly.mean()))*1000:.2f} g/kg")
Bulan q_2m (g/kg) q_850 (g/kg) Selisih (g/kg)
--------------------------------------------------
Jan 18.77 12.87 5.89
Feb 19.01 12.77 6.24
Mar 18.89 12.82 6.07
Apr 19.68 13.25 6.44
Mei 19.33 12.63 6.70
Jun 18.59 12.20 6.40
Jul 17.03 11.19 5.85
Agu 16.93 10.80 6.12
Sep 17.15 11.49 5.66
Okt 17.52 11.41 6.12
Nov 18.77 12.27 6.50
Des 18.47 12.75 5.73
Rata-rata tahunan q_2m : 18.35 g/kg
Rata-rata tahunan q_850 : 12.20 g/kg
Gradien vertikal rata-rata: 6.14 g/kg
Output di atas menunjukkan pola yang konsisten dengan fisika atmosfer tropis: \(q\) di permukaan selalu lebih tinggi daripada di 850 hPa karena kandungan uap air berkurang seiring ketinggian (peningkatan tekanan uap dikompensasi oleh penurunan tekanan total, tapi profil kelembaban aktual turun tajam). Nilai absolut \(q\) permukaan tertinggi memang muncul pada musim basah (Desember–Maret), tetapi selisih vertikal permukaan-minus-850 hPa justru relatif datar sepanjang tahun, berkisar 5,7–6,7 g/kg, dengan nilai terbesar pada April–Mei serta November dan terendah pada Desember dan September. Pola ini mengindikasikan bahwa profil kelembaban 850 hPa di Jakarta mengikuti fluktuasi permukaan cukup dekat, sehingga gradien tidak melebar tajam saat musim basah.
Visualisasi Kelembaban Vertikal
Perbandingan tabel berguna, tapi plot timeseries memperlihatkan struktur temporal yang lebih kaya — termasuk variabilitas intraseasonal yang mungkin terkait dengan MJO maupun siklus diurnal yang teraverage ke skala bulanan.
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots(figsize=(12, 5))
# Konversi ke g/kg untuk tampilan yang lebih intuitif
times = pd.to_datetime(q_2m_monthly.time.values)
y_2m = q_2m_monthly.values * 1000
y_850 = q_850_monthly.values * 1000
ax.plot(times, y_2m, color="#1565C0", linewidth=2.0,
marker="o", markersize=5, label="$q$ permukaan (2 m, derivasi Bolton)")
ax.plot(times, y_850, color="#E64A19", linewidth=2.0,
marker="s", markersize=5, label="$q$ 850 hPa (ERA5 langsung)")
ax.fill_between(times, y_850, y_2m, alpha=0.12, color="#1565C0",
label="Selisih vertikal")
ax.set_title("Kelembaban Spesifik Jakarta 2024 — Permukaan vs 850 hPa",
fontsize=13, fontweight="bold", pad=12)
ax.set_xlabel("Bulan", fontsize=11)
ax.set_ylabel("Specific humidity (g kg$^{-1}$)", fontsize=11)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.grid(True, linestyle="--", alpha=0.5)
ax.legend(fontsize=10)
ax.set_xlim(times[0], times[-1])
# Anotasi musim basah
ax.axvspan(times[0], times[2], alpha=0.07, color="blue", label="_DJF")
ax.axvspan(times[11], times[-1], alpha=0.07, color="blue")
ax.text(times[1], ax.get_ylim()[0] + 0.3, "DJF", color="blue", fontsize=9, alpha=0.7)
plt.tight_layout()
plt.savefig("snippet-5-1.png", dpi=150, bbox_inches="tight")
print("Plot disimpan ke snippet-5-1.png")
Plot memperlihatkan separasi vertikal yang konsisten antara \(q\) permukaan dan \(q\) 850 hPa sepanjang tahun 2024. Pola musiman terlihat jelas: nilai \(q\) lebih tinggi pada bulan-bulan musim basah (Desember–Maret) ketika adveksi uap air dari Samudra Hindia dan Samudra Pasifik meningkat bersamaan dengan konveksi yang lebih aktif. Area biru transparan di antara kedua kurva menggambarkan gradien kelembaban vertikal — semakin lebar area ini, semakin besar penurunan kelembaban dari permukaan ke 850 hPa.
Penutup
Tutorial ini menunjukkan dua jalur untuk mendapatkan specific humidity dari ERA5. Jalur pertama adalah pembacaan langsung: untuk pressure level (500, 850 hPa, dst.), ERA5 sudah menyediakan \(q\) dalam kg/kg sehingga kita cukup memuat file dan memilih dimensi yang diinginkan. Jalur kedua adalah derivasi: untuk permukaan, ERA5 menyimpan d2m dan sp (atau MSLP sebagai pendekatan), lalu kita terapkan formula Bolton untuk \(e_s\) dan formula standard AMS untuk \(q\).
Nilai \(q\) sangat berguna karena bersifat konservatif dalam proses adiabatik — artinya, jika kita mengikuti parsel udara yang naik secara adiabatik tanpa kondensasi, nilai \(q\) tidak berubah. Sifat ini menjadikan \(q\) variabel tracer yang baik untuk analisis moisture budget monsun, trajectory analisis massa udara, dan perbandingan profil vertikal kelembaban antar wilayah. Ketika analisis MJO melibatkan budget uap air, \(q\) pada berbagai pressure level menjadi kuantitas kunci.
Langkah lanjutan yang bisa dicoba: perluas analisis ke grid spasial penuh dan buat peta rata-rata \(q\) per musim, bandingkan \(q\) ERA5 dengan produk satellite moisture seperti MODIS Total Precipitable Water, atau gunakan profil \(q\) sebagai input dalam perhitungan CAPE dan CIN. Coba juga ekstraksi titik yang berbeda di wilayah pegunungan seperti Bandung atau Manado untuk melihat bagaimana topografi mempengaruhi profil kelembaban vertikal, dan bandingkan hasilnya dengan mixing ratio jenuh untuk analisis stabilitas.
Eksplorasi artikel meteorologi lainnya di meteo.my.id melalui tautan resmi di https://meteo.my.id.
Referensi
- ERA5: data documentation — ECMWF (2026). Copernicus Knowledge Base, ECMWF Confluence Wiki.
- Discrepancies between the humidity at the 2m level and on pressure levels — ECMWF (n.d.). Copernicus User Support Forum, ECMWF Confluence Wiki.
- IFS Documentation CY41R2, Part IV Physical Processes — ECMWF (2016). ECMWF e-Library.
- ERA5 hourly data on pressure levels from 1940 to present — Copernicus Climate Change Service (2026). Copernicus Climate Data Store. DOI: 10.24381/cds.bd0915c6.
- Specific humidity — AMS (n.d.). Glossary of Meteorology, American Meteorological Society.
Tidak ada komentar:
Posting Komentar