Memproses Citra Satelit MODIS dengan xarray

Sumber: NASA Earth Observatory, Lauren Dauphin — Indonesian Fires Return in 2023

Pendahuluan

Musim kemarau 2023 meninggalkan jejak yang bisa dilihat dari luar angkasa. El Niño yang kuat mendorong kekeringan ekstrem di Kalimantan, memicu kebakaran lahan yang membakar lebih dari 267.900 hektare sepanjang Oktober 2023. Asapnya melintasi batas provinsi, mengganggu penerbangan, dan mendorong kualitas udara ke level berbahaya di kota-kota besar Kalimantan.

Yang merekam semua itu dari orbit adalah MODIS — Moderate Resolution Imaging Spectroradiometer, sensor yang sudah mengorbit sejak akhir 1999 dan kini menjadi salah satu sumber data paling produktif untuk monitoring lingkungan dan meteorologi tropis.

Tutorial ini menunjukkan cara kerja praktis: dari search dan download data MOD11A1 (Land Surface Temperature harian) menggunakan library earthaccess, membuka file format HDF-EOS2 dengan rioxarray, mengkonversi nilai raw ke suhu Celsius, hingga memvisualisasikan hasilnya sebagai peta raster.

Produk MODIS untuk Meteorologi dan Monitoring Lingkungan

MODIS dipasang pada dua satelit NASA: Terra (diluncurkan 18 Desember 1999, melintas ekuator pukul 10:30 local time saat descending) dan Aqua (diluncurkan 4 Mei 2002, pukul 13:30 saat ascending). Keduanya mengorbit pada ketinggian 705 km dengan swath selebar 2.330 km, menghasilkan cakupan global setiap 1–2 hari.

Sensor ini merekam 36 band spektral dari 0,4 hingga 14,4 µm, dengan resolusi spasial yang bervariasi: 250 m untuk band 1–2, 500 m untuk band 3–7, dan 1 km untuk band 8–36. Dari data ini, NASA memproduksi lebih dari 40 produk geofisika dengan prefix MOD (Terra), MYD (Aqua), dan MCD (kombinasi keduanya).

Untuk keperluan meteorologi dan monitoring lingkungan di Indonesia, beberapa produk paling relevan adalah:

• MOD11A1 — Land Surface Temperature (LST) harian resolusi 1 km. Wajib untuk monitoring panas permukaan, urban heat island, dan identifikasi daerah kering.
• MOD14A1 — Fire hotspots harian 1 km. Produk ini yang merekam kebakaran Kalimantan 2023.
• MOD04_L2 — Aerosol Optical Depth (AOD) resolusi 10 km. Proxy untuk sebaran asap dan partikulat.
• MOD13Q1 — NDVI 16-hari resolusi 250 m. Pemantauan kondisi vegetasi dan drought stress.
• MOD35_L2 — Cloud mask swath. Digunakan sebagai input QC di produk turunan lain.

Catatan penting soal distribusi data: LP DAAC Data Pool, platform distribusi tradisional MODIS, resmi dihentikan pada 30 Juni 2025. Untuk akses data sekarang, gunakan earthaccess (Python library resmi NASA Earthdata) atau portal AppEEARS.

Setup

Kita butuh beberapa library. rioxarray adalah ekstensi xarray yang menambahkan dukungan format raster geospasial (termasuk HDF-EOS2 lewat GDAL), sementara earthaccess menangani autentikasi dan download dari NASA Earthdata.

pip install earthaccess rioxarray xarray numpy matplotlib

Selain itu, kita perlu akun NASA Earthdata. Daftar gratis di urs.earthdata.nasa.gov — prosesnya sekitar dua menit, dan akun yang sama berlaku untuk semua produk NASA Earthdata.

Mencari dan Mengunduh Data MOD11A1 dengan earthaccess

earthaccess menyederhanakan workflow yang sebelumnya butuh navigasi manual di Earthdata Search. Tiga fungsi utama: login() untuk autentikasi, search_data() untuk query granul, dan download() untuk menarik file ke lokal.

Contoh berikut mencari granul MOD11A1 untuk wilayah Indonesia (bounding box dari barat Sumatra hingga timur Papua) selama September–Oktober 2023 — periode puncak kebakaran Kalimantan.

import earthaccess

earthaccess.login()

results = earthaccess.search_data(
    short_name="MOD11A1",
    version="061",
    temporal=("2023-09-01", "2023-10-31"),
    bounding_box=(95, -11, 141, 6)   # (lon_min, lat_min, lon_max, lat_max)
)

