TANGERANG SELATAN WEATHER

Sabtu, 04 Juli 2026

Menginterpolasi Data ERA5 ke Ketinggian Tetap

Juli 04, 2026 data science , meteorologi 0 Comments

Gradien vertikal dari biru dingin di troposfer atas ke oranye hangat di permukaan, merepresentasikan struktur pressure levels ERA5

Pengenalan

ERA5, produk reanalysis global dari ECMWF, menyimpan hampir semua variabel atmosfernya di pressure levels — 37 level standar dari 1 hPa di stratosfer atas sampai 1000 hPa mendekati permukaan. Pilihan koordinat tekanan ini punya alasan yang kuat: setiap lapisan tekanan mengandung massa atmosfer yang konstan, sehingga persamaan hidrostatik dan termodinamik dalam NWP menjadi lebih sederhana dan komputasi lebih efisien. ERA5 sendiri diproduksi secara native di 137 hybrid sigma-pressure model levels, lalu diinterpolasi ke 37 pressure levels oleh post-processor FULL-POS dalam IFS sebelum diarsipkan.

Masalah muncul ketika kita butuh nilai variabel atmosfer di ketinggian geometris tertentu — misalnya 3000 m di atas permukaan laut. Ini situasi yang sering terjadi dalam praktik:

  • Validasi terhadap rawinsonde: balon cuaca melaporkan profil T, RH, dan angin di berbagai ketinggian geometris; membandingkannya langsung dengan ERA5 butuh interpolasi ke ketinggian yang sama.
  • Analisis stabilitas: menghitung CAPE dan CIN memerlukan profil suhu dan kelembaban pada ketinggian tertentu dalam lapisan campuran (mixed layer) — bukan di pressure level tetap.
  • Perbandingan multi-model: menyandingkan ERA5 dengan output WRF atau GFS lebih bersih bila keduanya sudah diinterpolasi ke grid ketinggian yang seragam.

Tutorial ini memandu prosesnya dari awal sampai akhir: kita download data ERA5 pressure levels lewat cdsapi, konversi geopotensial ke geopotential height menggunakan konstanta gravitasi standar WMO, lalu interpolasi suhu ke 3000 m dengan linear interpolation. Di akhir kita visualisasikan profil vertikal hasilnya.

Mengunduh dan Membaca Data ERA5

ERA5 mengarsipkan geopotensial sebagai variabel z (parameter ID 129, satuan m² s⁻²) dan suhu sebagai variabel t, masing-masing tersedia di setiap pressure level. Kita butuh dua file terpisah — satu untuk geopotensial, satu untuk suhu — yang masing-masing mencakup pressure levels 500 dan 850 hPa, resolusi harian 00 UTC, wilayah Indonesia, tahun 2024. Sebelum menjalankan snippet ini, pastikan kamu sudah mendaftar dan mengkonfigurasi kredensial CDS di cds.climate.copernicus.eu.

Download berjalan sekali; run berikutnya melewati blok if dan langsung membaca dari file lokal:

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

OUT_Z = "era5_z_pl500-850_indonesia_2024_d.nc"
OUT_T = "era5_t_pl500-850_indonesia_2024_d.nc"

if not os.path.exists(OUT_Z) or not os.path.exists(OUT_T):
    c = cdsapi.Client(quiet=True)
    if not os.path.exists(OUT_Z):
        c.retrieve(
            "reanalysis-era5-pressure-levels",
            {
                "product_type": "reanalysis",
                "variable": ["geopotential"],
                "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_Z,
        )
    if not os.path.exists(OUT_T):
        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,
        )

ds_z = xr.open_dataset(OUT_Z)
ds_t = xr.open_dataset(OUT_T)
print("Geopotential:")
print(ds_z)
print("\nTemperature:")
print(ds_t)

