Plot deret waktu suhu udara 2 meter (T2m) ERA5 selama satu tahun untuk satu titik di Jakarta tampak seperti fluktuasi acak yang sulit dibaca. Tapi di balik "noise" itu tersembunyi sinyal periodik yang terstruktur — dengan siklus diurnal (24 jam) sebagai puncak yang jauh dominan, diikuti puncak semi-tahunan (~183 hari) yang mencerminkan fisika tropis. Fast Fourier Transform (FFT) mengubah domain waktu menjadi domain frekuensi sehingga setiap sinyal periodik muncul sebagai puncak tajam yang langsung bisa diidentifikasi.
Tutorial ini menerapkan numpy.fft.rfft dan scipy.signal.welch pada data ERA5 T2m 6-hourly tahun 2024 untuk wilayah Indonesia — dari download data di Copernicus CDS sampai memplot dan menafsirkan power spectrum dengan anotasi puncak diurnal dan semi-tahunan.
Mengenali Pola Temperatur Berkala
Setiap hari, permukaan bumi menyerap radiasi matahari saat siang dan memancarkan radiasi gelombang panjang secara kontinu. Selisih antara pemanasan siang dan pendinginan malam menciptakan osilasi 24-jam pada T2m yang kita sebut siklus diurnal. Menurut NASA, sebagian besar pemanasan matahari terjadi di permukaan sementara sebagian besar pendinginan radiatif terjadi di atmosfer — perbedaan inilah yang memicu naik-turunnya suhu permukaan secara periodik setiap harinya.
Di Indonesia — kawasan tropis dengan intensitas radiasi matahari yang konsisten sepanjang tahun — siklus diurnal sangat teratur dan amplitudonya bisa mencapai 5–10°C untuk titik di atas daratan. Di atas laut, kapasitas panas yang besar meredam osilasi ini hingga di bawah 1°C, sehingga grid point darat seperti Jakarta memberikan puncak diurnal yang jauh lebih tajam dibanding grid point oseanik.
Di kawasan tropis dalam seperti Indonesia, matahari melintas dekat zenith dua kali setahun (sekitar Maret dan September/Oktober), menghasilkan dua puncak insolation per tahun. Inilah yang memunculkan sinyal semi-tahunan (~183 hari) ketimbang siklus tahunan penuh dalam spektrum T2m; amplitudo musiman (<2°C) jauh lebih kecil dari amplitudo diurnal. FFT dapat memisahkan semua komponen periodik ini sekaligus.
Data yang kita gunakan: ERA5 reanalysis 6-hourly untuk Indonesia 2024, variabel 2 m temperature (T2m) pada grid \(0{,}25° \times 0{,}25°\). Kita ekstrak satu grid point terdekat ke Jakarta, sekitar \(-6{,}2°\)LS, \(106{,}8°\)BT. Kita download file satu kali dari CDS — daftar akun di cds.climate.copernicus.eu jika belum punya.
import os
import cdsapi
import xarray as xr
import numpy as np
OUT = "era5_t2m_indonesia_2024_6h.nc"
if not os.path.exists(OUT):
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,
)
ds = xr.open_dataset(OUT)
# ERA5 menyimpan variabel sebagai 't2m' atau '2m_temperature'
var_name = "t2m" if "t2m" in ds.data_vars else list(ds.data_vars)[0]
t2m_da = ds[var_name]
# Koordinat waktu bisa 'valid_time' atau 'time'
time_coord = "valid_time" if "valid_time" in t2m_da.dims else "time"
# Pilih grid point terdekat ke Jakarta
LAT_JKT, LON_JKT = -6.2, 106.8
t2m_jkt = t2m_da.sel(latitude=LAT_JKT, longitude=LON_JKT, method="nearest")
t2m_c = t2m_jkt.values - 273.15 # K → °C
lat_actual = float(t2m_jkt.latitude)
lon_actual = float(t2m_jkt.longitude)
N = len(t2m_c)
dt_day = 0.25 # 6 jam = 0,25 hari
t_first = str(t2m_jkt[time_coord].values[0])[:10]
t_last = str(t2m_jkt[time_coord].values[-1])[:10]
print(f"File : {OUT}")
print(f"Variabel : {var_name}")
print(f"Grid Jakarta: lat={lat_actual:.2f}°, lon={lon_actual:.2f}°")
print(f"N : {N} sampel (dt = 6 jam = {dt_day} hari)")
print(f"Periode : {t_first} – {t_last} UTC")
print(f"T2m : min={t2m_c.min():.1f}°C, mean={t2m_c.mean():.1f}°C, max={t2m_c.max():.1f}°C")
File : era5_t2m_indonesia_2024_6h.nc
Variabel : t2m
Grid Jakarta: lat=-6.25°, lon=106.75°
N : 1464 sampel (dt = 6 jam = 0.25 hari)
Periode : 2024-01-01 – 2024-12-31 UTC
T2m : min=21.6°C, mean=27.7°C, max=34.7°C
Output di atas mengonfirmasi jumlah sampel, resolusi temporal, dan rentang suhu pada grid point Jakarta. Nilai mean T2m sekitar 27–28°C tipikal untuk dataran rendah tropis. Fluktuasi di sekitar nilai mean itulah yang akan kita analisis secara spektral di langkah berikutnya.
Prinsip FFT dan Frekuensi Nyquist
FFT menguraikan deret waktu diskret dengan \(N\) sampel menjadi \(N/2 + 1\) komponen frekuensi unik (untuk input real, menggunakan rfft). Setiap komponen adalah bilangan kompleks; besarnya dikuadratkan menghasilkan power spectrum \(P(f) = |A(f)|^2\) dalam satuan K².
Frekuensi sampling data ERA5 6-hourly adalah \(f_s = 4\) cycles/day. Maka batas frekuensi tertinggi yang bisa direpresentasikan — Nyquist frequency — adalah:
$$f_{Nyquist} = \frac{f_s}{2} = 2\ \text{cycles/day}$$
yang setara dengan periode minimum \(T_{min} = 1/f_{Nyquist} = 0{,}5\) hari = 12 jam. Ini berarti siklus semi-diurnal (periode 12 jam, harmonik kedua siklus diurnal) jatuh tepat di batas Nyquist — zona rawan di mana sinyal bisa alias dan tidak bisa diinterpretasi dengan andal. Ini adalah keterbatasan data 6-hourly: untuk mengakses harmonik di atas 1 cycles/day dengan andal, diperlukan ERA5 hourly dari CDS yang mendorong Nyquist ke 12 cycles/day.
Resolusi frekuensi bergantung pada panjang deret waktu total:
$$\Delta f = \frac{f_s}{N} = \frac{4}{N}\ \text{cycles/day}$$
Tahun 2024 adalah tahun kabisat (366 hari × 4 = 1464 sampel), sehingga \(\Delta f \approx 0{,}00273\) cycles/day. Siklus tahunan berada di frekuensi \(\approx 0{,}00273\) cycles/day (bin pertama non-DC), dan siklus diurnal di tepat 1 cycle/day (jauh dari tetangganya). Keduanya bisa diidentifikasi tanpa ambiguitas.
Komputasi Spektrum Daya Mentah
Sebelum FFT, kita kurangi mean deret waktu untuk membuang komponen DC (level rata-rata, frekuensi nol). Selanjutnya numpy.fft.rfftfreq menghasilkan frekuensi dalam cycles/day langsung jika kita umpankan d=dt_day (spasi sampel dalam hari). Power tiap bin dihitung sebagai \(|A|^2\).
Alur komputasi FFT: dari deret waktu T2m mentah hingga ranking puncak power spectrum.
# Variabel t2m_c, N, dt_day diperoleh dari snippet-1
# FFT: subtract mean untuk buang DC
fft_vals = np.fft.rfft(t2m_c - t2m_c.mean())
freqs = np.fft.rfftfreq(N, d=dt_day) # cycles/day
power = np.abs(fft_vals) ** 2 # K²
# Konversi frekuensi ke periode (hindari bagi-nol pada DC bin)
with np.errstate(divide="ignore"):
periods = np.where(freqs > 0, 1.0 / freqs, np.inf)
# Top-5 puncak power (lewati DC di indeks 0)
top5_idx = np.argsort(power[1:])[::-1][:5] + 1
print(f"N={N}, fs={1/dt_day:.1f} cyc/day, Nyquist={0.5/dt_day:.1f} cyc/day")
print(f"Δf = {freqs[1]:.5f} cyc/day → resolusi terpanjang: {periods[1]:.0f} hari")
print()
print(f"{'Rank':<6}{'Periode (hari)':<18}{'Freq (cyc/day)':<22}{'Power (K²)'}")
print("-" * 60)
for rank, idx in enumerate(top5_idx, 1):
print(f"{rank:<6}{periods[idx]:<18.2f}{freqs[idx]:<22.5f}{power[idx]:.1f}")
N=1464, fs=4.0 cyc/day, Nyquist=2.0 cyc/day
Δf = 0.00273 cyc/day → resolusi terpanjang: 366 hari
Rank Periode (hari) Freq (cyc/day) Power (K²)
------------------------------------------------------------
1 1.00 1.00000 4703261.5
2 0.50 2.00000 1750790.2
3 183.00 0.00546 157677.5
4 1.00 1.00273 96822.5
5 122.00 0.00820 53454.7
Output di atas menampilkan ranking power spectral berdasarkan data aktual. Puncak paling dominan jatuh di periode 1,00 hari — siklus diurnal mendominasi spektrum ini dengan selisih yang sangat besar (~4,7 juta K²). Puncak kedua di periode 0,50 hari (semi-diurnal) berada tepat di frekuensi Nyquist (2 cycles/day) dan tidak bisa diinterpretasi dengan andal karena aliasing menjadi risiko nyata di batas sampling ini. Puncak ketiga yang menarik ada di periode 183 hari — ini adalah sinyal semi-tahunan yang secara fisik masuk akal untuk kawasan tropis: matahari melintas dekat zenith dua kali setahun, menghasilkan dua insolation maksimum, bukan satu. Tidak ada puncak tahunan (~366 hari) yang menonjol dalam top-5.
Visualisasi Spektrum Daya dengan Anotasi
Tabel angka di atas informatif, tetapi plot spektrum daya vs periode dengan sumbu-x logaritmik jauh lebih intuitif: rentang dari 0,5 hari hingga 500 hari muat dalam satu panel, dan kedua puncak langsung terlihat sebagai tonjolan tajam di atas noise floor.
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# Gunakan variabel freqs, periods, power dari snippet-2
# Buang DC (indeks 0) dan batasi pada periode < 500 hari
mask = (freqs[1:] > 0) & (1.0 / freqs[1:] < 500)
freq_p = freqs[1:][mask]
period_p = 1.0 / freq_p
power_p = power[1:][mask]
# Temukan indeks puncak diurnal (~1 hari) dan semi-tahunan (~183 hari)
idx_diurnal = int(np.argmin(np.abs(period_p - 1.0)))
idx_semiannual = int(np.argmin(np.abs(period_p - 183.0)))
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(period_p, power_p, color="#1E5A8A", lw=0.8, alpha=0.85, label="Power spectrum (FFT raw)")
# Garis vertikal puncak
ax.axvline(period_p[idx_diurnal], color="#2CA02C", lw=1.8, ls="--",
label=f"Siklus diurnal (~{period_p[idx_diurnal]:.1f} hari)")
ax.axvline(period_p[idx_semiannual], color="#E25A1C", lw=1.8, ls="--",
label=f"Siklus semi-tahunan (~{period_p[idx_semiannual]:.0f} hari)")
# Anotasi teks dengan panah
ax.annotate(
"Diurnal\n(~1 hari)",
xy=(period_p[idx_diurnal], power_p[idx_diurnal]),
xytext=(3.0, power_p[idx_diurnal] * 0.70),
fontsize=10, color="#2CA02C",
arrowprops=dict(arrowstyle="->", color="#2CA02C", lw=1.2),
)
ax.annotate(
"Semi-tahunan\n(~183 hari)",
xy=(period_p[idx_semiannual], power_p[idx_semiannual]),
xytext=(40, power_p[idx_semiannual] * 3.0),
fontsize=10, color="#E25A1C",
arrowprops=dict(arrowstyle="->", color="#E25A1C", lw=1.2),
)
ax.set_xscale("log")
ax.set_xlabel("Periode (hari)", fontsize=12)
ax.set_ylabel("Power (K²)", fontsize=12)
ax.set_title("Spektrum Daya T2m ERA5 — Grid Point Jakarta, 2024", fontsize=13, weight="bold")
ax.legend(fontsize=10, loc="upper left")
ax.grid(True, which="both", ls=":", alpha=0.4)
ax.set_xlim([0.5, 500])
plt.tight_layout()
plt.savefig("spectrum_fft.png", dpi=150, bbox_inches="tight")
print("Saved: spectrum_fft.png")
Plot yang dihasilkan memperlihatkan betapa "bersihnya" puncak diurnal (~1 hari) di sisi kiri panel — mendominasi spektrum dengan margin yang sangat besar. Di sisi kanan, puncak semi-tahunan (~183 hari) terlihat lebih rendah namun tetap menonjol di atas noise floor. Tidak ada puncak tahunan (~366 hari) yang signifikan: hal ini konsisten secara fisik — amplitudo naik-turun suhu harian di dataran rendah tropis (3–8°C) jauh melebihi variasi musiman T2m di Indonesia yang relatif kecil (<2°C). Puncak semi-tahunan yang muncul justru mencerminkan dua kali lintasan matahari dekat zenith per tahun, ciri khas spektrum suhu kawasan tropis dalam.
Pemulusan dengan Metode Welch
Raw FFT punya satu kelemahan: spectral leakage — energi dari puncak tajam "bocor" ke bin frekuensi tetangga karena ketidaksinambungan di ujung deret waktu yang panjangnya terbatas. Bocoran ini mempersulit identifikasi puncak yang lebih lemah di dekat puncak dominan.
Metode Welch (scipy.signal.welch) mengatasinya dengan tiga langkah:
- Potong deret waktu menjadi segmen-segmen yang saling tumpang tindih (overlap 50%)
- Terapkan Hann window — fungsi taper yang memuluskan tepi setiap segmen ke nol agar tidak ada diskontinuitas buatan. Hann window adalah pilihan default yang direkomendasikan dengan trade-off yang seimbang antara estimasi power yang akurat dan leakage yang rendah
- Rata-ratakan periodogram setiap segmen menjadi satu power spectral density (PSD) dengan variance yang lebih rendah
Parameter detrend="constant" membuang mean tiap segmen sebelum FFT, mencegah spike DC mendominasi ujung bawah spektrum.
from scipy.signal import welch
# Gunakan t2m_c dan dt_day dari snippet-1; power dan periods dari snippet-2
fs_day = 1.0 / dt_day # 4,0 cycles/day
# Panjang segmen 365 titik; resolusi terpanjang Welch = nperseg/fs = 365/4 ≈ 91 hari
nperseg = 365
noverlap = nperseg // 2 # 50% overlap
f_w, psd_w = welch(
t2m_c - t2m_c.mean(),
fs=fs_day,
window="hann",
nperseg=nperseg,
noverlap=noverlap,
detrend="constant",
)
# Lewati DC (f=0)
valid_w = f_w > 0
f_w = f_w[valid_w]
psd_w = psd_w[valid_w]
period_w = 1.0 / f_w
# Top-3 puncak Welch
top3_w = np.argsort(psd_w)[::-1][:3]
print(f"Welch: nperseg={nperseg}, noverlap={noverlap}, window=hann")
print(f"Δf Welch = {f_w[1] - f_w[0]:.5f} cyc/day (vs FFT: {freqs[1]:.5f})")
print()
print(f"{'Rank':<6}{'Periode Welch (hari)':<24}{'Freq (cyc/day)':<20}{'PSD (K²/[cyc/day])'}")
print("-" * 66)
for rank, idx in enumerate(top3_w, 1):
print(f"{rank:<6}{period_w[idx]:<24.2f}{f_w[idx]:<20.5f}{psd_w[idx]:.4f}")
# Bandingkan puncak #1 dan #2: FFT raw vs Welch
fft_top2 = np.argsort(power[1:])[::-1][:2] + 1
print()
print("--- Perbandingan puncak dominan ---")
print(f"{'Metode':<12}{'Puncak 1 (hari)':<22}{'Puncak 2 (hari)'}")
print(f"{'FFT raw':<12}{periods[fft_top2[0]]:<22.1f}{periods[fft_top2[1]]:.1f}")
print(f"{'Welch':<12}{period_w[top3_w[0]]:<22.2f}{period_w[top3_w[1]]:.2f}")
Welch: nperseg=365, noverlap=182, window=hann
Δf Welch = 0.01096 cyc/day (vs FFT: 0.00273)
Rank Periode Welch (hari) Freq (cyc/day) PSD (K²/[cyc/day])
------------------------------------------------------------------
1 1.00 0.99726 264.7197
2 0.99 1.00822 139.1834
3 0.50 1.99452 78.3605
--- Perbandingan puncak dominan ---
Metode Puncak 1 (hari) Puncak 2 (hari)
FFT raw 1.0 0.5
Welch 1.00 0.99
Output di atas mengonfirmasi bahwa tiga puncak teratas Welch semuanya berada di band diurnal (1,00 hari, 0,99 hari, dan 0,50 hari) — konsisten dengan hasil FFT raw yang menempatkan diurnal sebagai puncak jauh dominan. Ini bukan kebetulan: dengan nperseg=365 dan fs=4 cyc/day, periode terpanjang yang bisa diresolvasi Welch adalah \(365/4 = 91{,}25\) hari — siklus tahunan (~366 hari) jatuh jauh di luar jendela resolusi ini. Pelajaran penting: panjang segmen (nperseg) menentukan batas bawah frekuensi yang bisa diresolvasi, bukan panjang deret waktu total. Untuk resolusi tahunan penuh, gunakan nperseg ≥ 4 × 366 ≈ 1464. Perhatikan perbedaan satuan: power FFT dalam K² (energi total), sedangkan PSD Welch dalam K²/(cyc/day) (densitas per unit frekuensi).
Satu trade-off penting: dengan nperseg=365, resolusi frekuensi Welch adalah \(f_s/nperseg \approx 0{,}011\) cyc/day — sekitar empat kali lebih kasar dari FFT penuh (\(\approx 0{,}00273\) cyc/day). Ini harga yang wajar dibayar untuk mendapatkan spektrum yang lebih halus dan variance yang lebih rendah, terutama berguna ketika kita mencari puncak lemah di dekat puncak dominan yang kuat.
Apa Selanjutnya
Lewat tutorial ini kita berhasil mengungkap sinyal periodik yang tersembunyi dalam deret waktu T2m ERA5 Jakarta 2024: siklus diurnal (24 jam) adalah puncak yang jauh mendominasi, diikuti puncak semi-tahunan (~183 hari) yang mencerminkan dua kali lintasan matahari dekat zenith di kawasan tropis. Siklus semi-diurnal (0,5 hari) muncul tepat di Nyquist dan harus ditafsirkan dengan hati-hati. Tidak ada puncak tahunan yang signifikan — konsisten dengan amplitudo musiman yang kecil (<2°C) di dataran rendah tropis. FFT mengungkap semua ini dalam hitungan milidetik, dan metode Welch mengonfirmasi dominasi band diurnal sekaligus mengajarkan trade-off antara resolusi frekuensi dan panjang segmen.
Beberapa ekstensi yang layak dicoba berikutnya:
- Multi-grid mapping: jalankan FFT untuk seluruh grid Indonesia dan petakan amplitudo diurnal (\(\sqrt{2 \cdot P_{diurnal}} / N\)) per titik. Titik daratan Sumatera dan Jawa biasanya menunjukkan amplitudo diurnal jauh lebih besar dibanding titik oseanik — kontras yang menarik untuk divisualisasikan dengan cartopy.
- Wavelet transform: FFT mengasumsikan stasioneritas.
pywtatauscipyCWT memungkinkan spektrum yang berevolusi dalam waktu — berguna untuk melihat apakah amplitudo diurnal berubah antar musim. - Variabel lain: tekanan permukaan laut, kecepatan angin, dan curah hujan ERA5 menyimpan sinyal periodik berbeda yang bisa dianalisis dengan pola kode yang sama persis.
- ERA5 hourly: data 6-hourly menempatkan siklus semi-diurnal (12 jam) tepat di Nyquist. ERA5 hourly dari CDS mendorong Nyquist ke 12 cycles/day dan membuka akses ke harmonik yang lebih tinggi.
Eksplorasi artikel meteorologi lainnya di meteo.my.id — dari analisis MJO, pola tekanan monsun, hingga tutorial cartopy dan xarray untuk data Indonesia.
Referensi
- Discrete Fourier Transform (numpy.fft) — NumPy Documentation — dokumentasi
rfftdanrfftfreq: untuk input real, FFT menghasilkan \(N/2+1\) komponen unik; power dihitung sebagai \(|A(f)|^2\). - scipy.signal.welch — SciPy Documentation — Welch's method: averaging periodogram berwindow Hann dengan overlap 50% menghasilkan PSD dengan variance lebih rendah dan leakage yang berkurang.
- ERA5 Hourly Data on Single Levels from 1940 to Present — Copernicus Climate Data Store — deskripsi dataset ERA5: grid \(0{,}25°\), 4D-Var data assimilation, lisensi CC-BY, DOI 10.24381/cds.adbb2d47.
- Climate and Earth's Energy Budget — NASA Earth Observatory — anggaran energi Bumi: pemanasan permukaan siang hari oleh radiasi surya dan emisi gelombang panjang kontinu sebagai mekanisme fisik siklus diurnal T2m.
- Annual Migration of the Tropical Rain Belt — NOAA Climate.gov — intensitas pemanasan matahari yang konsisten di kawasan tropis sebagai dasar fisik dominannya puncak diurnal dalam spektrum ERA5 Indonesia.
Tidak ada komentar:
Posting Komentar