files = earthaccess.download(results[:3], local_path="./modis_data/")
print(f"Downloaded {len(files)} file(s)")

Snippet ini memerlukan koneksi internet dan akun Earthdata yang valid. earthaccess.login() membaca kredensial dari ~/.netrc atau meminta input interaktif. Setelah download selesai, file .hdf tersimpan di folder ./modis_data/.

Membuka File HDF-EOS2 dengan rioxarray

File MODIS menggunakan format HDF-EOS2 — varian HDF4 dengan metadata geospasial tambahan. rioxarray dapat membacanya langsung lewat GDAL, meskipun hasilnya berupa list of DataArray, bukan satu Dataset tunggal (karena setiap layer HDF-EOS2 diperlakukan sebagai subdataset terpisah).

import rioxarray as rxr
import xarray as xr

layers = rxr.open_rasterio(
    "modis_data/MOD11A1.A2023275.h27v08.061.hdf",
    masked=True
)

# layers adalah list DataArray, satu per HDF group/subdataset
lst_layer = layers[0]   # LST_Day_1km biasanya layer pertama
print(lst_layer.dims, lst_layer.shape)
print(lst_layer.attrs)

Tile h27v08 mencakup bagian selatan Sumatra dan barat Kalimantan — lokasi tepat untuk memonitor kebakaran 2023. Nama tile MODIS menggunakan grid sinusoidal global; untuk referensi: h28v08 mencakup tengah/timur Kalimantan, h27v09 mencakup Jawa, dan h28v09 mencakup Sulawesi.

Setelah open_rasterio(), kita perlu identifikasi layer yang benar dari lst_layer.attrs['long_name'] karena urutannya bisa berbeda antar granul. Layer LST siang (LST_Day_1km) dan malam (LST_Night_1km) keduanya ada di file yang sama.

Pemrosesan Data LST

Nilai raw dalam MOD11A1 disimpan sebagai uint16. Konversi ke suhu Celsius menggunakan rumus linear:

$$T(\text{°C}) = \text{raw} \times 0{,}02 - 273{,}15$$

Nilai raw = 0 adalah fill value (piksel berawan atau tidak valid) yang harus dimasking sebelum analisis.

Snippet berikut membuat data sintetik yang merepresentasikan karakteristik nilai raw MOD11A1 khas — nilai sekitar 14.950–15.600 raw setara suhu permukaan 26–39 °C di wilayah tropis dataran rendah (Kalimantan musim kemarau). Kita juga menyimulasikan sebagian piksel berawan (raw = 0) yang lazim muncul di citra MODIS wilayah tropis.

import numpy as np
import xarray as xr

rng = np.random.default_rng(42)

# Buat data sintetik 100x100 piksel menyerupai MOD11A1 raw uint16
# Suhu target 26–39 °C → raw ≈ (26+273.15)/0.02=14957 hingga (39+273.15)/0.02=15607
raw_data = rng.integers(14950, 15600, size=(100, 100), dtype=np.uint16)

# Simulasikan ~15% piksel berawan (fill value = 0)
cloud_mask = rng.random((100, 100)) < 0.15
raw_data[cloud_mask] = 0

# Bungkus ke xarray DataArray
lat = np.linspace(-1.0, -1.9, 100)   # contoh: sekitar Kalimantan Tengah
lon = np.linspace(113.5, 114.4, 100)
da_raw = xr.DataArray(
    raw_data,
    dims=["lat", "lon"],
    coords={"lat": lat, "lon": lon},
    attrs={"long_name": "LST_Day_1km (synthetic)", "units": "uint16 raw"}
)

# Konversi ke Celsius: T(°C) = raw * 0.02 - 273.15
# Masking piksel berawan (raw == 0) sebelum konversi
da_celsius = da_raw.where(da_raw != 0).astype(float) * 0.02 - 273.15

# Simpan ke variabel global untuk snippet-5
lst_celsius = da_celsius

# Statistik ringkas
valid = da_celsius.values[~np.isnan(da_celsius.values)]
print(f"Total piksel      : {da_raw.size}")
print(f"Piksel valid (non-cloud): {valid.size}")
print(f"Piksel berawan    : {da_raw.size - valid.size}")
print(f"LST min           : {valid.min():.2f} °C")
print(f"LST max           : {valid.max():.2f} °C")
print(f"LST mean          : {valid.mean():.2f} °C")
print(f"LST std           : {valid.std():.2f} °C")

Total piksel      : 10000
Piksel valid (non-cloud): 8511
Piksel berawan    : 1489
LST min           : 25.85 °C
LST max           : 38.83 °C
LST mean          : 32.29 °C
LST std           : 3.75 °C