Geopotential:
<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:
    z               (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:08 GRIB to CDM+CF via cfgrib-0.9.1...

Temperature:
<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 mengonfirmasi struktur dataset: dua pressure levels (500 dan 850 hPa), 366 timestep daily 00 UTC untuk seluruh tahun 2024 (tahun kabisat), dan grid spasial 0,25° × 0,25° yang mencakup kotak Indonesia.

Konversi Geopotensial ke Ketinggian Geopotensial

ERA5 menyimpan geopotensial \(z\) dalam satuan m² s⁻², bukan ketinggian langsung. Geopotensial merepresentasikan kerja yang diperlukan untuk mengangkat satu satuan massa dari permukaan laut ke ketinggian tersebut melawan gravitasi. Untuk mendapatkan geopotential height \(H\) dalam satuan metre, kita bagi dengan konstanta gravitasi standar WMO:

$$H = \frac{z}{g_0}$$

Nilai \(g_0 = 9{,}80665\ \text{m s}^{-2}\) adalah konstanta yang digunakan dalam IFS ECMWF — uniform untuk semua lintang dan semua ketinggian.

Geopotential height \(H\) berbeda sedikit dari ketinggian geometris (altitude) nyata. Hubungan konversinya adalah:

$$\text{alt} = \frac{R_e \cdot H}{R_e - H}$$

dengan \(R_e = 6\,371\,229\ \text{m}\) (radius bumi yang digunakan ECMWF dalam model IFS). Koreksi ini kecil di troposfer: pada 10 km, perbedaannya sekitar 0,16% atau ~16 m. Di level 850 hPa (sekitar 1,5 km), selisihnya kurang dari 1 metre. Untuk keperluan tutorial ini kita hitung keduanya supaya angkanya terlihat langsung.

Mengapa memakai geopotential height alih-alih ketinggian geometris? Karena dengan \(H\), nilai \(|g|\) bisa diperlakukan sebagai konstanta \(g_0\) dalam seluruh persamaan atmosfer, meskipun secara fisik gravitasi sedikit berkurang dengan ketinggian. Ini yang membuat persamaan termodinamik atmosfer lebih bersih — dan itulah alasan ERA5 serta hampir semua model NWP bekerja dalam koordinat geopotensial.

g0 = 9.80665        # konstanta gravitasi standar WMO (m s⁻²)
Re = 6_371_229.0    # radius bumi ECMWF IFS (m)

# Grid point terdekat Jakarta
lat0, lon0 = -6.2, 106.8

# Ambil geopotensial di titik itu, timestep pertama (2024-01-01 00 UTC)
z_pt = ds_z["z"].sel(latitude=lat0, longitude=lon0, method="nearest").isel(valid_time=0)

# Konversi ke geopotential height: H = z / g0
H = z_pt / g0  # satuan: metre

H_850 = float(H.sel(pressure_level=850))
H_500 = float(H.sel(pressure_level=500))

# Ketinggian geometris sebagai perbandingan
alt_850 = Re * H_850 / (Re - H_850)
alt_500 = Re * H_500 / (Re - H_500)

print(f"Grid point terdekat Jakarta ({lat0}°, {lon0}°) — 2024-01-01 00 UTC")
print(f"  850 hPa:  H = {H_850:.1f} m   |  alt = {alt_850:.1f} m   (selisih {alt_850-H_850:.2f} m)")
print(f"  500 hPa:  H = {H_500:.1f} m   |  alt = {alt_500:.1f} m   (selisih {alt_500-H_500:.2f} m)")

Grid point terdekat Jakarta (-6.2°, 106.8°) — 2024-01-01 00 UTC
  850 hPa:  H = 1512.1 m   |  alt = 1512.5 m   (selisih 0.36 m)
  500 hPa:  H = 5896.5 m   |  alt = 5902.0 m   (selisih 5.46 m)

Nilai yang tercetak mengonfirmasi dua hal: level 850 hPa di wilayah Jakarta pada 1 Januari 2024 berada di sekitar 1400–1600 m, dan level 500 hPa di sekitar 5700–6000 m — kisaran yang kita harapkan untuk kondisi tropis. Selisih geometric vs geopotential height di level-level ini hanya beberapa metre; untuk sebagian besar analisis meteorologi, koreksi geometris bisa diabaikan dan kita cukup memakai \(H\) langsung.

Interpolasi Vertikal ke Ketinggian Tetap

Dengan geopotential height di kedua level sudah diketahui, interpolasi suhu ke ketinggian target adalah langkah yang langsung. Metode paling sederhana adalah linear interpolation dalam ruang ketinggian: kita asumsikan suhu berubah secara linier antara dua titik yang diketahui.

Alternatif yang lebih akurat secara fisik adalah interpolasi linear-dalam-log-tekanan — menggunakan \(\ln(p)\) sebagai sumbu vertikal alih-alih \(H\). Pendekatan ini bekerja lebih baik karena suhu dan kebanyakan variabel termodinamik lain cenderung linier terhadap log-tekanan di troposfer bebas. Untuk tutorial ini dengan hanya dua level (850 dan 500 hPa), perbedaan hasil antara kedua metode kecil; kita gunakan linear-in-height untuk kesederhanaan. Di bagian penutup kita bahas cara beralih ke metode log-tekanan.

Target kita adalah 3000 m — ketinggian yang berada di antara \(H_{850}\) dan \(H_{500}\), jadi tidak ada ekstrapolasi. Kita pakai scipy.interpolate.interp1d yang bisa diperluas ke lebih banyak level tanpa mengubah struktur kode secara signifikan.

from scipy.interpolate import interp1d

# Suhu di grid point dan timestep yang sama
T_pt = ds_t["t"].sel(latitude=lat0, longitude=lon0, method="nearest").isel(valid_time=0)

T_850_K = float(T_pt.sel(pressure_level=850))
T_500_K = float(T_pt.sel(pressure_level=500))

# Konversi ke Celsius
T_850 = T_850_K - 273.15
T_500 = T_500_K - 273.15

target_H = 3000.0  # ketinggian target (m)

# Susun array height dan temperature untuk interpolasi
heights = np.array([H_850, H_500])
temps   = np.array([T_850, T_500])

# Interpolasi linear dalam ruang ketinggian
f_interp = interp1d(heights, temps, kind="linear", bounds_error=True)
T_target = float(f_interp(target_H))

print("Suhu di level ERA5 (dekat Jakarta, 2024-01-01 00 UTC):")
print(f"  850 hPa  —  H = {H_850:.1f} m   →  T = {T_850:.2f} °C")
print(f"  500 hPa  —  H = {H_500:.1f} m   →  T = {T_500:.2f} °C")
print(f"\nInterpolasi ke {target_H:.0f} m:")
print(f"  T = {T_target:.2f} °C")

Suhu di level ERA5 (dekat Jakarta, 2024-01-01 00 UTC):
  850 hPa  —  H = 1512.1 m   →  T = 18.61 °C
  500 hPa  —  H = 5896.5 m   →  T = -3.84 °C

Interpolasi ke 3000 m:
  T = 10.99 °C

Nilai yang tercetak adalah estimasi suhu di 3000 m berdasarkan dua pressure levels terdekat. Di wilayah tropis seperti Indonesia, suhu di ketinggian sekitar 3000 m biasanya berada di rentang 5–12 °C tergantung kondisi sinoptik — jauh lebih dingin dari permukaan namun masih di atas suhu beku. Akurasi interpolasi meningkat bila kita menambahkan pressure level 700 hPa, karena asumsi linearitas berlaku di interval ketinggian yang lebih pendek.

Visualisasi Profil Vertikal

Profil vertikal dalam bentuk grafik mempermudah kita memverifikasi bahwa interpolasi berjalan masuk akal — apakah titik hasil interpolasi tepat berada pada garis yang menghubungkan dua data point ERA5, dan apakah gradien suhu terlihat realistis untuk kondisi tropis.

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(9, 10))

# Garis profil dan data point ERA5
ax.plot(temps, heights, "o-", color="#1E5A8A", linewidth=2.5, markersize=10,
        label="Data ERA5 (850 & 500 hPa)")

# Titik interpolasi
ax.plot(T_target, target_H, "*", color="crimson", markersize=18, zorder=5,
        label=f"Interpolasi ke {target_H:.0f} m: {T_target:.1f} °C")

# Garis referensi horizontal
ax.axhline(H_850,    color="#aaa", linestyle="--", linewidth=1.2,
           label=f"850 hPa (H = {H_850:.0f} m)")
ax.axhline(H_500,    color="#888", linestyle="--", linewidth=1.2,
           label=f"500 hPa (H = {H_500:.0f} m)")
ax.axhline(target_H, color="crimson", linestyle=":", linewidth=1.8,
           alpha=0.75, label="Target 3000 m")

ax.set_xlabel("Suhu (°C)", fontsize=13)
ax.set_ylabel("Geopotential Height (m)", fontsize=13)
ax.set_title(
    "Profil Suhu Vertikal ERA5\n"
    "Dekat Jakarta — 2024-01-01 00 UTC",
    fontsize=14,
)
ax.legend(fontsize=11, loc="lower left")
ax.grid(True, alpha=0.35)
fig.tight_layout()
fig.savefig("profil_suhu_vertikal.png", dpi=150, bbox_inches="tight")
print("Plot tersimpan: profil_suhu_vertikal.png")