Output di atas menampilkan statistik suhu permukaan dari data sintetik: jumlah piksel valid setelah cloud masking, rentang LST sekitar 26–39 °C, dan rata-rata sekitar 32 °C — representatif untuk dataran rendah Kalimantan pada musim kemarau. Piksel berawan (sekitar 15% dari total) tidak ikut dihitung dalam statistik ini.

Visualisasi LST pada Peta Raster

Setelah masking dan konversi, langkah terakhir adalah visualisasi sebagai peta raster. Colormap RdYlBu_r (merah untuk panas, biru untuk dingin) adalah pilihan konvensional untuk LST karena intuitif dan mudah dibaca.

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

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

# Ambil data dan koordinat dari DataArray
data = lst_celsius.values
lats = lst_celsius.coords["lat"].values
lons = lst_celsius.coords["lon"].values

im = ax.imshow(
    data,
    extent=[lons.min(), lons.max(), lats.min(), lats.max()],
    origin="upper",
    cmap="RdYlBu_r",
    aspect="auto",
    vmin=np.nanpercentile(data, 2),
    vmax=np.nanpercentile(data, 98),
    interpolation="nearest"
)

cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
cbar.set_label("LST (°C)", fontsize=12)

ax.set_title("LST Siang (MOD11A1) — Data Sintetik\nKalimantan Tengah, September 2023",
             fontsize=13, fontweight="bold")
ax.set_xlabel("Longitude (°E)", fontsize=11)
ax.set_ylabel("Latitude (°S)", fontsize=11)
ax.tick_params(labelsize=10)

# Tandai piksel berawan sebagai putih
cmap_obj = plt.cm.RdYlBu_r.copy()
cmap_obj.set_bad(color="lightgray")
im.set_cmap(cmap_obj)

fig.tight_layout()
fig.savefig("lst_modis_synthetic.png", dpi=150, bbox_inches="tight")
plt.close(fig)
print("Plot saved: lst_modis_synthetic.png")

Peta yang dihasilkan menampilkan distribusi suhu permukaan siang hari dengan piksel berawan ditandai abu-abu terang. Pada data nyata, pola spasial LST mencerminkan kombinasi tutupan lahan (hutan primer vs lahan terdegradasi), elevasi, dan kelembaban tanah — komponen penting untuk analisis drought dan fire risk.

Langkah Berikutnya

Artikel ini baru menyentuh satu hari satu tile. Potensi analisis yang lebih dalam dari data MODIS mencakup:

Temporal stacking dan anomali LST — gabungkan 30–60 granul MOD11A1 dari rentang beberapa bulan ke satu xarray.Dataset menggunakan xr.open_mfdataset(). Dari sana, kita bisa hitung anomali relatif terhadap climatology multi-tahun untuk deteksi heatwave atau drought onset.

Kombinasi MOD14A1 + MOD11A1 — overlay fire hotspot dengan LST untuk memahami apakah area dengan LST tinggi memang berkorelasi dengan kejadian api beberapa hari kemudian. Ini dasar dari sistem early-warning berbasis MODIS.

MOD04_L2 untuk monitoring kualitas udara — AOD dari MODIS merupakan proxy yang baik untuk konsentrasi PM2.5 selama periode kebakaran. Kombinasi dengan ground station BMKG atau sensor PurpleAir menghasilkan analisis yang lebih kuat.

Cloud streaming — earthaccess mendukung akses langsung ke file di cloud NASA Earthdata lewat fsspec tanpa perlu download lokal. Untuk dataset besar atau analisis multi-tahun, ini menghemat storage secara signifikan.

Eksplorasi artikel meteorologi lainnya di meteo.my.id — kunjungi meteo.my.id untuk teknik analisis data atmosfer lanjutan.

Referensi

• MODIS Instrument — NASA Earthdata — Deskripsi resmi sensor MODIS, karakteristik spektral, dan dua platform satelit Terra dan Aqua.
• earthaccess — NASA Earthdata — Dokumentasi resmi library Python untuk search, download, dan streaming data NASA Earthdata.
• LP DAAC to Discontinue Data Pool Distribution for MODIS Datasets — NASA Earthdata — Pengumuman resmi penghentian LP DAAC Data Pool per 30 Juni 2025.
• MODIS/Terra Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V061 — NASA Earthdata — Halaman katalog produk MOD11A1, termasuk spesifikasi algoritma dan metadata.
• Indonesian Fires Return in 2023 — NASA Earth Observatory — Laporan NASA Earth Observatory tentang kebakaran Kalimantan Oktober 2023 beserta citra MODIS Terra.