snippet-4

Pada grafik di atas, sumbu \(y\) menampilkan geopotential height (bukan tekanan), sehingga jarak vertikal antar level sebanding dengan ketinggian fisik di atmosfer. Garis biru menghubungkan dua data point ERA5, dan bintang merah menandai hasil interpolasi di 3000 m — tepat berada pada garis interpolasi linier, seperti yang diharapkan. Gradien suhu yang terlihat (suhu menurun dengan ketinggian) konsisten dengan lapse rate normal di troposfer tropis.

Langkah Selanjutnya

Teknik yang kita demonstrasikan di sini bersifat modular dan bisa diperluas ke skenario yang lebih kompleks:

Tambah lebih banyak pressure levels. Dengan menambahkan level 700, 600, atau 750 hPa di request cdsapi, kita mendapatkan profil yang lebih padat. Interpolasi linier di antara level yang lebih rapat memberi hasil lebih akurat karena asumsi linearitas berlaku di interval yang jauh lebih pendek.

Terapkan ke seluruh grid spasial. Kode di atas hanya memproses satu grid point. Untuk menginterpolasi seluruh domain Indonesia sekaligus, gunakan xarray.apply_ufunc atau loop sederhana di atas dimensi latitude-longitude — hasilnya adalah field suhu di ketinggian tetap untuk seluruh wilayah.

Coba interpolasi log-tekanan. Ganti heights dengan np.log(np.array([850, 500])) sebagai sumbu vertikal dan sesuaikan target menjadi nilai log-tekanan yang setara. Metode ini lebih akurat secara fisik karena sebagian besar variabel atmosfer lebih linier terhadap \(\ln(p)\) di troposfer bebas.

Hitung indeks stabilitas. Setelah T, q, dan angin tersedia di ketinggian tetap, kita punya semua input untuk menghitung CAPE, CIN, dan indeks stabilitas lain — analisis yang tidak bisa dilakukan langsung dari pressure-level coordinates tanpa konversi ini.

Interpolasi vertikal ini juga merupakan komponen kunci dalam data assimilation di NWP, di mana observasi rawinsonde pada ketinggian geometris tertentu perlu dibandingkan dengan model first-guess di pressure coordinates. Memahami konversi geopotensial dan interpolasi ini adalah fondasi yang solid untuk melangkah ke analisis atmosfer yang lebih dalam.

Eksplorasi artikel meteorologi lainnya di meteo.my.id — kunjungi meteo.my.id untuk artikel selanjutnya.

Referensi

By di
Label: ,

Tidak ada komentar:

Posting Komentar

Langganan: Posting Komentar (Atom)