BAB 26: Hutan dalam Film Pendek: Deret Waktu, Perubahan, dan Klasifikasi dengan Skrip

Studi kasus: "Apa yang terjadi sejak dua tahun lalu?" #

Kepala Seksi menaruh dua puluh empat citra satelit KPH Contoh di meja Anda. Pertanyaannya bertubi-tubi. Di mana hutan hilang, kapan, dan seberapa luas? Mana yang terbakar, dan benarkah titik panas itu? Kalau semua dijawab lewat klik, Anda mengulang pekerjaan dua puluh empat kali. Dengan skrip, tiap pertanyaan menjadi satu berkas Python yang bisa dijalankan ulang kapan saja.

Konsep: deret waktu dalam tiga kalimat #

Deret waktu itu seperti buku lipat bergambar: tiap lembar hampir sama, dan gerakan baru terlihat saat lembarnya dibalik berurutan. Satu citra hanya memotret satu hari. Puluhan citra yang ditumpuk menurut tanggal memberi tiap piksel sebuah riwayat angka, dan riwayat itu bisa dibaca seperti detak jantung.

Lima istilah yang perlu Anda kenal sekarang:

• Deret waktu: kumpulan citra lokasi yang sama pada tanggal berbeda. Di bab ini ada 24 tanggal, dari 5 Januari 2024 sampai 25 Desember 2025.
• NDVI: angka kehijauan tumbuhan, dihitung dari cahaya merah dan inframerah dekat. Hutan sehat bernilai sekitar 0,85. Awan, tanah terbuka, dan air bernilai jauh lebih rendah.
• SCL: peta penanda yang menyertai tiap citra. Isinya kode kelas tiap piksel, misalnya vegetasi, air, atau awan.
• Fitur: angka ringkas yang menggambarkan perilaku satu piksel. Fitur menjadi bahan belajar bagi mesin.
• Validasi: pemeriksaan apakah peta hasil boleh dipercaya.

Bab ini panjang karena isinya empat belas skrip. Bacalah per kelompok, dan jangan terburu-buru ingin menghafal. Tiap kelompok punya Konsep dan gambar sendiri.

Tiap skrip bernomor Skrip 26.n. Di dalam berkasnya mungkin tertulis nomor lain, misalnya SKRIP 1.1, mengikuti penomoran Buku M3. Bab ini selalu menyebut nama berkas aslinya.

Persiapan: folder, data, dan cara menjalankan #

Semua skrip bab ini memakai satu paket data dan satu folder hasil. Ketiganya harus Anda atur sebelum menjalankan skrip pertama.

Bahan. Paket data paket-m3. Isinya 24 citra Sentinel-2 buatan (berkas _refl.tif berisi enam band, berkas _SCL.tif berisi kode kelas), folder acuan (kunci jawaban), folder vektor (titik panas, titik latih, kawasan hutan tanaman), dan daftar tanggal Daftar_Citra.csv. Semua data sintetis, jadi hasilnya bukan kondisi hutan nyata. Luas area 400 m kali 400 m, atau 40 kali 40 piksel, dengan piksel 10 m.

1. Buat satu folder untuk skrip, misalnya D:/Latihan_M3/skrip. Salin ke sana keempat belas skrip dan dua berkas pembantu (m3_umum.py, hutan_mini.py). Semuanya harus satu folder.
2. Letakkan paket-m3 di tempat yang Anda tahu, misalnya D:/Latihan_M3/paket-m3.
3. Tentukan folder hasil, misalnya D:/Latihan_M3/hasil. Skrip membuatnya sendiri bila belum ada.
4. Beri tahu skrip lokasi paket dan hasil lewat dua variabel lingkungan (pengaturan di komputer yang bisa dibaca program): M3_PAKET dan M3_HASIL. Cara kerjanya dibahas di Bab 8.

Dari Konsol Python QGIS (Bab 3), langkah 4 dan menjalankan skrip pertama kira-kira seperti ini. Ketik tiap blok, lalu Enter. Konsol ini belum saya coba di jendela QGIS; urutannya saya ambil dari pola Bab 3 [CEK].

import os, sys
os.environ["M3_PAKET"] = r"D:/Latihan_M3/paket-m3"
os.environ["M3_HASIL"] = r"D:/Latihan_M3/hasil"
sys.path.append(r"D:/Latihan_M3/skrip")
exec(open(r"D:/Latihan_M3/skrip/1_1_susun_ndvi.py", encoding="utf8").read())

Urutan itu penting. Kedua variabel harus diisi sebelum skrip pertama dijalankan, karena m3_umum.py membacanya saat pertama kali diimpor. Bila Python sudah terlanjur mengimpornya, mulai ulang Konsol.

Bila Anda menjalankan lewat terminal, pakai pola "menyalakan QGIS tanpa jendela" dari Bab 3. Semua hasil di bab ini saya uji dengan pola itu, pada QGIS 4.0.2, dengan paket data di folder kerja saya. Karena itu, jalur yang muncul di blok HASIL saya tuliskan sebagai D:/Latihan_M3/hasil. Jalur di layar Anda akan mengikuti folder Anda sendiri.

Urutan menjalankan. Beberapa skrip memerlukan hasil skrip lain. Jalankan sesuai urutan di tabel ini, dan baca baris pertama keluarannya sebelum lanjut.

Dua kebiasaan akan menyelamatkan Anda sepanjang bab ini. Pertama, baca baris pertama keluaran. Kedua, jangan percaya skrip yang tidak berteriak: skrip yang tidak dijaga bisa tetap jalan dengan masukan yang salah dan menghasilkan angka yang keliru tanpa pesan galat. Skrip 26.4 memberi contoh penjagaannya: ia berhenti dengan pesan bila masukannya tidak terbaca. Bagian Jebakan Skrip 26.4 menjelaskannya.

Kotak perkakas: m3_umum.py #

Empat belas skrip bab ini memulai dengan baris import m3_umum as U. Berkas m3_umum.py adalah kotak perkakas buatan penulis: kumpulan fungsi kecil yang dipakai bersama, supaya tiap skrip tidak mengulang pekerjaan yang sama. Bayangkan bengkel dengan satu lemari alat. Semua tukang mengambil obeng dari lemari itu, dan tidak ada yang membuat obeng sendiri.

Kata U hanyalah nama pendek (alias) untuk kotak itu. U.daftar_citra() berarti: "ambil perkakas daftar_citra dari kotak U". Cara kerja import dan alias dibahas di Bab 8, dan cara membuat fungsi sendiri di Bab 7.

Tujuan. Menyediakan jalur data, daftar citra, pembaca dan penulis GeoTIFF, rumus NDVI dan NBR, serta penambal lubang awan, untuk dipakai semua skrip.

Bahan. Berkas ini tidak dijalankan sendiri. Ia diimpor oleh skrip lain dan harus satu folder dengan mereka. Ia memakai numpy dan osgeo.gdal (dibahas di Bab 12).

Skrip utuh:

# -*- coding: utf-8 -*-
"""m3_umum.py: fungsi bantu bersama untuk skrip Seri M3. Penulis: Badar Mubarok Yogaswara
Letakkan berkas ini satu folder dengan skrip lain. Ubah PAKET dan HASIL sesuai komputer Anda,
atau isi variabel lingkungan M3_PAKET dan M3_HASIL."""
import csv
import datetime as dt
import os

import numpy as np
from osgeo import gdal

gdal.UseExceptions()

PAKET = os.environ.get("M3_PAKET", r"D:/Latihan_M3/paket-m3")     # folder paket data sintetis
HASIL = os.environ.get("M3_HASIL", r"D:/Latihan_M3/hasil")        # folder keluaran Anda
os.makedirs(HASIL, exist_ok=True)

BAND = {"B02": 0, "B03": 1, "B04": 2, "B08": 3, "B11": 4, "B12": 5}   # urutan band di berkas *_refl.tif
BAIK = (4, 5, 6)                                                       # kode SCL: vegetasi, bukan vegetasi, air


def daftar_citra():
    """Baca Daftar_Citra.csv: kembalikan daftar (tanggal, jalur_refl, jalur_scl)."""
    hasil = []
    with open(os.path.join(PAKET, "citra", "Daftar_Citra.csv"), encoding="utf8") as f:
        for r in csv.DictReader(f):
            hasil.append((dt.date.fromisoformat(r["tanggal"]),
                          os.path.join(PAKET, "citra", r["berkas_refl"]),
                          os.path.join(PAKET, "citra", r["berkas_scl"])))
    return hasil


def baca_deret():
    """Kembalikan tanggal (daftar), refl (T,6,H,W), scl (T,H,W), geotransform, proyeksi."""
    tgl, refl, scl = [], [], []
    for d, pr, ps in daftar_citra():
        a = gdal.Open(pr)
        refl.append(a.ReadAsArray().astype("float32"))
        scl.append(gdal.Open(ps).ReadAsArray())
        gt, prj = a.GetGeoTransform(), a.GetProjection()
        tgl.append(d)
    return tgl, np.array(refl), np.array(scl), gt, prj


def nilai_ndvi(refl):
    nir, red = refl[:, BAND["B08"]], refl[:, BAND["B04"]]
    return (nir - red) / (nir + red + 1e-9)


def nilai_nbr(refl):
    nir, swir = refl[:, BAND["B08"]], refl[:, BAND["B12"]]
    return (nir - swir) / (nir + swir + 1e-9)


def bersih(scl):
    """True di piksel yang boleh dipakai (bukan awan, bayangan, atau cirrus)."""
    return np.isin(scl, BAIK)


def tulis_tif(path, arr, gt, prj, nodata=None, deskripsi=None, tipe=None):
    """Tulis larik (H,W) atau (B,H,W) ke GeoTIFF."""
    arr = np.asarray(arr)
    if arr.ndim == 2:
        arr = arr[None]
    if tipe is None:
        tipe = gdal.GDT_Float32 if arr.dtype.kind == "f" else gdal.GDT_Byte
    os.makedirs(os.path.dirname(path), exist_ok=True)
    ds = gdal.GetDriverByName("GTiff").Create(path, arr.shape[2], arr.shape[1], arr.shape[0], tipe, ["COMPRESS=LZW"])
    ds.SetGeoTransform(gt)
    ds.SetProjection(prj)
    for i in range(arr.shape[0]):
        b = ds.GetRasterBand(i + 1)
        b.WriteArray(arr[i])
        if nodata is not None:
            b.SetNoDataValue(nodata)
        if deskripsi:
            b.SetDescription(deskripsi[i])
    ds = None


def baca_tif(path):
    ds = gdal.Open(path)
    return ds.ReadAsArray(), ds.GetGeoTransform(), ds.GetProjection()


def piksel_ke_xy(gt, baris, kolom):
    return gt[0] + (kolom + 0.5) * gt[1], gt[3] + (baris + 0.5) * gt[5]


def xy_ke_piksel(gt, x, y):
    return int((gt[3] - y) / -gt[5]), int((x - gt[0]) / gt[1])


def isi_waktu(x, baik, hari):
    """Interpolasi linear sepanjang waktu. x dan baik berbentuk (T,H,W); hari = hari ke- tiap tanggal."""
    T, H, W = x.shape
    hasil = np.full(x.shape, np.nan, dtype="float64")
    for i in range(H):
        for j in range(W):
            ok = baik[:, i, j]
            if ok.any():
                hasil[:, i, j] = np.interp(hari, hari[ok], x[ok, i, j])
    return hasil

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Perkakas ini bisa Anda coba langsung di Konsol, setelah dua variabel lingkungan terisi:

import m3_umum as U
d = U.daftar_citra()
print(len(d), d[0][0])
tgl, refl, scl, gt, prj = U.baca_deret()
print(refl.shape, scl.shape)
print(U.xy_ke_piksel(gt, 312100, 9996100))

Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

24 2024-01-05
(24, 6, 40, 40) (24, 40, 40)
(30, 10)

Artinya: ada 24 tanggal, tiap tanggal enam band berukuran 40 kali 40, dan titik 312100, 9996100 jatuh di baris 30, kolom 10.

Bagian yang boleh diubah.

• Baris 14-15, jalur cadangan. Ganti dengan folder di komputer Anda, supaya tidak perlu mengisi variabel lingkungan tiap kali.
• Baris 19, BAIK. Tuple kelas SCL yang dianggap bersih. Menghapus angka 6 membuang piksel air dari analisis. Skrip 26.4 punya ekspresi sendiri, jadi ubahan ini tidak memengaruhinya.
• Baris 68, COMPRESS=LZW. Pilihan kompresi berkas hasil. Hapus bila Anda ingin berkas tanpa kompresi.

Jebakan.

• Folder dibuat saat impor. Baris 16 berjalan begitu import m3_umum dilakukan. Saya mengujinya: dengan M3_HASIL menunjuk folder yang belum ada, folder itu langsung muncul setelah impor, walau belum ada hasil apa pun. Salah ketik jalur akan menghasilkan folder baru yang salah.
• Tipe data diam-diam berubah. Baris 65-66 hanya mengenal Float32 dan Byte. Larik bilangan bulat besar terpotong tanpa pesan, kecuali Anda menyebut tipe sendiri. Kasus nyatanya ada di Skrip 26.9, yang memilih gdal.GDT_Int32. Buktinya saya tunjukkan di sana.
• Jalur cadangan menunjuk komputer penulis. Bila Anda lupa mengisi M3_PAKET, galatnya FileNotFoundError saat membuka Daftar_Citra.csv. Saya mengujinya dengan jalur paket yang salah, dan galat itu muncul dari baris 25.
• Semua citra dianggap sebidang. Baris 40 hanya mengingat geotransform citra terakhir. Bila citra Anda berbeda ukuran atau lokasi, hasilnya keliru tanpa peringatan.

Latihan variasi.

1. Cetak tanggal dan jalur berkas SCL untuk tanggal terakhir.
2. Hitung berapa persen piksel bersih (menurut bersih) pada seluruh deret.

Contoh jawaban:

print(d[-1][0], d[-1][2])
print(round(float(U.bersih(scl).mean()), 3))

Baris pertama mencetak tanggal terakhir dan jalur SCL-nya. Baris kedua mencetak 0.747 pada data ini: sekitar 75 persen piksel-tanggal bersih, sisanya (25 persen) tertutup awan, bayangan, atau cirrus.

Kelompok A: Menyusun deret, dari 24 berkas menjadi satu tumpukan #

Konsep: tumpukan bertanggal dalam tiga kalimat #

Bayangkan 24 foto yang dijilid menjadi album bertanggal: satu halaman satu tanggal. Di komputer, album itu disebut tumpukan (stack): satu berkas raster dengan 24 band, dan tiap band adalah satu tanggal. Bila Anda menusuk album itu dengan jarum di satu titik, nilai-nilai yang tertusuk membentuk profil waktu titik itu.

Lima istilah yang perlu Anda kenal sekarang:

• Tumpukan (stack): satu berkas raster yang berisi banyak band.
• Band: satu lembar dalam tumpukan. Di bab ini, satu band berarti satu tanggal.
• Profil waktu: deretan nilai satu piksel sepanjang tanggal.
• Kalkulator raster (Raster Calculator): alat QGIS yang menghitung rumus di atas raster, piksel demi piksel.
• processing.run: perintah Python untuk memanggil alat QGIS, seperti menekan tombol alat itu dari skrip (dibahas di Bab 21).

Kelompok ini terdiri dari tiga skrip. Skrip 26.1 membuat tumpukannya. Skrip 26.2 membacanya dan menggambar grafik. Skrip 26.3 menempelkan tanggal supaya QGIS bisa memutarnya seperti film.

Skrip 26.1: Menyusun NDVI menjadi tumpukan 24 band #

Berkas 1_1_susun_ndvi.py. Di dalam berkasnya tertulis "SKRIP 1.1", mengikuti penomoran Buku M3.

Tujuan. Menghitung NDVI untuk tiap dari 24 citra, lalu menumpuk semuanya menjadi satu berkas 24 band yang tiap bandnya diberi tanggal.

Bahan. Paket data paket-m3 (folder citra) dan kotak perkakas m3_umum.py. Jalankan dari Konsol Python QGIS, atau dari terminal dengan pola Bab 3. Waktu jalan sekitar sepuluh detik.

Skrip ini bekerja dalam enam tahap.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 1.1: Susun deret waktu NDVI dari 24 citra Sentinel-2 sintetis (jalankan di Python Console QGIS)
# Penulis: Badar Mubarok Yogaswara
# Hasil: NDVI per tanggal (hasil/ndvi/NDVI_YYYYMMDD.tif) dan satu tumpukan 24 band: hasil/NDVI_deret.tif
import os
import processing
from qgis.core import QgsRasterLayer, QgsProject
from osgeo import gdal
import m3_umum as U

os.makedirs(os.path.join(U.HASIL, "ndvi"), exist_ok=True)
daftar = U.daftar_citra()
keluaran = []
for tgl, pr, ps in daftar:
    nama = "S2_" + tgl.strftime("%Y%m%d") + "_refl"
    lyr = QgsRasterLayer(pr, nama)
    QgsProject.instance().addMapLayer(lyr, False)            # tanpa menampilkan di kanvas
    out = os.path.join(U.HASIL, "ndvi", "NDVI_%s.tif" % tgl.strftime("%Y%m%d"))
    ekspresi = '("%s@4" - "%s@3") / ("%s@4" + "%s@3")' % (nama, nama, nama, nama)   # band 4 = B08, band 3 = B04
    processing.run("native:rastercalc", {"LAYERS": [lyr], "EXPRESSION": ekspresi, "OUTPUT": out})
    keluaran.append(out)

# gabungkan menjadi satu berkas: tiap berkas menjadi satu band (SEPARATE)
stack = os.path.join(U.HASIL, "NDVI_deret.tif")
processing.run("gdal:merge", {"INPUT": keluaran, "SEPARATE": True, "DATA_TYPE": 5, "OUTPUT": stack})
# beri nama band sesuai tanggal
ds = gdal.Open(stack, gdal.GA_Update)
for i, (tgl, _, _) in enumerate(daftar):
    ds.GetRasterBand(i + 1).SetDescription(tgl.isoformat())
ds = None
lyr = QgsRasterLayer(stack, "NDVI_deret")
print("Tumpukan:", stack, "valid:", lyr.isValid(), "band:", lyr.bandCount())

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Tumpukan: D:/Latihan_M3/hasil\NDVI_deret.tif valid: True band: 24

Dari hasil ini Anda tahu: tumpukan terbaca (valid: True) dan memuat 24 band. Folder ndvi kini berisi 24 berkas NDVI, satu per tanggal.

Bagian yang boleh diubah.

• Baris 19, nomor band. Ganti @3 menjadi @6 di seluruh baris untuk menghitung NBR (inframerah dekat dan SWIR2) alih-alih NDVI. Dengan angka 6, rata-rata band pertama berubah dari 0.726 menjadi 0.586, jadi indeksnya memang lain.
• Baris 25, DATA_TYPE. Angka 5 berarti Float32. Saya memeriksanya dengan bantuan alat (qgis_process help gdal:merge) dan dengan gdalinfo pada hasilnya. Angka lain memberi tipe lain, dan desimal NDVI bisa terpotong.
• Baris 11 dan 18, nama folder dan berkas. Ganti bila Anda ingin menyimpan di tempat lain.

Jebakan.

• Nama layer harus sama persis dengan nama di rumus. Saya mencoba mengganti satu nama di baris 19 dengan nama lain. Hasilnya galat QgsProcessingException: Error occurred while performing calculation. Pesannya tidak menyebut nama mana yang salah, jadi periksa ejaan di baris 15 dan 19.
• import processing hanya berjalan sendiri di Konsol QGIS. Saya mencobanya dari terminal biasa tanpa persiapan: layar menampilkan ModuleNotFoundError: No module named 'processing'. Nyalakan QGIS tanpa jendela lebih dulu (Bab 3).
• Menjalankan ulang menimpa berkas lama. Tidak ada pertanyaan "yakin?". Pindahkan hasil lama bila ingin menyimpannya.
• Kotak ps tidak dipakai. Python tidak protes terhadap kotak yang menganggur. Itu tanda skrip ini hanya memerlukan dua dari tiga isi tuple.

Latihan variasi.

1. Hitung NBR (inframerah dekat dan SWIR2) alih-alih NDVI.
2. Cetak juga nilai rata-rata band pertama tumpukan, untuk membandingkan dua indeks.

Contoh jawaban, hanya bagian yang berubah:

    ekspresi = '("%s@4" - "%s@6") / ("%s@4" + "%s@6")' % (nama, nama, nama, nama)   # band 6 = B12

dan satu baris tambahan di akhir skrip:

print("rata-rata band 1:", round(float(gdal.Open(stack).ReadAsArray()[0].mean()), 3))

Saya menjalankan versi NDVI dan versi NBR pada salinan data. Layar menampilkan rata-rata band 1: 0.726 untuk NDVI dan rata-rata band 1: 0.586 untuk NBR.

Skrip 26.2: Menggambar profil waktu empat titik #

Berkas 1_2_profil_waktu.py. Di dalam berkasnya tertulis "SKRIP 1.2".

Tujuan. Menggambar grafik NDVI sepanjang 24 tanggal untuk empat titik (hutan alam, sawah, bukaan, bekas bakar), mencetak ringkasannya, dan menyimpan grafik sebagai gambar.

Bahan. Hasil Skrip 26.1 (NDVI_deret.tif), m3_umum.py, dan matplotlib (perkakas membuat grafik, dibahas di Bab 12). Jalankan dari Konsol atau terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 1.2: Profil waktu satu titik (NDVI sepanjang 24 tanggal) dari tumpukan NDVI_deret.tif
# Penulis: Badar Mubarok Yogaswara
import os
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from qgis.core import QgsRasterLayer, QgsPointXY
import m3_umum as U

lyr = QgsRasterLayer(os.path.join(U.HASIL, "NDVI_deret.tif"), "NDVI_deret")
dp = lyr.dataProvider()
tanggal = [t for t, _, _ in U.daftar_citra()]
titik = {"Hutan alam": (312100, 9996100), "Sawah (Beta)": (312365, 9996025),
         "Bukaan (Gamma)": (312100, 9996330), "Bekas bakar (Beta)": (312270, 9996125)}
fig, ax = plt.subplots(figsize=(9, 4))
for nama, (x, y) in titik.items():
    nilai = []
    for b in range(1, lyr.bandCount() + 1):
        v, ok = dp.sample(QgsPointXY(x, y), b)
        nilai.append(v if ok else np.nan)
    print("%-20s NDVI awal %.2f, minimum %.2f, akhir %.2f" % (nama, nilai[0], np.nanmin(nilai), nilai[-1]))
    ax.plot(tanggal, nilai, marker="o", ms=3, label=nama)
ax.set_ylabel("NDVI")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.autofmt_xdate()
fig.savefig(os.path.join(U.HASIL, "profil_waktu.png"), dpi=130, bbox_inches="tight")
print("Grafik tersimpan di", os.path.join(U.HASIL, "profil_waktu.png"))

Penjelasan baris demi baris. Hanya hal baru yang dijelaskan. Perkakas yang sudah dikenal di Skrip 26.1 tidak diulang.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Hutan alam           NDVI awal 0.87, minimum 0.00, akhir 0.88
Sawah (Beta)         NDVI awal 0.74, minimum 0.06, akhir 0.72
Bukaan (Gamma)       NDVI awal 0.86, minimum 0.03, akhir 0.04
Bekas bakar (Beta)   NDVI awal 0.87, minimum 0.01, akhir 0.01
Grafik tersimpan di D:/Latihan_M3/hasil\profil_waktu.png

Berkas profil_waktu.png berisi grafik berikut. Ini keluaran sungguhan skrip, bukan skema.

Perhatikan garis biru, hutan alam. Ia hampir datar di sekitar 0,85, lalu jatuh ke nol pada beberapa tanggal. Hutannya tidak hilang. Pada tanggal itu, awan menutupi piksel. Itulah mengapa kelompok B ada.

Bagian yang boleh diubah.

• Baris 15-16, koordinat titik. Ganti pasangan angka dengan titik Anda sendiri, dalam sistem koordinat yang sama dengan raster. Tambah titik baru dengan menambah satu pasangan nama dan koordinat di kamus.
• Baris 17, figsize. Ubah ukuran gambar. (12, 5) membuat grafik lebih lebar.
• Baris 29, dpi. Naikkan menjadi 200 untuk gambar yang lebih tajam.

Jebakan.

• Titik di luar raster tidak menghasilkan galat. Saya mengganti koordinat hutan alam menjadi (0, 0). Layar menampilkan Hutan alam NDVI awal nan, minimum nan, akhir nan dan peringatan RuntimeWarning: All-NaN axis encountered. Skrip tetap selesai. Selalu baca angka ringkasan, bukan hanya gambarnya.
• Skrip 26.1 harus sudah dijalankan. Saya menjalankan skrip ini di folder hasil yang kosong. Layar menampilkan nan untuk titik pertama, lalu berhenti dengan ValueError: x and y must have same first dimension, but have shapes (24,) and (1,). Artinya tumpukan tidak ditemukan.
• Urutan baris 6-8 sengaja begitu. Aturan matplotlib: pilih mode lebih dulu, baru panggil pyplot. Jangan memindahkan baris 7 ke bawah baris 8.
• matplotlib harus terpasang di Python QGIS. Skrip tidak memeriksanya (Bab 12). Di QGIS 4.0.2 milik saya, matplotlib sudah ada.

Latihan variasi.

1. Tambahkan titik Air di koordinat (312295, 9996195).
2. Pastikan hasilnya masuk ke ringkasan dan grafik.

Contoh jawaban, hanya bagian yang berubah (baris 15-16):

titik = {"Hutan alam": (312100, 9996100), "Sawah (Beta)": (312365, 9996025),
         "Bukaan (Gamma)": (312100, 9996330), "Bekas bakar (Beta)": (312270, 9996125),
         "Air": (312295, 9996195)}

Layar menampilkan satu baris tambahan: Air NDVI awal -0.31, minimum -0.52, akhir 0.04. NDVI air negatif, tanda khas air.

Skrip 26.3: Memberi tanggal agar bisa diputar seperti film #

Berkas 1_3_atur_waktu.py. Di dalam berkasnya tertulis "SKRIP 1.3".

Tujuan. Memberi rentang tanggal pada tiap band NDVI_deret.tif supaya QGIS bisa memutarnya di panel Temporal Controller.

Bahan. Hasil Skrip 26.1 dan m3_umum.py. Skrip ini hanya bermakna di dalam proyek QGIS yang terbuka, jadi jalankan dari Konsol.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 1.3: Beri tanggal pada tiap band NDVI_deret.tif agar bisa dimainkan di Temporal Controller
# Penulis: Badar Mubarok Yogaswara
# Mode yang dipakai: Fixed Time Range Per Band (tiap band punya rentang tanggal sendiri).
import os
from qgis.core import QgsRasterLayer, QgsProject, QgsDateTimeRange, Qgis
from qgis.PyQt.QtCore import QDateTime, QDate, QTime
import m3_umum as U

lyr = QgsRasterLayer(os.path.join(U.HASIL, "NDVI_deret.tif"), "NDVI_deret")
rentang = {}
for b, (tgl, _, _) in enumerate(U.daftar_citra(), start=1):
    awal = QDateTime(QDate(tgl.year, tgl.month, tgl.day), QTime(0, 0, 0))
    rentang[b] = QgsDateTimeRange(awal, awal.addDays(1))          # satu hari untuk tiap band
tp = lyr.temporalProperties()
tp.setMode(Qgis.RasterTemporalMode.FixedRangePerBand)
tp.setFixedRangePerBand(rentang)
tp.setIsActive(True)
QgsProject.instance().addMapLayer(lyr)
print("Aktif:", tp.isActive(), "| mode:", tp.mode(), "| band bertanggal:", len(tp.fixedRangePerBand()))
print("Band 9 =", tp.fixedRangePerBand()[9].begin().toString("yyyy-MM-dd"))

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Aktif: True | mode: 3 | band bertanggal: 24
Band 9 = 2024-09-16

Band 9 bertanggal 2024-09-16, sesuai urutan di Daftar_Citra.csv. Angka mode: 3 adalah nomor mode FixedRangePerBand. Saya memeriksanya dengan int(Qgis.RasterTemporalMode.FixedRangePerBand), yang bernilai 3. Memutar layer di panel Temporal Controller saya tandai [CEK]: pengujian saya berjalan tanpa jendela QGIS.

Bagian yang boleh diubah.

• Baris 14, addDays(1). Panjang rentang tiap band. Saya mengganti menjadi addDays(10), dan tanggal akhir band 9 berubah dari 2024-09-17 menjadi 2024-09-26. Rentang yang lebih panjang membuat band bertahan lebih lama saat film diputar.
• Baris 10, nama layer. Ubah bila Anda memakai tumpukan lain.

Jebakan.

• Pengaturan hidup di proyek, bukan di berkas. Tanggal ini tidak tertulis ke NDVI_deret.tif. Simpan proyek QGIS (Bab 14) bila ingin menyimpannya.
• Menjalankan dua kali menambah layer kembar. Saya mengujinya: setelah dua kali jalan, proyek berisi dua layer bernama NDVI_deret. Hapus yang lama sebelum mengulang.
• Nama Qgis.RasterTemporalMode mengikuti gaya QGIS baru. Saya mengujinya hanya di QGIS 4.0.2. Di versi lain, nama ini mungkin berbeda [CEK].

Latihan variasi.

1. Cetak juga tanggal akhir band ke-9.
2. Buat tiap band berlaku 10 hari, lalu cetak tanggal akhir band ke-9 lagi.

Contoh jawaban, tambahkan di akhir skrip:

print("Band 9 berakhir:", tp.fixedRangePerBand()[9].end().toString("yyyy-MM-dd"))

Dengan addDays(1), layar menampilkan Band 9 berakhir: 2024-09-17. Setelah baris 14 diganti menjadi addDays(10), layar menampilkan Band 9 berakhir: 2024-09-26.

Kelompok B: Membersihkan awan #

Konsep: tiga cara menyingkirkan awan #

Awan itu seperti kaca buram di depan kamera: apa pun yang tampak di baliknya bukan permukaan bumi. Anda punya tiga cara. Pertama, membuang piksel yang buram, memakai peta penanda. Kedua, menggabungkan beberapa tanggal menjadi satu citra bersih. Ketiga, menambal lubangnya dengan nilai dari tanggal sebelum dan sesudah. Ketiganya bisa dipakai bersama.

Lima istilah yang perlu Anda kenal sekarang:

• Mask: peta penanda piksel mana yang dipakai dan mana yang dibuang.
• NoData: nilai khusus yang berarti "piksel ini kosong". Piksel kosong tidak ikut dihitung.
• Median: nilai tengah bila semua nilai diurutkan. Lawan kata rata-rata.
• Komposit median: citra baru yang tiap pikselnya berisi median dari beberapa tanggal.
• Interpolasi linear: menarik garis lurus antara dua nilai yang diketahui, untuk menebak nilai di antaranya.

Median lebih tahan banting daripada rata-rata. Satu nilai ekstrem, misalnya awan yang lolos dari penyaringan, tidak ikut menarik hasilnya.

Skrip 26.4: Membuang awan dengan band SCL #

Berkas 2_1_masker_awan.py. Di dalam berkasnya tertulis "SKRIP 2.1".

Tujuan. Membuang piksel berawan, berbayang awan, atau berkabut tipis dari tiap NDVI memakai band SCL. Piksel buruk menjadi NoData.

Bahan. Hasil Skrip 26.1 (folder ndvi), berkas _SCL.tif dari paket data, dan m3_umum.py. Jalankan dari Konsol atau terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 2.1: Buang awan dan bayangan memakai band SCL, per tanggal, dengan Raster Calculator
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 1.1 sudah dijalankan (hasil/ndvi/NDVI_*.tif ada). Piksel buruk menjadi NoData.
import os
import processing
from qgis.core import QgsRasterLayer, QgsProject
from osgeo import gdal
import numpy as np
import m3_umum as U

os.makedirs(os.path.join(U.HASIL, "ndvi_bersih"), exist_ok=True)
daftar = U.daftar_citra()
keluaran = []
for tgl, pr, ps in daftar:
    k = tgl.strftime("%Y%m%d")
    n_ndvi, n_scl = "NDVI_" + k, "SCL_" + k
    l_ndvi = QgsRasterLayer(os.path.join(U.HASIL, "ndvi", "NDVI_%s.tif" % k), n_ndvi)
    l_scl = QgsRasterLayer(ps, n_scl)
    if not l_ndvi.isValid() or not l_scl.isValid():
        raise SystemExit("Layer tidak terbaca (jalankan Skrip 1.1 dulu): " + n_ndvi)
    QgsProject.instance().addMapLayers([l_ndvi, l_scl], False)
    # piksel baik: SCL 4 (vegetasi), 5 (bukan vegetasi), 6 (air). Nilai baik dibagi 1; nilai buruk dibagi 0 -> NoData
    ekspresi = '"%s@1" / (("%s@1" = 4) OR ("%s@1" = 5) OR ("%s@1" = 6))' % (n_ndvi, n_scl, n_scl, n_scl)
    out = os.path.join(U.HASIL, "ndvi_bersih", "NDVI_bersih_%s.tif" % k)
    processing.run("native:rastercalc", {"LAYERS": [l_ndvi, l_scl], "EXPRESSION": ekspresi, "OUTPUT": out})
    keluaran.append(out)
    a = gdal.Open(out)
    arr = a.GetRasterBand(1).ReadAsArray()
    nd = a.GetRasterBand(1).GetNoDataValue()
    valid = np.isfinite(arr) if nd is None else (arr != nd) & np.isfinite(arr)
    print(tgl, "piksel valid: %5.1f%%" % (100 * valid.mean()))
print("selesai,", len(keluaran), "berkas")

Penjelasan baris demi baris. Hanya hal baru yang diuraikan lebih panjang.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

2024-01-05 piksel valid:  93.5%
2024-02-04 piksel valid:  61.5%
2024-03-10 piksel valid:  41.6%
2024-04-09 piksel valid:  70.3%
2024-05-09 piksel valid:  81.3%
2024-06-13 piksel valid: 100.0%
2024-07-13 piksel valid: 100.0%
2024-08-12 piksel valid:  93.2%
2024-09-16 piksel valid:  66.9%
2024-10-16 piksel valid:  90.5%
2024-11-15 piksel valid:  75.8%
2024-12-20 piksel valid:  52.1%
2025-01-14 piksel valid:  46.1%
2025-02-13 piksel valid:  67.4%
2025-03-15 piksel valid:  51.0%
2025-04-14 piksel valid:  64.3%
2025-05-19 piksel valid:  88.6%
2025-06-18 piksel valid: 100.0%
2025-07-18 piksel valid:  93.2%
2025-08-17 piksel valid:  83.3%
2025-09-21 piksel valid: 100.0%
2025-10-21 piksel valid:  75.6%
2025-11-20 piksel valid:  57.9%
2025-12-25 piksel valid:  39.5%
selesai, 24 berkas

Bacalah angka-angkanya. Pada musim kemarau (Juni dan Juli 2024), 100 persen piksel valid. Pada Desember 2025, hanya 39,5 persen. Tiap tanggal punya porsi awan sendiri, dan itulah yang akan Anda atasi di dua skrip berikutnya.

Bagian yang boleh diubah.

• Baris 24, kode SCL. Hapus bagian OR ("%s@1" = 6) untuk membuang piksel air, dan kurangi satu n_scl di tuple isian. Saya mencobanya: pada 5 Januari 2024, persen valid turun dari 93.5 menjadi 89.4.
• Baris 12 dan 25, nama folder dan berkas. Bebas Anda ganti.
• Baris 21, teks pesan. Bebas Anda ubah, misalnya dengan menyebut folder yang Anda pakai. Pesan itu hanya tampil bila ada layer yang gagal dibuka.

Jebakan.

• Pengaman di baris 20-21 menjaga urutan kerja. Saya menjalankan skrip di folder hasil yang kosong, sebelum Skrip 26.1. Skrip langsung berhenti dengan pesan Layer tidak terbaca (jalankan Skrip 1.1 dulu): NDVI_20240105 dan tidak menulis satu berkas pun (hanya folder kosong ndvi_bersih yang sudah sempat dibuat). Di berkasnya, Skrip 26.1 tertulis "Skrip 1.1". Tanpa dua baris itu, skrip tetap selesai tanpa galat dan menulis 24 berkas sampah. Pengaman ini hanya memeriksa bahwa layer bisa dibuka, bukan bahwa isinya benar. Tetap cocokkan angka keluaran dengan angka yang seharusnya.
• Trik "bagi nol menjadi NoData" bergantung pada kalkulator raster QGIS. Saya memeriksanya pada 10 Maret 2024: ada 935 piksel dengan kode SCL buruk, dan tepat 935 piksel NoData di hasil, di posisi yang sama. Pada QGIS 4.0.2 trik ini bekerja. Pada versi lain, periksa dulu [CEK].
• Nilai NoData bukan -9999. Berkas hasil memakai bilangan terkecil Float32 (sekitar -3.4e+38). Karena itu skrip mengambil nd dari berkas, tidak menebaknya.
• SCL adalah tebakan otomatis. Tepi awan yang sangat tipis bisa lolos sebagai "bersih". Itu sebabnya ada komposit median di skrip berikutnya.

Latihan variasi.

1. Hanya terima piksel vegetasi dan bukan vegetasi (kode 4 dan 5). Air dibuang.
2. Lihat persen valid pada tiga tanggal pertama.

Contoh jawaban, hanya baris yang berubah (baris 24):

    ekspresi = '"%s@1" / (("%s@1" = 4) OR ("%s@1" = 5))' % (n_ndvi, n_scl, n_scl)

Layar menampilkan tiga tanggal pertama: 2024-01-05 piksel valid: 89.4%, 2024-02-04 piksel valid: 59.7%, 2024-03-10 piksel valid: 41.2%. Semuanya lebih rendah daripada hasil asli (93,5; 61,5; 41,6), karena piksel air ikut dibuang.

Skrip 26.5: Komposit median dari beberapa tanggal #

Berkas 2_2_komposit_median.py. Di dalam berkasnya tertulis "SKRIP 2.2".

Tujuan. Membuat citra komposit median NDVI dan peta jumlah pengamatan valid untuk tiga periode: sebelum, kemarau 2024, dan sesudah.

Bahan. Hasil Skrip 26.4 (folder ndvi_bersih) dan m3_umum.py. Jalankan dari Konsol atau terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 2.2: Komposit median dan jumlah pengamatan valid dengan Cell Statistics
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 2.1 sudah dijalankan.
import os
import numpy as np
import processing
from osgeo import gdal
import m3_umum as U

daftar = U.daftar_citra()
berkas = [os.path.join(U.HASIL, "ndvi_bersih", "NDVI_bersih_%s.tif" % t.strftime("%Y%m%d")) for t, _, _ in daftar]
periode = {"sebelum": (0, 8), "kemarau2024": (5, 9), "sesudah": (20, 24)}   # indeks [awal, akhir) di daftar tanggal
for nama, (a, b) in periode.items():
    for kode, stat in (("median", 3), ("jumlah", 1)):
        out = os.path.join(U.HASIL, "komposit_ndvi_%s_%s.tif" % (nama, kode))
        processing.run("native:cellstatistics", {"INPUT": berkas[a:b], "STATISTIC": stat, "IGNORE_NODATA": True,
                                                  "REFERENCE_LAYER": berkas[a], "OUTPUT_NODATA_VALUE": -9999, "OUTPUT": out})
    med = gdal.Open(os.path.join(U.HASIL, "komposit_ndvi_%s_median.tif" % nama)).ReadAsArray()
    jml = gdal.Open(os.path.join(U.HASIL, "komposit_ndvi_%s_jumlah.tif" % nama)).ReadAsArray()
    print("%-12s tanggal %s sampai %s | pengamatan valid per piksel: min %d, rata-rata %.1f | piksel tanpa data: %d"
          % (nama, daftar[a][0], daftar[b - 1][0], jml.min(), jml.mean(), int((med == -9999).sum())))

Penjelasan baris demi baris.

Irisan berkas[a:b] adalah cara Python mengambil bagian daftar: dari nomor a sampai sebelum nomor b (Bab 5).

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

sebelum      tanggal 2024-01-05 sampai 2024-08-12 | pengamatan valid per piksel: min 3, rata-rata 6.4 | piksel tanpa data: 0
kemarau2024  tanggal 2024-06-13 sampai 2024-09-16 | pengamatan valid per piksel: min 2, rata-rata 3.6 | piksel tanpa data: 0
sesudah      tanggal 2025-09-21 sampai 2025-12-25 | pengamatan valid per piksel: min 1, rata-rata 2.7 | piksel tanpa data: 0

Periode "sebelum" punya rata-rata 6,4 pengamatan valid per piksel. Periode "sesudah" hanya 2,7, dan ada piksel yang hanya punya satu pengamatan. Median dari satu angka sama dengan angka itu sendiri, jadi komposit "sesudah" di piksel tersebut paling rapuh. Tidak ada piksel tanpa data (0).

Bagian yang boleh diubah.

• Baris 13, indeks periode. Tambah atau ubah periode. Saya menambahkan "kemarau2025": (17, 21). Layar menampilkan satu baris baru: kemarau2025 tanggal 2025-06-18 sampai 2025-09-21 | pengamatan valid per piksel: min 2, rata-rata 3.8 | piksel tanpa data: 0.
• Baris 15, nomor statistik. Angka adalah kode, bukan nama. Menurut bantuan alat: 0 penjumlahan (Sum), 1 banyaknya data (Count), 2 rata-rata (Mean), 3 median (Median). Saya memeriksanya lewat qgis_process help native:cellstatistics.
• Baris 18, -9999. Tanda kosong. Ubah juga pada baris 22 bila Anda mengubahnya.

Jebakan.

• Indeks akhir tidak termasuk. (5, 9) memakai tanggal ke-5, 6, 7, dan 8, bukan ke-9. Itu sebabnya tanggal akhir di baris 22 memakai b - 1.
• Nomor statistik mudah tertukar. Kode 1 bukan "satu" dalam arti terkecil, melainkan banyak pengamatan. Selalu cek hasil: peta "jumlah" harus berisi bilangan bulat kecil.
• Periode di luar daftar menghasilkan galat. Saya mengganti periode "sesudah" menjadi (24, 28). Layar menampilkan IndexError: list index out of range pada baris 18, karena berkas[24] tidak ada.
• Skrip 26.4 harus sudah jalan. Saya menjalankannya di folder kosong. Layar menampilkan Incorrect parameter value for INPUT.

Latihan variasi.

1. Tambahkan periode kemarau2025 untuk tanggal ke-17 sampai ke-20.
2. Pastikan hasilnya muncul sebagai baris keempat.

Contoh jawaban, hanya baris 13:

periode = {"sebelum": (0, 8), "kemarau2024": (5, 9), "sesudah": (20, 24), "kemarau2025": (17, 21)}

Layar menampilkan empat baris, dan baris terakhir sama seperti yang saya sebut di atas.

Skrip 26.6: Menambal lubang awan, lalu menguji tambalannya #

Berkas 2_3_isi_lubang.py. Di dalam berkasnya tertulis "SKRIP 2.3".

Tujuan. Mengisi lubang akibat awan pada deret NDVI dengan interpolasi linear sepanjang waktu, lalu menguji ketepatannya dengan menyembunyikan 20 persen pengamatan bersih.

Bahan. m3_umum.py dan paket data. Skrip ini langsung membaca citra mentah, jadi tidak memerlukan skrip lain. Jalankan dari Konsol atau terminal.

Skrip ini memuat satu fungsi buatan sendiri dan satu uji. Gambar berikut menunjukkan urutannya.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 2.3: Isi lubang akibat awan dengan interpolasi linear sepanjang waktu, lalu uji kebenarannya
# Penulis: Badar Mubarok Yogaswara
# Cara uji: sembunyikan 20% pengamatan bersih secara acak, isi dengan interpolasi, bandingkan dengan nilai aslinya.
import os
import numpy as np
import m3_umum as U

tgl, refl, scl, gt, prj = U.baca_deret()
ndvi = U.nilai_ndvi(refl)
baik = U.bersih(scl)
hari = np.array([(t - tgl[0]).days for t in tgl], dtype="float64")     # sumbu waktu dalam hari (selang tidak sama)


def isi(ndvi, baik):
    T, H, W = ndvi.shape
    hasil = np.empty_like(ndvi)
    for i in range(H):
        for j in range(W):
            ok = baik[:, i, j]
            if ok.sum() == 0:
                hasil[:, i, j] = np.nan
            else:
                # np.interp: di luar rentang, nilai ujung dipakai (tidak ekstrapolasi)
                hasil[:, i, j] = np.interp(hari, hari[ok], ndvi[ok, i, j])
    return hasil


lengkap = isi(ndvi, baik)
U.tulis_tif(os.path.join(U.HASIL, "NDVI_isi.tif"), lengkap.astype("float32"), gt, prj,
            deskripsi=[t.isoformat() for t in tgl])
print("Piksel-tanggal awan sebelum diisi: %.1f%%, sesudah: %.1f%%" % (100 * (~baik).mean(), 100 * np.isnan(lengkap).mean()))

# uji sembunyi-isi
rng = np.random.default_rng(1)
sembunyi = baik & (rng.random(baik.shape) < 0.20)
sembunyi[0], sembunyi[-1] = False, False                   # tanggal ujung tidak diuji (tanpa tetangga di satu sisi)
uji = isi(ndvi, baik & ~sembunyi)
galat = (uji - ndvi)[sembunyi]
print("Uji sembunyi-isi: %d pengamatan; RMSE %.3f; galat mutlak median %.3f" % (sembunyi.sum(), np.sqrt(np.nanmean(galat ** 2)), np.nanmedian(np.abs(galat))))
# galat per kelas acuan
dinamika, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
nama = {1: "Hutan alam", 2: "Hutan tanaman", 3: "Pertanian", 4: "Air", 5: "Terbuka", 6: "Deforestasi", 7: "Terbakar", 8: "Panen HTI"}
for k in range(1, 9):
    m = sembunyi & (dinamika == k)[None]
    e = (uji - ndvi)[m]
    print("  %-14s n=%4d  RMSE %.3f" % (nama[k], m.sum(), np.sqrt(np.nanmean(e ** 2))))

Penjelasan baris demi baris.

RMSE adalah ukuran besar galat: selisih dikuadratkan, dirata-ratakan, lalu diakarkan. Selisih besar dihukum lebih berat daripada selisih kecil. Makin kecil, makin tepat.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Piksel-tanggal awan sebelum diisi: 25.3%, sesudah: 0.0%
Uji sembunyi-isi: 5336 pengamatan; RMSE 0.078; galat mutlak median 0.014
  Hutan alam     n=2572  RMSE 0.038
  Hutan tanaman  n= 304  RMSE 0.039
  Pertanian      n= 625  RMSE 0.158
  Air            n= 277  RMSE 0.093
  Terbuka        n= 136  RMSE 0.020
  Deforestasi    n= 438  RMSE 0.117
  Terbakar       n= 520  RMSE 0.055
  Panen HTI      n= 464  RMSE 0.075

Awan sebelumnya 25,3 persen, sesudah ditambal 0,0 persen. Median galat hanya 0,014, jadi tambalan umumnya sangat dekat. RMSE total 0,078 lebih besar karena beberapa kelas sulit ditebak. Pertanian punya RMSE 0,158, empat kali hutan alam. Itu masuk akal: sawah naik turun cepat, dan garis lurus antar dua titik tidak menangkap lonjakan di antaranya.

Bagian yang boleh diubah.

• Baris 36, 0.20. Porsi yang disembunyikan. Saya menggantinya dengan 0.40: jumlah pengamatan menjadi 10700 dan RMSE naik menjadi 0.091. Makin banyak yang disembunyikan, makin sulit tebakannya.
• Baris 35, benih 1. Ganti ke 7. Jumlah pengamatan menjadi 5209 dan RMSE 0.076. Angka bergeser sedikit, tetapi gambaran besarnya sama.
• Baris 25, fungsi interpolasi. Fungsi lain bisa dipakai, asal mengembalikan nilai untuk semua hari.

Jebakan.

• Fungsi isi menduplikasi U.isi_waktu. Isi keduanya hampir sama: fungsi di kotak perkakas (m3_umum.py baris 94-103) dipakai Skrip 26.12. Saya menjalankan keduanya pada data yang sama. Hasilnya np.allclose benar, selisih terbesar 2.98e-08, dan posisi nan sama. Selisih sekecil itu hanya karena tipe angka (satu 32 bit, satu 64 bit).
• Putaran per piksel lambat untuk raster besar. Data contoh hanya 1600 piksel. Raster 10.000 kali 10.000 punya 62.500 kali lebih banyak.
• Di luar rentang, np.interp memakai nilai ujung. Ia tidak menarik garis ke luar (tidak ekstrapolasi). Karena itu tanggal ujung dikeluarkan dari uji (baris 37).
• Ukuran dan luas tidak dikunci di skrip ini. T, H, W diambil dari data. Skrip lain di bab ini mengunci 40 kali 40 piksel dan 0,01 hektare per piksel, dan itu dicatat di tempatnya.

Latihan variasi.

1. Sembunyikan 40 persen pengamatan bersih.
2. Bandingkan RMSE totalnya dengan hasil asli (0,078).

Contoh jawaban, hanya baris 36:

sembunyi = baik & (rng.random(baik.shape) < 0.40)

Layar menampilkan Uji sembunyi-isi: 10700 pengamatan; RMSE 0.091; galat mutlak median 0.015. RMSE naik dari 0,078 menjadi 0,091.

Kelompok C: Mengukur irama tiap tutupan lahan #

Konsep: irama dalam tiga kalimat #

Tiap tutupan lahan punya sidik jari musiman, seperti detak jantung tiap orang berbeda: hutan berdetak tenang, sawah berdetak cepat dan tinggi. Anda bisa meringkas detak itu dalam beberapa angka per piksel: rata-rata, seberapa tinggi naik turunnya, ke mana arah umumnya, dan kapan puncaknya. Angka ringkas ini disebut metrik deret waktu.

Lima istilah yang perlu Anda kenal sekarang:

• Amplitudo: selisih nilai tertinggi dan terendah dalam satu deret.
• Simpangan baku: ukuran seberapa jauh nilai biasanya menyimpang dari rata-ratanya.
• Tren: kemiringan garis lurus yang paling cocok dengan deret. Di bab ini dinyatakan sebagai perubahan NDVI per tahun.
• Hari puncak: hari ke-berapa (sejak tanggal pertama) NDVI mencapai maksimum.
• Anomali: selisih nilai satu tanggal dari rata-rata piksel itu sendiri.

Skrip 26.7: Metrik deret waktu per piksel #

Berkas 3_1_metrik_deret.py. Di dalam berkasnya tertulis "SKRIP 3.1".

Tujuan. Menghitung tujuh metrik NDVI tiap piksel, meringkasnya per kelas tutupan lahan, dan menghitung anomali tiap tanggal.

Bahan. Hasil Skrip 26.6 (NDVI_isi.tif), peta acuan acuan/Acuan_Dinamika.tif (delapan kelas tutupan), dan m3_umum.py. Skrip ini hanya memakai NumPy dan GDAL, jadi bisa dijalankan dari Konsol maupun terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 3.1: Metrik deret waktu per piksel: rata-rata, minimum, maksimum, amplitudo, simpangan baku, tren per tahun, tanggal puncak
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 2.3 sudah dijalankan (hasil/NDVI_isi.tif). Menghasilkan hasil/metrik_ndvi.tif (7 band).
import os
import numpy as np
import m3_umum as U

tgl = [t for t, _, _ in U.daftar_citra()]
hari = np.array([(t - tgl[0]).days for t in tgl], dtype="float64")
x, gt, prj = U.baca_tif(os.path.join(U.HASIL, "NDVI_isi.tif"))          # (24, H, W)
rata, mini, maks = x.mean(0), x.min(0), x.max(0)
amplitudo = maks - mini
sb = x.std(0)
# tren: kemiringan garis lurus terhadap waktu dalam HARI, lalu dikali 365 -> perubahan NDVI per tahun
hc = hari - hari.mean()
kemiringan = (hc[:, None, None] * (x - rata)).sum(0) / (hc ** 2).sum()
tren_tahun = kemiringan * 365.0
puncak = hari[x.argmax(0)]                                              # hari ke- (sejak tanggal pertama) saat NDVI tertinggi
nama = ["rata", "minimum", "maksimum", "amplitudo", "simpangan_baku", "tren_per_tahun", "hari_puncak"]
U.tulis_tif(os.path.join(U.HASIL, "metrik_ndvi.tif"), np.array([rata, mini, maks, amplitudo, sb, tren_tahun, puncak]).astype("float32"), gt, prj, deskripsi=nama)

din, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
kls = {1: "Hutan alam", 2: "Hutan tanaman", 3: "Pertanian", 4: "Air", 5: "Terbuka", 6: "Deforestasi", 7: "Terbakar", 8: "Panen HTI"}
print("%-14s %6s %9s %9s %12s" % ("Kelas", "rata", "amplitudo", "simp.baku", "tren/tahun"))
for k, n in kls.items():
    m = din == k
    print("%-14s %6.2f %9.2f %9.2f %12.3f" % (n, rata[m].mean(), amplitudo[m].mean(), sb[m].mean(), tren_tahun[m].mean()))

# anomali: selisih NDVI tiap tanggal terhadap rata-rata piksel itu
anomali = x - rata[None]
U.tulis_tif(os.path.join(U.HASIL, "anomali_ndvi.tif"), anomali.astype("float32"), gt, prj, deskripsi=[t.isoformat() for t in tgl])
i_sep = tgl.index(__import__("datetime").date(2024, 9, 16))
m = din == 7
print("Anomali NDVI rata-rata di bekas terbakar pada 2024-09-16: %.3f; di hutan stabil: %.3f" % (anomali[i_sep][m].mean(), anomali[i_sep][din == 1].mean()))

Penjelasan baris demi baris.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Kelas            rata amplitudo simp.baku   tren/tahun
Hutan alam       0.86      0.10      0.02       -0.004
Hutan tanaman    0.77      0.10      0.02       -0.006
Pertanian        0.50      0.46      0.15       -0.027
Air             -0.33      0.29      0.07        0.003
Terbuka          0.11      0.06      0.02       -0.001
Deforestasi      0.64      0.59      0.22       -0.281
Terbakar         0.82      0.23      0.06       -0.022
Panen HTI        0.61      0.51      0.17       -0.200
Anomali NDVI rata-rata di bekas terbakar pada 2024-09-16: -0.115; di hutan stabil: -0.015

Bacalah tabelnya per kelas. Hutan alam stabil: rata-rata 0,86, amplitudo hanya 0,10, hampir tidak punya tren. Pertanian naik turun: amplitudo 0,46. Deforestasi dan panen hutan tanaman punya tren turun tajam (-0,281 dan -0,200 per tahun). Baris terakhir menunjukkan tanda awal kebakaran: pada 16 September 2024, NDVI di bekas terbakar turun 0,115 dari kebiasaannya, sedangkan hutan stabil hanya 0,015.

Bagian yang boleh diubah.

• Baris 18, 365.0. Pengali satuan tren. Saya mengubahnya menjadi 30.0, yaitu perubahan per 30 hari. Tren deforestasi di tabel berubah dari -0.281 menjadi -0.023.
• Baris 33, tanggal anomali. Ganti ke tanggal lain yang ada di Daftar_Citra.csv. Saya mengganti menjadi date(2025, 8, 17). Layar menampilkan Anomali NDVI rata-rata di bekas terbakar pada 2024-09-16: -0.031; di hutan stabil: -0.020.
• Baris 12-14, statistik tambahan. Tambah x.std(0) atau np.median(x, axis=0) untuk metrik baru, lalu tambahkan namanya ke daftar baris 20.

Jebakan.

• Teks di print tidak ikut berubah. Pada contoh di atas, tanggal diganti tetapi teksnya masih tertulis 2024-09-16. Angka yang tampil memang milik 17 Agustus 2025, tetapi labelnya salah. Setelah mengganti tanggal, ganti juga tulisannya di baris 35. Pada tanggal itu pun angkanya mencampur dua kebakaran, sebab kelas 7 mencakup bekas bakar tahun 2024 dan 2025.
• Tanggal harus ada di daftar. Saya mengganti baris 33 menjadi date(2024, 9, 17) dan layar menampilkan ValueError: datetime.date(2024, 9, 17) is not in list.
• Skrip 26.6 harus sudah jalan. Di folder kosong, layar menampilkan RuntimeError: ...NDVI_isi.tif: No such file or directory.
• Selang antar tanggal tidak sama. Karena itu tren dihitung terhadap hari, bukan nomor urut. Tren pada data dua tahun hanya cukup untuk belajar, belum untuk laporan resmi [CEK].

Latihan variasi.

1. Cetak juga hari puncak rata-rata tiap kelas.
2. Bandingkan kelas yang kemudian berubah (deforestasi, terbakar) dengan hutan alam.

Contoh jawaban, tambahkan di akhir skrip:

for k, n in kls.items():
    print("%-14s hari puncak rata-rata: %5.0f" % (n, puncak[din == k].mean()))

Saya menjalankannya. Layar menampilkan delapan baris tambahan. Tiga di antaranya: Hutan alam hari puncak rata-rata: 303, Deforestasi hari puncak rata-rata: 113, dan Terbakar hari puncak rata-rata: 87. Kelas yang kemudian berubah mencapai puncaknya jauh lebih awal (sekitar hari ke-100) daripada hutan alam (hari ke-303). Wajar: setelah kejadian, NDVI-nya jatuh dan tidak naik lagi.

Kelompok D: Mendeteksi hutan yang hilang #

Konsep: dua foto lawan seluruh film #

Deteksi dua waktu itu seperti permainan "temukan perbedaan" antara foto sebelum dan sesudah: cepat dan mudah, tetapi hanya tahu bahwa sesuatu berubah, bukan kapan. Deteksi titik patah menonton seluruh film setiap piksel, dan mencatat saat nilainya jatuh lalu tidak kembali. Cara kedua bisa menyebut tanggal, tetapi lebih lambat.

Lima istilah yang perlu Anda kenal sekarang:

• dNDVI: NDVI sesudah dikurangi NDVI sebelum. Nilai sangat negatif berarti kehijauan hilang.
• Ambang: batas angka untuk menyebut sesuatu "berubah".
• Luas minimum: gumpalan piksel yang lebih kecil dari batas ini dibuang sebagai derau.
• Titik patah: saat sebuah deret berubah nilai secara menetap.
• Pengamatan berturut-turut: syarat bahwa penurunan harus bertahan beberapa tanggal, supaya awan yang lolos tidak menipu.

Kepala Seksi berpesan satu hal: panen hutan tanaman bukan deforestasi. Itu tebangan terencana. Kedua skrip di bawah membuang kawasan hutan tanaman memakai peta Kawasan_HTI.gpkg.

Skrip 26.8: Deforestasi dari dua waktu #

Berkas 4_1_deforestasi_dua_waktu.py. Di dalam berkasnya tertulis "SKRIP 4.1".

Tujuan. Mendeteksi deforestasi dari selisih komposit NDVI sebelum dan sesudah, dengan ambang dan saringan luas minimum. Lalu membuang kawasan hutan tanaman, dan memeriksa hasilnya terhadap peta acuan.

Bahan. Hasil Skrip 26.5 (komposit median sebelum dan sesudah), vektor/Kawasan_HTI.gpkg, acuan/Acuan_Perubahan.tif, acuan/Acuan_Dinamika.tif, dan m3_umum.py. Jalankan dari Konsol atau terminal.

Skrip ini panjang dan memuat sembilan tahap. Gambar berikut membantu Anda membaginya menjadi bagian kecil.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 4.1: Deteksi deforestasi dua waktu: selisih komposit NDVI (sesudah dikurangi sebelum), ambang, saringan luas minimum
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 2.2 sudah dijalankan. Ambang dan luas minimum bisa Anda ubah di bagian PARAMETER.
import os
import numpy as np
import processing
from qgis.core import QgsRasterLayer, QgsVectorLayer, QgsProject
import m3_umum as U

AMBANG = -0.25          # selisih NDVI di bawah nilai ini dianggap kehilangan tutupan
LUAS_MIN_PIKSEL = 5     # gumpalan lebih kecil dari ini (5 piksel = 0,05 ha) dibuang sebagai derau

sebelum = QgsRasterLayer(os.path.join(U.HASIL, "komposit_ndvi_sebelum_median.tif"), "sebelum")
sesudah = QgsRasterLayer(os.path.join(U.HASIL, "komposit_ndvi_sesudah_median.tif"), "sesudah")
QgsProject.instance().addMapLayers([sebelum, sesudah], False)
d_ndvi = os.path.join(U.HASIL, "dNDVI.tif")
processing.run("native:rastercalc", {"LAYERS": [sebelum, sesudah], "EXPRESSION": '"sesudah@1" - "sebelum@1"', "OUTPUT": d_ndvi})
# 1) ambang: hanya piksel yang hutan sebelumnya (NDVI sebelum >= 0.7) dan turun lewat AMBANG
kand = os.path.join(U.HASIL, "kandidat_deforestasi.tif")
processing.run("native:rastercalc", {"LAYERS": [sebelum, sesudah],
    "EXPRESSION": '(("sebelum@1" >= 0.7) AND (("sesudah@1" - "sebelum@1") < %s)) * 1' % AMBANG, "OUTPUT": kand})
# 2) saring gumpalan kecil
saring = os.path.join(U.HASIL, "kandidat_disaring.tif")
processing.run("gdal:sieve", {"INPUT": kand, "THRESHOLD": LUAS_MIN_PIKSEL, "EIGHT_CONNECTEDNESS": True, "NO_MASK": True, "OUTPUT": saring})
# 3) buang kawasan hutan tanaman (panen HTI itu terencana, bukan deforestasi)
hti = os.path.join(U.HASIL, "mask_hti.tif")
processing.run("gdal:rasterize", {"INPUT": os.path.join(U.PAKET, "vektor", "Kawasan_HTI.gpkg"), "BURN": 1, "UNITS": 1,
    "WIDTH": 10, "HEIGHT": 10, "EXTENT": "312000,312400,9996000,9996400 [EPSG:32749]", "INIT": 0, "NODATA": 255, "DATA_TYPE": 0, "OUTPUT": hti})
final = os.path.join(U.HASIL, "deforestasi_2waktu.tif")
l1, l2 = QgsRasterLayer(saring, "saring"), QgsRasterLayer(hti, "hti")
QgsProject.instance().addMapLayers([l1, l2], False)
processing.run("native:rastercalc", {"LAYERS": [l1, l2], "EXPRESSION": '("saring@1" = 1) AND ("hti@1" = 0)', "OUTPUT": final})

# pemeriksaan terhadap peta acuan
hasil, gt, _ = U.baca_tif(final)
hasil = hasil == 1
acuan, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Perubahan.tif"))
benar = acuan == 1
tp = (hasil & benar).sum(); fp = (hasil & ~benar).sum(); fn = (~hasil & benar).sum()
print("Terdeteksi: %d piksel = %.2f ha; acuan: %d piksel = %.2f ha" % (hasil.sum(), hasil.sum() * 0.01, benar.sum(), benar.sum() * 0.01))
print("Tepat %d, salah tambah %d, terlewat %d; presisi %.2f, kepekaan %.2f" % (tp, fp, fn, tp / max(tp + fp, 1), tp / max(tp + fn, 1)))
din, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
nama = {1: "Hutan alam", 2: "Hutan tanaman", 3: "Pertanian", 4: "Air", 5: "Terbuka", 6: "Deforestasi", 7: "Terbakar", 8: "Panen HTI"}
print("Piksel terdeteksi menurut kelas acuan:", {nama[k]: int((hasil & (din == k)).sum()) for k in range(1, 9) if (hasil & (din == k)).sum()})
rr, cc = np.mgrid[0:40, 0:40]
kph = np.where(rr < 20, "Gamma", np.where(cc < 20, "Alpha", "Beta"))
for k in ("Alpha", "Beta", "Gamma"):
    print("  %-6s deforestasi terdeteksi %.2f ha, acuan %.2f ha" % (k, (hasil & (kph == k)).sum() * 0.01, (benar & (kph == k)).sum() * 0.01))
poli = os.path.join(U.HASIL, "deforestasi_2waktu.gpkg")
if os.path.exists(poli):
    os.remove(poli)                  # polygonize ke GPKG yang sudah ada MENAMBAH fitur, bukan menimpa
processing.run("gdal:polygonize", {"INPUT": final, "BAND": 1, "FIELD": "DN", "EIGHT_CONNECTEDNESS": True, "OUTPUT": poli})
print("Poligon tersimpan:", poli)



# uji kepekaan ambang (tanpa saringan luas minimum, dengan dan tanpa peta HTI)
dn, _, _ = U.baca_tif(os.path.join(U.HASIL, "dNDVI.tif"))
sb, _, _ = U.baca_tif(os.path.join(U.HASIL, "komposit_ndvi_sebelum_median.tif"))
htim, _, _ = U.baca_tif(hti)
print("Ambang   dengan peta HTI (tepat/salah/lewat)   tanpa peta HTI (tepat/salah/lewat)")
for a in (-0.15, -0.25, -0.35, -0.45):
    k = (sb >= 0.7) & (dn < a)
    for_hti = k & (htim == 0)
    t1, f1, l1_ = (for_hti & benar).sum(), (for_hti & ~benar).sum(), (~for_hti & benar).sum()
    t2, f2, l2_ = (k & benar).sum(), (k & ~benar).sum(), (~k & benar).sum()
    print("%6.2f   %4d / %4d / %4d                    %4d / %4d / %4d" % (a, t1, f1, l1_, t2, f2, l2_))

Penjelasan baris demi baris. Baris 55-57 kosong, jadi tidak muncul di tabel.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Terdeteksi: 130 piksel = 1.30 ha; acuan: 130 piksel = 1.30 ha
Tepat 130, salah tambah 0, terlewat 0; presisi 1.00, kepekaan 1.00
Piksel terdeteksi menurut kelas acuan: {'Deforestasi': 130}
  Alpha  deforestasi terdeteksi 0.30 ha, acuan 0.30 ha
  Beta   deforestasi terdeteksi 0.00 ha, acuan 0.00 ha
  Gamma  deforestasi terdeteksi 1.00 ha, acuan 1.00 ha
Poligon tersimpan: D:/Latihan_M3/hasil\deforestasi_2waktu.gpkg
Ambang   dengan peta HTI (tepat/salah/lewat)   tanpa peta HTI (tepat/salah/lewat)
 -0.15    130 /   13 /    0                     130 /  133 /    0
 -0.25    130 /    1 /    0                     130 /   24 /    0
 -0.35    118 /    1 /   12                     118 /    1 /   12
 -0.45     26 /    0 /  104                      26 /    0 /  104

Hasilnya sempurna: 130 piksel (1,30 hektare) terdeteksi, semuanya benar, tidak ada yang terlewat. Rincian per KPH cocok dengan acuan. Data ini sintetis, jadi jangan berharap sebersih ini pada data nyata. Tabel uji ambang adalah bagian paling berguna. Pada ambang -0,15, kedua versi sama-sama menangkap semua 130 piksel. Namun versi dengan peta hutan tanaman salah 13 piksel, sedangkan versi tanpa peta salah 133 piksel. Itulah gunanya membuang kawasan hutan tanaman.

Bagian yang boleh diubah.

• Baris 11, AMBANG. Makin negatif, makin ketat. Saya mencoba -0.35: layar menampilkan Terdeteksi: 118 piksel = 1.18 ha dan Tepat 118, salah tambah 0, terlewat 12; presisi 1.00, kepekaan 0.91.
• Baris 12, LUAS_MIN_PIKSEL. Ubah menjadi 1: terdeteksi 131 piksel, dengan satu salah tambah (di kelas Pertanian). Ubah menjadi 60: terdeteksi hanya 100 piksel, dan bukaan kecil di Alpha yang 30 piksel ikut terbuang (terlewat 30, kepekaan 0,77). Saringan yang terlalu besar membuang kejadian nyata.
• Baris 22, 0.7. Ambang NDVI hutan sebelum perubahan.

Jebakan.

• Banyak hal dikunci untuk data contoh. EXTENT di baris 29 (312000,312400,9996000,9996400), ukuran piksel 10 meter, ukuran 40 kali 40 di baris 46, luas 0,01 hektare per piksel di baris 41 dan 49, serta pembagian KPH rr < 20 di baris 47. Pada data lain, semuanya harus diganti, atau hasilnya salah tanpa pesan galat.
• Alat Polygonize menambah, bukan menimpa. Itu sebabnya baris 51-52 menghapus deforestasi_2waktu.gpkg lebih dulu. Saya menjalankan skrip dua kali pada folder yang sama, dan jumlah poligon tetap 3 pada jalan kedua. Tanpa dua baris itu, jumlahnya naik menjadi 6. Bila Windows menolak penghapusan karena berkas itu masih terbuka sebagai layer di QGIS, os.remove gagal dengan PermissionError: [WinError 32] (diuji; layer di kanvas jendela QGIS tidak diuji). Hapus layernya dari proyek, lalu jalankan ulang.
• Nama layer harus sama persis dengan nama di rumus. "sebelum@1" mengacu pada layer bernama sebelum.
• Skrip 26.5 harus sudah jalan. Di folder kosong, layar menampilkan QgsProcessingException: Error creating output file. Pesan itu menyesatkan: yang hilang adalah masukan, bukan keluaran.
• Skrip ini panjang dan mencampur tiga pekerjaan: mendeteksi, memeriksa, dan menguji ambang. Memecahnya menjadi tiga berkas akan lebih rapi, tetapi buku ini membahasnya apa adanya.

Latihan variasi.

1. Ubah ambang menjadi -0.35, lalu bandingkan hasilnya dengan ambang asli.
2. Ubah luas minimum menjadi 1, lalu cari tahu piksel salah tambah itu masuk kelas apa.

Contoh jawaban, hanya baris yang berubah:

AMBANG = -0.35          # selisih NDVI di bawah nilai ini dianggap kehilangan tutupan
LUAS_MIN_PIKSEL = 1     # gumpalan lebih kecil dari ini dibuang sebagai derau

Untuk ambang -0.35, layar menampilkan 118 piksel (turun dari 130), dengan 12 terlewat. Untuk luas minimum 1, layar menampilkan Piksel terdeteksi menurut kelas acuan: {'Pertanian': 1, 'Deforestasi': 130}, jadi piksel salah itu ada di kelas Pertanian. Tiap variasi saya jalankan terpisah, dengan salah satu konstanta saja yang berubah.

Skrip 26.9: Titik patah, kapan hutan mulai hilang #

Berkas 4_2_titik_patah.py. Di dalam berkasnya tertulis "SKRIP 4.2".

Tujuan. Mencari kapan tiap piksel hutan mulai hilang: NDVI turun jauh dari dasarnya pada beberapa pengamatan bersih berturut-turut. Lalu menyimpan tanggalnya dan membandingkannya dengan tanggal sebenarnya.

Bahan. Citra mentah (lewat m3_umum.py), peta acuan (Acuan_Dinamika.tif dan Acuan_Perubahan.tif), dan mask_hti.tif dari Skrip 26.8. Skrip ini tidak memakai QGIS. Jalankan dari Konsol atau terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 4.2: Deteksi titik patah pada deret waktu NDVI: kapan hutan mulai hilang?
# Penulis: Badar Mubarok Yogaswara
# Aturan: dasar = median 6 pengamatan bersih pertama. Perubahan = NDVI turun lebih dari AMBANG dari dasar
# pada K pengamatan bersih berturut-turut. Tanggal perubahan = pengamatan pertama dalam rangkaian itu.
import datetime as dt
import os
import numpy as np
from osgeo import gdal
import m3_umum as U

AMBANG, K, N_DASAR = 0.25, 3, 6
tgl, refl, scl, gt, prj = U.baca_deret()
ndvi = U.nilai_ndvi(refl)
bersih = U.bersih(scl)
T, H, W = ndvi.shape
hari = np.array([(t - tgl[0]).days for t in tgl])


def deteksi(pakai_mask, ambang=AMBANG, k=K):
    ok = bersih if pakai_mask else np.ones_like(bersih)
    hasil = np.full((H, W), -1)            # indeks tanggal awal perubahan, -1 = tidak ada
    for i in range(H):
        for j in range(W):
            idx = np.where(ok[:, i, j])[0]
            if len(idx) < N_DASAR + k:
                continue
            dasar = np.median(ndvi[idx[:N_DASAR], i, j])
            if dasar < 0.7:                # hanya piksel yang awalnya berhutan
                continue
            deret = ndvi[idx[N_DASAR:], i, j]
            turun = deret < dasar - ambang
            for n in range(len(turun) - k + 1):
                if turun[n:n + k].all():
                    hasil[i, j] = idx[N_DASAR + n]
                    break
    return hasil


din, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
per, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Perubahan.tif"))
hti, _, _ = U.baca_tif(os.path.join(U.HASIL, "mask_hti.tif"))     # dibuat Skrip 4.1
benar = per == 1
print("K  mask-awan  tepat salah(di luar HTI) lewat")
for pakai in (False, True):
    for k in (1, 2, 3):
        h = deteksi(pakai, k=k) >= 0
        h = h & (hti == 0)
        print("%d  %-9s  %5d %5d %5d" % (k, "ya" if pakai else "tidak", (h & benar).sum(), (h & ~benar).sum(), (~h & benar).sum()))

idx = deteksi(True)
flag = (idx >= 0) & (hti == 0)
keluar = np.where(flag, np.array([int(tgl[i].strftime("%Y%m%d")) if i >= 0 else 0 for i in idx.ravel()]).reshape(H, W), 0).astype("int32")
U.tulis_tif(os.path.join(U.HASIL, "tanggal_deforestasi.tif"), keluar, gt, prj, nodata=0, tipe=gdal.GDT_Int32)
# tanggal sebenarnya (acuan): bukaan 1 pada 2024-10-28, bukaan 2 pada 2025-06-05
benar_tgl = {"Bukaan 1 (Gamma)": (dt.date(2024, 10, 28), benar & (np.mgrid[0:H, 0:W][0] < 20)),
             "Bukaan 2 (Alpha)": (dt.date(2025, 6, 5), benar & (np.mgrid[0:H, 0:W][0] >= 20))}
for nama, (tg, m) in benar_tgl.items():
    terdeteksi = [tgl[idx[i, j]] for i, j in zip(*np.where(m & flag))]
    tunda = [(d - tg).days for d in terdeteksi]
    pertama = min(terdeteksi)
    print("%s: tanggal sebenarnya %s | terdeteksi %d piksel | tanggal terdeteksi paling awal %s, paling akhir %s | tunda rata-rata %.0f hari"
          % (nama, tg, len(terdeteksi), pertama, max(terdeteksi), np.mean(tunda)))

Penjelasan baris demi baris. Fungsi deteksi adalah inti skrip ini, jadi diuraikan lebih rinci.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

K  mask-awan  tepat salah(di luar HTI) lewat
1  tidak        128   873     2
2  tidak        128   593     2
3  tidak        128   156     2
1  ya           130    92     0
2  ya           130     8     0
3  ya           130     0     0
Bukaan 1 (Gamma): tanggal sebenarnya 2024-10-28 | terdeteksi 100 piksel | tanggal terdeteksi paling awal 2024-11-15, paling akhir 2025-01-14 | tunda rata-rata 24 hari
Bukaan 2 (Alpha): tanggal sebenarnya 2025-06-05 | terdeteksi 30 piksel | tanggal terdeteksi paling awal 2025-06-18, paling akhir 2025-06-18 | tunda rata-rata 13 hari

Bacalah tabel enam baris itu. Tanpa mask awan dan dengan K = 1, ada 873 piksel salah: satu awan saja sudah dianggap hutan hilang. Dengan mask awan dan K = 3, salahnya nol dan semua 130 piksel benar. Dua syarat itu, yaitu membuang awan dan menunggu beberapa pengamatan berturut-turut, mengubah hasil total. Tundaan rata-rata 24 hari untuk bukaan Gamma kemungkinan berasal dari selang antar citra yang sekitar 30 hari dan dari tiga pengamatan yang harus ditunggu [kemungkinan].

Bagian yang boleh diubah.

• Baris 12, K. Ubah ke 2. Hasil akhir bertambah dari 130 menjadi 138 piksel bertanggal, karena delapan piksel salah ikut masuk (sama seperti baris 2 ya di tabel uji). Saya menambahkan print("Piksel bertanggal:", int(flag.sum())) di akhir untuk melihat angkanya.
• Baris 12, AMBANG. Makin kecil, makin peka tetapi makin mudah tertipu.
• Baris 29, 0.7. Batas NDVI hutan untuk dasar.

Jebakan.

Tipe data harus disebut sendiri (baris 54). Fungsi U.tulis_tif memilih Byte untuk semua larik non-desimal (m3_umum.py baris 65-66). Byte hanya memuat 0 sampai 255, jadi bilangan seperti 20241115 akan terpotong menjadi 255. Karena itu baris 54 menyebut tipe=gdal.GDT_Int32. Saya memeriksa berkas hasil dengan potongan ini:

from osgeo import gdal
import numpy as np
a = gdal.Open(r"D:/Latihan_M3/hasil/tanggal_deforestasi.tif")
print(gdal.GetDataTypeName(a.GetRasterBand(1).DataType), np.unique(a.ReadAsArray()))

Int32 [       0 20241115 20241220 20250114 20250618]

Lima nilai itu benar: nol untuk tidak ada perubahan, ditambah empat tanggal. Tanpa tipe=, berkas yang sama terbaca Byte [ 0 255]: semua tanggal rusak menjadi 255, dan skrip tidak mengeluh. Ringkasan di layar (baris 58-63) tidak terpengaruh karena memakai larik di memori, bukan berkas. Hanya berkas .tif yang salah. Jebakan lainnya:

• min() pada daftar kosong menimbulkan galat. Saya mengganti AMBANG menjadi 0.99, sehingga tidak ada piksel terdeteksi. Layar menampilkan ValueError: min() iterable argument is empty di baris 61.
• Skrip 26.8 harus sudah jalan. Di folder kosong, layar menampilkan RuntimeError: ...mask_hti.tif: No such file or directory.
• Tanggal sebenarnya dan pembagian wilayah dikunci di kode. Baris 56-57 hanya cocok dengan data contoh.
• Putaran per piksel lambat untuk raster besar, seperti di Skrip 26.6.

Latihan variasi.

1. Ubah K menjadi 2 dan lihat berapa piksel bertanggal.
2. Hapus tipe=gdal.GDT_Int32 dari baris 54, jalankan ulang, lalu periksa tipe berkasnya.

Contoh jawaban nomor 2, satu perubahan pada baris 54. Buang tipe:

U.tulis_tif(os.path.join(U.HASIL, "tanggal_deforestasi.tif"), keluar, gt, prj, nodata=0)

Layar tidak menampilkan galat. Memeriksa berkas dengan potongan di Jebakan, layar menampilkan Byte [ 0 255]: semua tanggal terpotong menjadi 255. Kembalikan tipe=gdal.GDT_Int32 setelah mencoba. Angka K=2 untuk nomor 1 adalah 138.

Kelompok E: Jejak api #

Konsep: arang yang menyerap hijau #

Tumbuhan sehat memantulkan banyak inframerah dekat dan sedikit inframerah gelombang pendek. Setelah terbakar, keadaannya terbalik, seperti arang yang menelan cahaya. Indeks NBR menangkap pembalikan itu. Selisih NBR antara sebelum dan sesudah kebakaran, yaitu dNBR, menunjukkan seberapa parah terbakarnya.

Lima istilah yang perlu Anda kenal sekarang:

• SWIR: inframerah gelombang pendek. Di Sentinel-2, SWIR2 adalah band B12.
• NBR: (inframerah dekat dikurangi SWIR2) dibagi (inframerah dekat ditambah SWIR2).
• dNBR: NBR sebelum dikurangi NBR sesudah. Makin besar, makin parah.
• Titik panas (hotspot): titik deteksi suhu tinggi oleh satelit. Titik panas adalah petunjuk api, bukan peta bekas bakar.
• Penyangga (buffer): daerah selebar jarak tertentu di sekeliling sebuah objek.

Mengapa NBR, bukan NDVI saja? Pada contoh grafik di atas, NDVI hanya turun sekitar 0,3, sedangkan NBR turun hampir 0,6. Gerakan yang lebih besar lebih mudah dibedakan dari derau.

Skrip 26.10: Bekas kebakaran dengan dNBR #

Berkas 5_1_dnbr.py. Di dalam berkasnya tertulis "SKRIP 5.1".

Tujuan. Memetakan bekas kebakaran dua kejadian (September 2024 di KPH Beta, Agustus 2025 di KPH Gamma) dengan NBR dan dNBR, menguji beberapa ambang, dan menyimpan peta akhir.

Bahan. Citra mentah (berkas _refl.tif dan _SCL.tif), acuan/Acuan_Dinamika.tif, acuan/Tutupan_Awal.tif, dan m3_umum.py. Skrip ini tidak memerlukan skrip lain. Jalankan dari Konsol atau terminal.

Skrip ini panjang. Gambar berikut membaginya menjadi delapan tahap.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 5.1: Bekas kebakaran dengan NBR dan dNBR (dua kejadian: September 2024 dan Agustus 2025)
# Penulis: Badar Mubarok Yogaswara
# NBR = (NIR - SWIR2) / (NIR + SWIR2); untuk Sentinel-2: NIR = B08 (band 4 di berkas), SWIR2 = B12 (band 6 di berkas).
# dNBR = NBR sebelum - NBR sesudah. Makin besar, makin parah terbakar.
import os
import numpy as np
import processing
from qgis.core import QgsRasterLayer, QgsProject
import m3_umum as U

AMBANG_BAKAR = 0.10          # dNBR di atas ini dianggap bekas terbakar (lihat uji ambang di akhir keluaran)
daftar = U.daftar_citra()
os.makedirs(os.path.join(U.HASIL, "nbr_bersih"), exist_ok=True)
nbr = []
for tgl, pr, ps in daftar:
    k = tgl.strftime("%Y%m%d")
    a, s = QgsRasterLayer(pr, "R" + k), QgsRasterLayer(ps, "S" + k)
    QgsProject.instance().addMapLayers([a, s], False)
    ek = '(("R%s@4" - "R%s@6") / ("R%s@4" + "R%s@6")) / (("S%s@1" = 4) OR ("S%s@1" = 5) OR ("S%s@1" = 6))' % ((k,) * 7)
    out = os.path.join(U.HASIL, "nbr_bersih", "NBR_bersih_%s.tif" % k)
    processing.run("native:rastercalc", {"LAYERS": [a, s], "EXPRESSION": ek, "OUTPUT": out})
    nbr.append(out)


def median(idx, nama):
    out = os.path.join(U.HASIL, nama)
    processing.run("native:cellstatistics", {"INPUT": [nbr[i] for i in idx], "STATISTIC": 3, "IGNORE_NODATA": True,
                                              "REFERENCE_LAYER": nbr[idx[0]], "OUTPUT_NODATA_VALUE": -9999, "OUTPUT": out})
    return out


tg = [t for t, _, _ in daftar]
kejadian = {"2024": ([5, 6, 7], [8, 9], "Beta"), "2025": ([17, 18], [19, 20], "Gamma")}   # indeks tanggal sebelum, sesudah, KPH
din, gt, prj = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
awal, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Tutupan_Awal.tif"))
rr, cc = np.mgrid[0:40, 0:40]
kph = np.where(rr < 20, "Gamma", np.where(cc < 20, "Alpha", "Beta"))
sapuan = {}
for nama, (pre, post, wilayah) in kejadian.items():
    a, b = median(pre, "nbr_pre_%s.tif" % nama), median(post, "nbr_post_%s.tif" % nama)
    la, lb = QgsRasterLayer(a, "pre"), QgsRasterLayer(b, "post")
    QgsProject.instance().addMapLayers([la, lb], False)
    d = os.path.join(U.HASIL, "dNBR_%s.tif" % nama)
    processing.run("native:rastercalc", {"LAYERS": [la, lb], "EXPRESSION": '"pre@1" - "post@1"', "OUTPUT": d})
    dn, _, _ = U.baca_tif(d)
    pre_n, _, _ = U.baca_tif(a)
    dn = np.where(dn < -100, np.nan, dn)                      # -9999 = tidak ada data
    acuan = (din == 7) & (kph == wilayah)                     # bekas terbakar kejadian ini (KPH tempat kejadian)
    # hanya hutan dan hutan tanaman menurut peta tutupan awal, yang NBR-nya tinggi (air dan sawah dibuang: NBR-nya tidak stabil)
    layak = (pre_n >= 0.5) & np.isin(awal, (1, 2))
    print("Kejadian %s (KPH %s): sebelum %s s.d. %s, sesudah %s s.d. %s" % (nama, wilayah, tg[pre[0]], tg[pre[-1]], tg[post[0]], tg[post[-1]]))
    print("  dNBR rata-rata di bekas terbakar acuan %.2f; di hutan stabil %.2f; acuan %d piksel (%.2f ha)"
          % (np.nanmean(dn[acuan]), np.nanmean(dn[din == 1]), acuan.sum(), acuan.sum() * 0.01))
    print("  ambang   terdeteksi(ha)  tepat salah lewat  presisi kepekaan")
    for amb in (0.05, 0.10, 0.20, 0.27, 0.44):
        bk = layak & (dn > amb)
        tp, fp, fn = (bk & acuan).sum(), (bk & ~acuan).sum(), (~bk & acuan).sum()
        print("  %5.2f   %8.2f      %4d %5d %4d    %.2f    %.2f" % (amb, bk.sum() * 0.01, tp, fp, fn, tp / max(tp + fp, 1), tp / max(tp + fn, 1)))
    sapuan[nama] = layak & (dn > AMBANG_BAKAR)
akhir = sapuan["2024"] | sapuan["2025"]
U.tulis_tif(os.path.join(U.HASIL, "bekas_bakar.tif"), akhir.astype("uint8"), gt, prj, nodata=255)
for k in ("Alpha", "Beta", "Gamma"):
    print("KPH %-6s terbakar terdeteksi %.2f ha (acuan %.2f ha)" % (k, (akhir & (kph == k)).sum() * 0.01, ((din == 7) & (kph == k)).sum() * 0.01))

Penjelasan baris demi baris.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Kejadian 2024 (KPH Beta): sebelum 2024-06-13 s.d. 2024-08-12, sesudah 2024-09-16 s.d. 2024-10-16
  dNBR rata-rata di bekas terbakar acuan 0.35; di hutan stabil -0.01; acuan 128 piksel (1.28 ha)
  ambang   terdeteksi(ha)  tepat salah lewat  presisi kepekaan
   0.05       1.38       128    10    0    0.93    1.00
   0.10       1.29       128     1    0    0.99    1.00
   0.20       1.14       114     0   14    1.00    0.89
   0.27       0.91        91     0   37    1.00    0.71
   0.44       0.30        30     0   98    1.00    0.23
Kejadian 2025 (KPH Gamma): sebelum 2025-06-18 s.d. 2025-07-18, sesudah 2025-08-17 s.d. 2025-09-21
  dNBR rata-rata di bekas terbakar acuan 0.42; di hutan stabil 0.02; acuan 23 piksel (0.23 ha)
  ambang   terdeteksi(ha)  tepat salah lewat  presisi kepekaan
   0.05       0.30        23     7    0    0.77    1.00
   0.10       0.23        23     0    0    1.00    1.00
   0.20       0.22        22     0    1    1.00    0.96
   0.27       0.19        19     0    4    1.00    0.83
   0.44       0.09         9     0   14    1.00    0.39
KPH Alpha  terbakar terdeteksi 0.01 ha (acuan 0.00 ha)
KPH Beta   terbakar terdeteksi 1.28 ha (acuan 1.28 ha)
KPH Gamma  terbakar terdeteksi 0.23 ha (acuan 0.23 ha)

Baca per kejadian. Pada ambang 0,10, kejadian 2024 mendapat 128 piksel tepat, hanya satu salah, dan tidak ada yang terlewat. Kejadian 2025 mendapat 23 piksel, semuanya tepat. Makin tinggi ambang, makin banyak yang terlewat: pada 0,44 kepekaan tinggal 0,23 dan 0,39. Angka uji 0,27 dan 0,44 berasal dari kelas keparahan literatur yang dimuat di Buku M3, dan literatur itu berasal dari hutan di wilayah lain. Pada data ini, ambang terbaik justru 0,10. Kalibrasi ambang dengan data Anda sendiri. Baris KPH Alpha mencatat satu piksel salah tambah (0,01 hektare), karena Alpha tidak punya bekas terbakar sama sekali.

Bagian yang boleh diubah.

• Baris 12, AMBANG_BAKAR. Saya mengubahnya menjadi 0.20. Layar menampilkan KPH Beta terbakar terdeteksi 1.14 ha (acuan 1.28 ha) dan KPH Gamma terbakar terdeteksi 0.22 ha (acuan 0.23 ha), dan Alpha kini 0.00 ha. Satu piksel salah di Alpha hilang, tetapi 14 piksel Beta terlewat.
• Baris 34, nomor tanggal. Pilih tanggal sebelum dan sesudah untuk kejadian Anda. Nomor mulai dari 0.
• Baris 51, syarat layak. Longgarkan atau ketatkan batas 0,5 dan kelas tutupan.
• Baris 56, daftar ambang uji. Tambah atau kurangi angkanya.

Jebakan.

• Rumus di baris 20 sangat panjang. Jumlah %s harus sama dengan jumlah isian. Saya mengubah * 7 menjadi * 6, dan layar menampilkan TypeError: not enough arguments for format string. Pesan ini bagus: ia langsung menunjuk masalahnya.
• Banyak hal dikunci. Nomor tanggal di baris 34, ukuran 40 kali 40 di baris 37, pembagian KPH di baris 38, dan luas 0,01 hektare per piksel di baris 54, 59, dan 64. Semuanya hanya benar untuk data contoh.
• Nilai -9999 dianggap kosong. Baris 48 menganggap setiap dNBR di bawah -100 sebagai kosong. Bila data Anda memakai tanda kosong lain, ganti.
• Variabel a, b, dan k dipakai ulang. Nama a dan b di baris 41 berbeda artinya dari baris 18. Python mengizinkannya, tetapi mudah membingungkan.
• Pembaca dari terminal. Skrip ini memanggil alat QGIS 30 kali, jadi dari terminal butuh persiapan Bab 3.

Latihan variasi.

1. Naikkan AMBANG_BAKAR menjadi 0.20.
2. Bandingkan luas terbakar per KPH dengan hasil asli.

Contoh jawaban, hanya baris 12:

AMBANG_BAKAR = 0.20          # dNBR di atas ini dianggap bekas terbakar (lihat uji ambang di akhir keluaran)

Tiga baris terakhir layar menjadi KPH Alpha terbakar terdeteksi 0.00 ha (acuan 0.00 ha), KPH Beta terbakar terdeteksi 1.14 ha (acuan 1.28 ha), dan KPH Gamma terbakar terdeteksi 0.22 ha (acuan 0.23 ha). Hasil asli (0,10) adalah 0,01; 1,28; dan 0,23 hektare. Ambang yang lebih tinggi menghilangkan salah tambah, tetapi menambah yang terlewat.

Skrip 26.11: Mencocokkan titik panas dan memperkirakan tanggal terbakar #

Berkas 5_2_titik_panas.py. Di dalam berkasnya tertulis "SKRIP 5.2".

Tujuan. Mengubah peta bekas bakar menjadi poligon, mencocokkan titik panas dengan poligon itu (dengan dan tanpa penyangga), memperkirakan tanggal terbakar tiap piksel, dan membandingkannya dengan tanggal titik panas.

Bahan. bekas_bakar.tif dari Skrip 26.10, vektor/Hotspot_Bertanggal.gpkg (12 titik panas, kolom ID_Titik dan TANGGAL), dan m3_umum.py. Jalankan dari Konsol atau terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 5.2: Cocokkan titik panas dengan bekas terbakar, perkirakan tanggal terbakar, hitung luas per KPH
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 5.1 sudah dijalankan (hasil/bekas_bakar.tif).
import datetime as dt
import os
import numpy as np
import processing
from qgis.core import QgsVectorLayer
import m3_umum as U

bakar = os.path.join(U.HASIL, "bekas_bakar.tif")
poli = os.path.join(U.HASIL, "bekas_bakar.gpkg")
if os.path.exists(poli):
    os.remove(poli)                  # polygonize ke GPKG yang sudah ada MENAMBAH fitur, bukan menimpa
processing.run("gdal:polygonize", {"INPUT": bakar, "BAND": 1, "FIELD": "DN", "EIGHT_CONNECTEDNESS": True, "OUTPUT": poli})
lp = QgsVectorLayer(poli + "|layername=bekas_bakar", "bekas_bakar")
lp.setSubsetString('"DN" = 1')                                       # hanya gumpalan terbakar (bukan latar)
print("Poligon bekas terbakar:", lp.featureCount())
hs = QgsVectorLayer(os.path.join(U.PAKET, "vektor", "Hotspot_Bertanggal.gpkg"), "hotspot")


def dalam(jarak):
    sumber = hs
    if jarak > 0:
        sumber = processing.run("native:buffer", {"INPUT": hs, "DISTANCE": jarak, "SEGMENTS": 8, "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    r = processing.run("native:extractbylocation", {"INPUT": sumber, "PREDICATE": [0], "INTERSECT": lp, "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    return sorted(f["ID_Titik"] for f in r.getFeatures())


for j in (0, 15, 30):
    ids = dalam(j)
    print("Titik panas di dalam bekas terbakar (penyangga %d m): %d dari %d -> %s" % (j, len(ids), hs.featureCount(), ", ".join(ids)))

# perkiraan tanggal terbakar per piksel: pengamatan bersih pertama dengan NBR turun lebih dari 0.2 dari dasar (3 pengamatan awal bersih)
tgl, refl, scl, gt, prj = U.baca_deret()
nbr = U.nilai_nbr(refl)
ok = U.bersih(scl)
b, _, _ = U.baca_tif(bakar)
tanggal_bakar = {}
for i, j in zip(*np.where(b == 1)):
    idx = np.where(ok[:, i, j])[0]
    dasar = np.median(nbr[idx[:3], i, j])
    for t in idx[3:]:
        if nbr[t, i, j] < dasar - 0.2:
            tanggal_bakar[(i, j)] = tgl[t]
            break
from collections import Counter
print("Perkiraan tanggal terbakar (jumlah piksel):", dict(sorted(Counter(tanggal_bakar.values()).items())))
# hubungan dengan titik panas
for f in hs.getFeatures():
    p = f.geometry().asPoint()
    i, j = U.xy_ke_piksel(gt, p.x(), p.y())
    tg_h = dt.date.fromisoformat(f["TANGGAL"])
    tb = tanggal_bakar.get((i, j))
    print("%s %s -> %s" % (f["ID_Titik"], tg_h, ("bekas terbakar, tanggal perkiraan %s (selisih %d hari)" % (tb, (tb - tg_h).days)) if tb else "tidak ada bekas terbakar di pikselnya"))

Penjelasan baris demi baris. Baris 11, 21-22, 29-30, dan 34 kosong.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Poligon bekas terbakar: 4
Titik panas di dalam bekas terbakar (penyangga 0 m): 4 dari 12 -> HS-006, HS-009, HS-010, HS-011
Titik panas di dalam bekas terbakar (penyangga 15 m): 7 dari 12 -> HS-002, HS-004, HS-005, HS-006, HS-009, HS-010, HS-011
Titik panas di dalam bekas terbakar (penyangga 30 m): 8 dari 12 -> HS-002, HS-004, HS-005, HS-006, HS-009, HS-010, HS-011, HS-012
Perkiraan tanggal terbakar (jumlah piksel): {datetime.date(2024, 9, 16): 122, datetime.date(2024, 10, 16): 2, datetime.date(2025, 8, 17): 23}
HS-001 2024-10-27 -> tidak ada bekas terbakar di pikselnya
HS-002 2025-08-11 -> tidak ada bekas terbakar di pikselnya
HS-003 2024-10-27 -> tidak ada bekas terbakar di pikselnya
HS-004 2024-09-03 -> tidak ada bekas terbakar di pikselnya
HS-005 2024-09-05 -> tidak ada bekas terbakar di pikselnya
HS-006 2024-09-04 -> bekas terbakar, tanggal perkiraan 2024-09-16 (selisih 12 hari)
HS-007 2024-10-28 -> tidak ada bekas terbakar di pikselnya
HS-008 2025-03-14 -> tidak ada bekas terbakar di pikselnya
HS-009 2024-09-04 -> bekas terbakar, tanggal perkiraan 2024-09-16 (selisih 12 hari)
HS-010 2025-08-12 -> bekas terbakar, tanggal perkiraan 2025-08-17 (selisih 5 hari)
HS-011 2024-09-04 -> bekas terbakar, tanggal perkiraan 2024-09-16 (selisih 12 hari)
HS-012 2025-03-14 -> tidak ada bekas terbakar di pikselnya

Dengan penyangga 0 meter, hanya 4 dari 12 titik panas jatuh di dalam bekas terbakar. Dengan penyangga 30 meter, 8 titik. Empat titik terakhir tetap di luar. Saya memeriksa kelas acuan di piksel tiap titik: HS-001 dan HS-007 jatuh di piksel deforestasi, dan bertanggal akhir Oktober 2024, saat bukaan Gamma terjadi (28 Oktober 2024). Jadi titik panas tidak selalu berarti kebakaran. Pembukaan lahan bisa ikut menyalakannya [kemungkinan]. Tanggal perkiraan terbakar selalu sesudah tanggal titik panas (selisih 12 dan 5 hari), karena tanggal itu adalah citra bersih pertama yang menunjukkan penurunan, sedangkan titik panas dideteksi saat api menyala.

Bagian yang boleh diubah.

• Baris 31, daftar penyangga. Saya mengganti menjadi (0, 50): layar menampilkan penyangga 50 m: 9 dari 12 (ID HS-002, 004, 005, 006, 008, 009, 010, 011, 012).
• Baris 45, 0.2. Ambang penurunan NBR untuk tanggal terbakar.
• Baris 43 dan 44, angka 3. Jumlah pengamatan dasar. Jangan lebih besar dari jumlah pengamatan bersih yang ada.

Jebakan.

• Alat Polygonize menambah, bukan menimpa. Itu sebabnya baris 14-15 menghapus bekas_bakar.gpkg lebih dulu. Saya menjalankan skrip dua kali di folder yang sama, dan kedua putaran mencetak Poligon bekas terbakar: 4. Tanpa dua baris itu, putaran kedua mencetak 8, karena poligon baru ditumpuk di atas yang lama. Bila Windows menolak penghapusan karena berkas itu masih terbuka sebagai layer di QGIS, os.remove gagal dengan PermissionError: [WinError 32] (diuji; layer di kanvas jendela QGIS tidak diuji). Hapus layernya dari proyek, lalu jalankan ulang.
• Skrip 26.10 harus sudah jalan. Di folder kosong, layar menampilkan QgsProcessingException: Unable to execute algorithm dan Could not load source layer for INPUT: ...bekas_bakar.tif not found.
• Nama j dipakai dua kali. Baris 31 dan baris 41 sama-sama memakai j. Keduanya berurutan, jadi aman. Namun menamai ulang salah satunya mencegah kekacauan bila skrip dikembangkan.
• Penyangga yang terlalu lebar tidak bermakna. Area contoh hanya 400 meter kali 400 meter. Saya mencoba (0, 100, 200): layar menampilkan 4, 10, dan 12 dari 12. Pada 200 meter semua titik "cocok", tanpa arti.
• Butuh minimal tiga pengamatan bersih. Baris 43 dan 44 mengandaikan idx sekurang-kurangnya sepanjang tiga.

Latihan variasi.

1. Ganti penyangga menjadi (0, 50).
2. Catat titik panas yang baru masuk pada 50 meter, dibandingkan 30 meter.

Contoh jawaban, hanya baris 31:

for j in (0, 50):

Layar menampilkan penyangga 0 m: 4 dari 12 dan penyangga 50 m: 9 dari 12 -> HS-002, HS-004, HS-005, HS-006, HS-008, HS-009, HS-010, HS-011, HS-012. Pada 30 meter daftarnya delapan ID. Titik baru pada 50 meter: HS-008. Setelah itu coba ingat: HS-008 bertanggal 14 Maret 2025, sedangkan kebakaran di data ini terjadi September 2024 dan Agustus 2025. Cocok dalam jarak, tidak cocok dalam waktu.

Kelompok F: Mesin yang belajar dari contoh #

Konsep: panitia juri dalam tiga kalimat #

Sebuah pohon keputusan itu seperti permainan tebak-tebakan bercabang: "NBR minimumnya di bawah 0,2? Ya, maka terbakar." Sebuah hutan acak adalah ratusan pohon yang masing-masing belajar dari contoh dan fitur yang sedikit berbeda, lalu memberi suara. Kelas dengan suara terbanyak menang, seperti keputusan panitia yang lebih dapat dipercaya daripada satu juri.

Lima istilah yang perlu Anda kenal sekarang:

• Titik latih: titik yang kelasnya sudah diketahui, dipakai mesin untuk belajar. Paket data punya 200 titik, 25 untuk tiap dari delapan kelas.
• Pohon keputusan: aturan bercabang "jika... maka...".
• Hutan acak (random forest): banyak pohon yang memberi suara.
• Uji silang: melatih mesin pada sebagian titik, lalu mengujinya pada titik yang tidak ikut dilatih.
• Kepentingan fitur: seberapa sering sebuah fitur dipakai pohon untuk memisahkan kelas.

Di bab ini, bahan belajar mesin adalah dua belas fitur per piksel, yang dibuat dari riwayat NDVI, NBR, dan SWIR1. Itu pekerjaan Skrip 26.12.

Skrip 26.12: Menyusun dua belas fitur #

Berkas 6_1_fitur_deret.py. Di dalam berkasnya tertulis "SKRIP 6.1".

Tujuan. Menyusun 12 fitur ringkas per piksel dari deret NDVI, NBR, dan SWIR1 yang sudah ditambal, sebagai bahan machine learning.

Bahan. Citra mentah lewat m3_umum.py. Skrip ini tidak memakai QGIS dan tidak memerlukan skrip lain. Jalankan dari Konsol atau terminal.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 6.1: Bangun tumpukan fitur (12 band) dari deret waktu untuk machine learning
# Penulis: Badar Mubarok Yogaswara
# Fitur = angka ringkas yang menggambarkan "perilaku" tiap piksel sepanjang waktu. Hasil: hasil/fitur_deret.tif
import os
import numpy as np
import m3_umum as U

tgl, refl, scl, gt, prj = U.baca_deret()
hari = np.array([(t - tgl[0]).days for t in tgl], dtype="float64")
baik = U.bersih(scl)
ndvi = U.isi_waktu(U.nilai_ndvi(refl), baik, hari)           # NDVI tanpa lubang awan
nbr = U.isi_waktu(U.nilai_nbr(refl), baik, hari)
b11 = U.isi_waktu(refl[:, U.BAND["B11"]], baik, hari)
hc = hari - hari.mean()


def tren_tahun(x):
    return (hc[:, None, None] * (x - x.mean(0))).sum(0) / (hc ** 2).sum() * 365.0


def awal_akhir(x):
    return np.median(x[-4:], axis=0) - np.median(x[:4], axis=0)      # median 4 tanggal terakhir dikurangi 4 tanggal pertama


fitur = {
    "ndvi_rata": ndvi.mean(0), "ndvi_min": ndvi.min(0), "ndvi_maks": ndvi.max(0), "ndvi_simpang": ndvi.std(0),
    "ndvi_tren": tren_tahun(ndvi), "ndvi_turun_maks": np.diff(ndvi, axis=0).min(0), "ndvi_akhir_awal": awal_akhir(ndvi),
    "nbr_rata": nbr.mean(0), "nbr_min": nbr.min(0), "nbr_turun_maks": np.diff(nbr, axis=0).min(0),
    "nbr_akhir_awal": awal_akhir(nbr), "swir1_rata": b11.mean(0),
}
nama = list(fitur)
U.tulis_tif(os.path.join(U.HASIL, "fitur_deret.tif"), np.array([fitur[n] for n in nama]).astype("float32"), gt, prj, deskripsi=nama)
print("Fitur (%d band): %s" % (len(nama), ", ".join(nama)))

Penjelasan baris demi baris.

Arti dua belas fitur itu:

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Fitur (12 band): ndvi_rata, ndvi_min, ndvi_maks, ndvi_simpang, ndvi_tren, ndvi_turun_maks, ndvi_akhir_awal, nbr_rata, nbr_min, nbr_turun_maks, nbr_akhir_awal, swir1_rata

Berkas fitur_deret.tif kini berisi 12 band. Urutan nama di layar adalah urutan band.

Bagian yang boleh diubah.

• Baris 26-31, isi kamus. Tambah fitur baru dengan satu entri. Saya menambahkan "ndvi_median": np.median(ndvi, axis=0) di akhir kamus (setelah swir1_rata). Layar menampilkan Fitur (13 band): ... swir1_rata, ndvi_median. Skrip 26.13 tetap berjalan: Titik latih: 200, fitur: 13.
• Baris 23, angka 4. Jumlah tanggal di awal dan akhir yang dibandingkan.

Jebakan.

• Urutan fitur tidak boleh diubah di tengah. Skrip 26.14 memakai tujuh fitur pertama (range(7)), yaitu fitur NDVI. Bila Anda menyisipkan fitur di depan, "tujuh fitur NDVI" itu berubah isi tanpa pesan. Tambahkan fitur baru di akhir.
• Rumus tren sedikit berbeda dari Skrip 26.7. Di sini memakai x.mean(0) di dalam fungsi, sedangkan di Skrip 26.7 memakai larik rata yang sudah dihitung. Penulisannya tidak persis sama, jadi periksa bila Anda menyalinnya.
• Penambalan per piksel lambat, seperti Skrip 26.6. Tiga kali pemanggilan U.isi_waktu menambah waktu.

Latihan variasi.

1. Tambahkan fitur ndvi_median (median NDVI).
2. Pastikan jumlah band menjadi 13.

Contoh jawaban, hanya baris 30:

    "nbr_akhir_awal": awal_akhir(nbr), "swir1_rata": b11.mean(0), "ndvi_median": np.median(ndvi, axis=0),

Layar menampilkan Fitur (13 band): ndvi_rata, ndvi_min, ... swir1_rata, ndvi_median.

Kotak perkakas kedua: hutan_mini.py #

Skrip berikutnya melatih hutan acak. Pustaka umumnya adalah scikit-learn, tetapi tidak semua QGIS memasangnya. Di QGIS 4.0.2 milik saya, import sklearn menghasilkan ModuleNotFoundError: No module named 'sklearn'. Karena itu penulis membuat hutan_mini.py: hutan acak ringkas yang hanya memakai NumPy.

Berkas ini memuat kelas buatan sendiri. Kelas adalah cetakan benda: ia punya data (disebut atribut) dan pekerjaan (disebut metode) sekaligus. Ada dua kelas di sini. _Pohon mewakili satu pohon keputusan. RandomForestClassifier mewakili hutan yang memegang banyak pohon. Nama terakhir sengaja sama dengan nama di scikit-learn, sehingga skrip lain bisa berganti pustaka tanpa mengubah pemanggilannya. Cara membuat kelas, arti self, dan __init__ dibahas di Bab 11.

Bagian yang paling menarik adalah rekursi: fungsi yang memanggil dirinya sendiri. Metode _simpul membagi data menjadi dua kelompok, lalu memanggil _simpul lagi untuk tiap kelompok. Pembagian berhenti di daun, yaitu kelompok yang kelasnya sudah murni atau datanya terlalu sedikit. Bayangkan Anda membagi tumpukan kartu menjadi dua, lalu membagi tiap tumpukan menjadi dua lagi, sampai tiap tumpukan berisi satu jenis kartu.

Tujuan. Menyediakan hutan acak ringkas bagi Skrip 26.13 dan 26.14, sekaligus membuka "kotak hitam" supaya Anda paham cara kerjanya.

Bahan. Berkas ini tidak dijalankan sendiri. Skrip 26.13 dan 26.14 mengimpornya bila sklearn tidak ada. Berkas harus satu folder dengan mereka. Ia hanya memakai numpy.

Skrip utuh:

# -*- coding: utf-8 -*-
"""hutan_mini.py: hutan acak (random forest) ringkas dengan NumPy saja. Penulis: Badar Mubarok Yogaswara
Tujuan: (1) membuka "kotak hitam" supaya Anda paham cara kerjanya, (2) tetap bisa jalan di QGIS yang belum memasang scikit-learn.
Antarmukanya sengaja dibuat mirip scikit-learn: RandomForestClassifier(...).fit(X, y), .predict(X), .feature_importances_.
Ini alat belajar; untuk pekerjaan besar pakai scikit-learn."""
import numpy as np


class _Pohon:
    def __init__(self, maks_fitur, kedalaman, daun_min, rng):
        self.mf, self.kd, self.dm, self.rng = maks_fitur, kedalaman, daun_min, rng
        self.fitur, self.ambang, self.kiri, self.kanan, self.nilai = [], [], [], [], []
        self.penting = None

    def _gini(self, hitung):
        n = hitung.sum()
        return 1.0 - ((hitung / n) ** 2).sum() if n else 0.0

    def _simpul(self, X, y, K, kdl):
        hitung = np.bincount(y, minlength=K).astype(float)
        idx = len(self.fitur)
        self.fitur.append(-1); self.ambang.append(0.0); self.kiri.append(-1); self.kanan.append(-1); self.nilai.append(hitung / hitung.sum())
        if kdl == self.kd or len(y) < 2 * self.dm or (hitung > 0).sum() == 1:
            return idx
        terbaik = (0.0, None, None)
        g0 = self._gini(hitung)
        n = len(y)
        for f in self.rng.choice(X.shape[1], size=self.mf, replace=False):
            urut = np.argsort(X[:, f], kind="stable")
            xs, ys = X[urut, f], y[urut]
            kiri = np.cumsum(np.eye(K)[ys], axis=0)             # jumlah kelas di sisi kiri untuk tiap titik potong
            kanan = hitung - kiri
            nk = np.arange(1, n + 1)[:, None]
            nr = n - nk
            gk = 1 - ((kiri / nk) ** 2).sum(1)
            gr = np.where(nr[:, 0] > 0, 1 - ((kanan / np.maximum(nr, 1)) ** 2).sum(1), 0)
            skor = g0 - (nk[:, 0] * gk + nr[:, 0] * gr) / n      # penurunan ketidakmurnian
            sah = (xs[:-1] < xs[1:]) & (nk[:-1, 0] >= self.dm) & (nr[:-1, 0] >= self.dm)
            if not sah.any():
                continue
            s = np.where(sah, skor[:-1], -1)
            k = int(s.argmax())
            if s[k] > terbaik[0]:
                terbaik = (s[k], f, (xs[k] + xs[k + 1]) / 2)
        if terbaik[1] is None:
            return idx
        gain, f, t = terbaik
        self.penting[f] += gain * n
        m = X[:, f] <= t
        self.fitur[idx], self.ambang[idx] = f, t
        self.kiri[idx] = self._simpul(X[m], y[m], K, kdl + 1)
        self.kanan[idx] = self._simpul(X[~m], y[~m], K, kdl + 1)
        return idx

    def latih(self, X, y, K):
        self.penting = np.zeros(X.shape[1])
        self._simpul(X, y, K, 0)
        self.fitur, self.ambang = np.array(self.fitur), np.array(self.ambang)
        self.kiri, self.kanan, self.nilai = np.array(self.kiri), np.array(self.kanan), np.array(self.nilai)

    def proba(self, X):
        simpul = np.zeros(len(X), dtype=int)
        while True:
            daun = self.fitur[simpul] < 0
            if daun.all():
                break
            f = self.fitur[simpul]
            ke_kiri = X[np.arange(len(X)), np.maximum(f, 0)] <= self.ambang[simpul]
            baru = np.where(ke_kiri, self.kiri[simpul], self.kanan[simpul])
            simpul = np.where(daun, simpul, baru)
        return self.nilai[simpul]


class RandomForestClassifier:
    def __init__(self, n_estimators=100, max_features="sqrt", max_depth=None, min_samples_leaf=1, random_state=0):
        self.n_estimators, self.max_features, self.max_depth = n_estimators, max_features, max_depth
        self.min_samples_leaf, self.random_state = min_samples_leaf, random_state

    def fit(self, X, y):
        X = np.asarray(X, dtype=float)
        self.classes_, yk = np.unique(y, return_inverse=True)
        K, (n, d) = len(self.classes_), X.shape
        mf = max(1, int(np.sqrt(d))) if self.max_features == "sqrt" else (d if self.max_features is None else int(self.max_features))
        rng = np.random.default_rng(self.random_state)
        self.pohon_, imp = [], np.zeros(d)
        for _ in range(self.n_estimators):
            b = rng.integers(0, n, n)                            # contoh bootstrap: ambil n baris dengan pengembalian
            p = _Pohon(mf, self.max_depth or 10 ** 6, self.min_samples_leaf, rng)
            p.latih(X[b], yk[b], K)
            self.pohon_.append(p)
            imp += p.penting
        self.feature_importances_ = imp / imp.sum() if imp.sum() else imp
        return self

    def predict_proba(self, X):
        X = np.asarray(X, dtype=float)
        return np.mean([p.proba(X) for p in self.pohon_], axis=0)

    def predict(self, X):
        return self.classes_[self.predict_proba(X).argmax(1)]

Penjelasan baris demi baris. Berkas ini padat, jadi dijelaskan per blok.

Hasil. Anda bisa mencobanya langsung dengan data mainan (bukan bagian berkas): dua fitur, enam titik.

from hutan_mini import RandomForestClassifier
X = [[0.86, 0.02], [0.84, 0.03], [0.88, 0.02], [0.52, 0.15], [0.48, 0.17], [0.55, 0.14]]
y = [1, 1, 1, 2, 2, 2]          # 1 = hutan, 2 = sawah
m = RandomForestClassifier(n_estimators=10, random_state=1).fit(X, y)
print("kelas dikenal:", m.classes_)
print("prediksi      :", m.predict([[0.85, 0.02], [0.50, 0.16]]))
print("kepentingan   :", m.feature_importances_.round(2))

Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

kelas dikenal: [1 2]
prediksi      : [1 2]
kepentingan   : [0.29 0.71]

Titik pertama ditebak hutan, titik kedua sawah. Fitur kedua (simpangan baku) lebih penting daripada fitur pertama pada data mainan ini.

Bagian yang boleh diubah.

• Baris 75, pengaturan hutan. n_estimators (jumlah pohon), max_depth (kedalaman maksimum), min_samples_leaf (anggota minimum tiap daun), random_state (benih). Skrip pemanggil yang menentukannya.
• Baris 83, "sqrt". Jumlah fitur yang dicoba tiap simpul.

Jebakan.

• Hasil belum tentu sama dengan scikit-learn. Meskipun benih sama, cara pengacakannya berbeda. Saya tidak bisa mengujinya karena sklearn tidak terpasang di QGIS saya. Bila Anda memasangnya, angka di Skrip 26.13 dan 26.14 bisa berbeda dari yang tertera di bab ini [CEK].
• Lambat untuk data besar. Ini alat belajar. Untuk pekerjaan besar, pakai scikit-learn.
• Label boleh teks. Saya mencoba label "hutan" dan "sawah" dengan max_depth=1: classes_ menjadi ['hutan' 'sawah'], tebakannya benar, dan pohon pertama hanya punya tiga simpul (satu akar, dua daun).

Latihan variasi.

1. Latih dengan max_depth=1 dan label teks.
2. Hitung simpul pohon pertama.

Contoh jawaban:

y = ["hutan", "hutan", "hutan", "sawah", "sawah", "sawah"]
m = RandomForestClassifier(n_estimators=10, max_depth=1, random_state=1).fit(X, y)
print(m.classes_, m.predict([[0.85, 0.02], [0.50, 0.16]]), len(m.pohon_[0].fitur))

Layar menampilkan ['hutan' 'sawah'] ['hutan' 'sawah'] 3.

Skrip 26.13: Melatih hutan acak dan memetakan seluruh area #

Berkas 6_2_latih_klasifikasi.py. Di dalam berkasnya tertulis "SKRIP 6.2".

Tujuan. Melatih hutan acak dari titik latih, membandingkan uji silang acak dengan uji silang blok ruang, lalu memetakan seluruh area.

Bahan. fitur_deret.tif dari Skrip 26.12, vektor/Titik_Latih.gpkg (kolom KELAS), acuan/Acuan_Dinamika.tif, dan m3_umum.py. Bila sklearn tidak ada, skrip memakai hutan_mini.py. Jalankan dari Konsol atau terminal.

Uji silang perlu dijelaskan sebelum skripnya. Titik latih dibagi menjadi empat kelompok (lipatan). Mesin dilatih pada tiga lipatan, lalu diuji pada satu lipatan yang tersisa. Proses diulang empat kali. Pada lipatan acak, titik dibagi sembarang, sehingga tetangga dekat bisa ada di kedua sisi. Pada lipatan blok, titik dibagi menurut lokasi (empat blok 200 meter), sehingga blok yang diuji benar-benar "asing" bagi mesin.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 6.2: Latih hutan acak dari titik latih, uji silang acak vs uji silang blok ruang, lalu petakan seluruh area
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 6.1 sudah dijalankan. Memakai scikit-learn bila terpasang; bila tidak, memakai hutan_mini.py.
import os
import numpy as np
import processing
from qgis.core import QgsVectorLayer
import m3_umum as U

try:
    from sklearn.ensemble import RandomForestClassifier
    print("Memakai scikit-learn")
except ImportError:
    from hutan_mini import RandomForestClassifier
    print("scikit-learn tidak ada; memakai hutan_mini.py")

fitur_tif = os.path.join(U.HASIL, "fitur_deret.tif")
titik = QgsVectorLayer(os.path.join(U.PAKET, "vektor", "Titik_Latih.gpkg"), "latih")
# ambil nilai 12 fitur pada tiap titik latih (Sample raster values)
s = processing.run("native:rastersampling", {"INPUT": titik, "RASTERCOPY": fitur_tif, "COLUMN_PREFIX": "f", "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
nama_f = [fl.name() for fl in s.fields() if fl.name().startswith("f") and fl.name() != "fid"]
X = np.array([[f[n] for n in nama_f] for f in s.getFeatures()], dtype=float)
y = np.array([f["KELAS"] for f in s.getFeatures()])
xy = np.array([[f.geometry().asPoint().x(), f.geometry().asPoint().y()] for f in s.getFeatures()])
print("Titik latih: %d, fitur: %d, kelas: %s" % (len(y), X.shape[1], dict(zip(*np.unique(y, return_counts=True)))))

# --- uji silang 1: acak 4 lipatan; uji silang 2: lipatan = blok ruang 200 m x 200 m (4 blok, area 400 m)
rng = np.random.default_rng(5)
lipat_acak = rng.permutation(len(y)) % 4
lipat_blok = (xy[:, 0] >= 312200).astype(int) + 2 * (xy[:, 1] >= 9996200).astype(int)


def uji_silang(lipat):
    benar = []
    for k in range(4):
        tr, te = lipat != k, lipat == k
        m = RandomForestClassifier(n_estimators=100, random_state=1).fit(X[tr], y[tr])
        benar.append((m.predict(X[te]) == y[te]).mean())
    return np.array(benar)


a, b = uji_silang(lipat_acak), uji_silang(lipat_blok)
print("Akurasi uji silang ACAK : %s rata-rata %.3f" % (np.round(a, 2), a.mean()))
print("Akurasi uji silang BLOK : %s rata-rata %.3f" % (np.round(b, 2), b.mean()))

# --- latih akhir dan petakan
model = RandomForestClassifier(n_estimators=200, random_state=1).fit(X, y)
urut = np.argsort(model.feature_importances_)[::-1]
from osgeo import gdal
ds = gdal.Open(fitur_tif)
nama_fitur = [ds.GetRasterBand(i + 1).GetDescription() for i in range(ds.RasterCount)]
print("Kepentingan fitur (5 teratas):", ", ".join("%s %.2f" % (nama_fitur[i], model.feature_importances_[i]) for i in urut[:5]))
st, gt, prj = U.baca_tif(fitur_tif)                  # (12, H, W)
nb, H, W = st.shape
peta = model.predict(st.reshape(nb, -1).T).reshape(H, W).astype("uint8")
U.tulis_tif(os.path.join(U.HASIL, "peta_rf.tif"), peta, gt, prj, nodata=0)
din, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
nama = {1: "Hutan alam", 2: "Hutan tanaman", 3: "Pertanian", 4: "Air", 5: "Terbuka", 6: "Deforestasi", 7: "Terbakar", 8: "Panen HTI"}
print("Luas hasil peta (ha) / luas acuan (ha):")
for k in range(1, 9):
    print("  %-14s %6.2f / %6.2f" % (nama[k], (peta == k).sum() * 0.01, (din == k).sum() * 0.01))
print("Kesesuaian peta dengan acuan di SEMUA piksel (bukan uji mandiri): %.3f" % (peta == din).mean())

Penjelasan baris demi baris.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2, tanpa scikit-learn, jadi memakai hutan_mini):

scikit-learn tidak ada; memakai hutan_mini.py
Titik latih: 200, fitur: 12, kelas: {np.int64(1): np.int64(25), np.int64(2): np.int64(25), np.int64(3): np.int64(25), np.int64(4): np.int64(25), np.int64(5): np.int64(25), np.int64(6): np.int64(25), np.int64(7): np.int64(25), np.int64(8): np.int64(25)}
Akurasi uji silang ACAK : [0.98 1.   1.   1.  ] rata-rata 0.995
Akurasi uji silang BLOK : [0.5  0.89 0.63 0.44] rata-rata 0.617
Kepentingan fitur (5 teratas): nbr_rata 0.13, ndvi_rata 0.12, ndvi_maks 0.12, swir1_rata 0.12, nbr_min 0.11
Luas hasil peta (ha) / luas acuan (ha):
  Hutan alam       7.99 /   8.00
  Hutan tanaman    0.96 /   0.96
  Pertanian        1.80 /   1.80
  Air              0.80 /   0.80
  Terbuka          0.43 /   0.43
  Deforestasi      1.30 /   1.30
  Terbakar         1.52 /   1.51
  Panen HTI        1.20 /   1.20
Kesesuaian peta dengan acuan di SEMUA piksel (bukan uji mandiri): 0.999

Baris {np.int64(1): np.int64(25), ...} adalah cara NumPy versi baru menampilkan angka. Artinya tetap: tiap kelas punya 25 titik.

Perhatikan dua baris akurasi. Pada uji silang acak, akurasi 0,995, nyaris sempurna. Pada uji silang blok, akurasi turun jauh menjadi 0,617, dan satu blok hanya 0,44. Mesin tampak hebat saat tetangga dekat boleh ikut belajar, tetapi jauh lebih lemah di tempat yang belum pernah ia lihat. Uji silang acak terlalu optimistis untuk data ruang. Peta akhir sendiri cocok 99,9 persen dengan acuan, tetapi itu bukan ujian mandiri: titik latih ikut dihitung, dan piksel itu juga bagian dari acuan.

Bagian yang boleh diubah.

• Baris 38, n_estimators. Saya mengubahnya menjadi 5: akurasi blok menjadi [0.5 0.62 0.88 0.44] (rata-rata 0,612) dan akurasi acak tetap 0,995. Pohon lebih sedikit, hasil sedikit bergeser.
• Baris 31, batas blok. Saya mengganti 312200 menjadi 312100. Akurasi blok menjadi [0.29 0.92 1. 0.63] dengan rata-rata 0,707. Blok yang berbeda memberi angka berbeda, jadi uji blok tergantung cara memotong.
• Baris 29, benih. Ganti 5 untuk undian lipatan acak yang lain.
• Baris 48, 200 pohon. Lebih banyak pohon, lebih lambat tetapi lebih stabil.

Jebakan.

• Pustaka bisa berganti tanpa Anda sadari. Baris pertama keluaran memberi tahu: scikit-learn tidak ada; memakai hutan_mini.py. Di komputer yang punya scikit-learn, angka akurasi dan peta bisa berbeda [CEK]. Selalu baca baris itu.
• Kolom fitur diambil dari awalan f. Saya mencetak daftar kolom hasil pengambilan nilai: ['fid', 'ID', 'KELAS', 'NAMA_KELAS', 'f1', ..., 'f12']. Daftar nama_f berisi tepat f1 sampai f12. Bila data Anda punya kolom lain berawalan f (huruf kecil), kolom itu ikut terambil dan merusak X.
• Angka empat dan batas blok dikunci. range(4), % 4, dan dua batas koordinat harus cocok. Saya mengganti range(4) menjadi range(3): layar hanya menampilkan tiga angka akurasi, dan satu lipatan tidak pernah diuji. Skrip tidak mengeluh.
• Skrip 26.12 harus sudah jalan. Di folder kosong, layar menampilkan QgsProcessingException: Unable to execute algorithm dan Could not load source layer for RASTERCOPY: ...fitur_deret.tif not found.
• from osgeo import gdal di tengah berkas. Boleh, tetapi biasanya impor ditaruh di atas.
• Luas 0,01 hektare per piksel dikunci di baris 62.

Latihan variasi.

1. Kurangi pohon uji silang menjadi 5, lalu bandingkan akurasi blok.
2. Geser batas blok barat-timur menjadi 312100.

Contoh jawaban, hanya baris yang berubah:

        m = RandomForestClassifier(n_estimators=5, random_state=1).fit(X[tr], y[tr])

lipat_blok = (xy[:, 0] >= 312100).astype(int) + 2 * (xy[:, 1] >= 9996200).astype(int)

Untuk nomor 1, layar menampilkan Akurasi uji silang BLOK : [0.5 0.62 0.88 0.44] rata-rata 0.612. Untuk nomor 2, Akurasi uji silang BLOK : [0.29 0.92 1. 0.63] rata-rata 0.707. Tiap variasi dijalankan terpisah.

Kelompok G: Seberapa boleh peta ini dipercaya? #

Konsep: ujian dengan soal baru #

Validasi itu seperti ujian dengan soal yang belum pernah dilihat murid. Nilai ujian yang soalnya sama dengan bahan belajar tidak membuktikan apa-apa. Anda menarik sampel acak titik di peta, memeriksa kebenarannya dengan data acuan yang lebih tepercaya, lalu menyusun matriks kebingungan. Dari matriks itu keluar akurasi, dan koreksi luas: perkiraan luas yang sudah memperhitungkan kesalahan peta.

Lima istilah yang perlu Anda kenal sekarang:

• Matriks kebingungan: tabel yang menyilangkan kelas peta (baris) dengan kelas acuan (kolom).
• Akurasi pengguna (UA): dari semua piksel yang peta sebut kelas X, berapa bagian yang benar kelas X. Mengukur kesalahan "terlalu banyak menandai".
• Akurasi produsen (PA): dari semua piksel yang sebenarnya kelas X, berapa bagian yang peta temukan. Mengukur kesalahan "terlewat".
• Sampel acak berstrata: sampel dengan jumlah tetap pada tiap kelas peta, supaya kelas langka ikut terwakili.
• Luas terkoreksi: luas yang disesuaikan dengan kesalahan peta, disertai rentang ketidakpastian (di sini 95 persen).

Skrip 26.14: Validasi dua peta dengan sampel berstrata #

Berkas 7_1_validasi.py. Di dalam berkasnya tertulis "SKRIP 7.1".

Tujuan. Menguji dua peta klasifikasi dengan sampel acak berstrata: matriks kebingungan, akurasi keseluruhan, akurasi pengguna dan produsen, luas terkoreksi dengan rentang 95 persen, dan kestabilan terhadap pengundian sampel.

Bahan. fitur_deret.tif (Skrip 26.12), peta_rf.tif (Skrip 26.13), acuan/Acuan_Dinamika.tif, vektor/Titik_Latih.gpkg, dan m3_umum.py. Dua peta diuji. Peta A memakai 12 fitur. Peta B hanya memakai tujuh fitur NDVI. Dalam praktik, acuan berasal dari tafsir citra resolusi tinggi atau survei lapangan. Di data ini, acuannya kunci jawaban.

Skrip ini yang paling panjang di bab ini, dengan dua fungsi besar.

Skrip utuh:

# -*- coding: utf-8 -*-
# SKRIP 7.1: Validasi peta: sampel acak berstrata, matriks kebingungan, akurasi, dan estimasi luas terkoreksi
# Penulis: Badar Mubarok Yogaswara
# Dua peta diuji: (A) hutan acak dengan 12 fitur (hasil Skrip 6.2) dan (B) hutan acak yang hanya memakai fitur NDVI.
# Acuan: Acuan_Dinamika.tif. Dalam praktik, acuan berasal dari tafsir citra resolusi tinggi atau survei lapangan.
import os
import numpy as np
from qgis.core import QgsVectorLayer
import m3_umum as U

try:
    from sklearn.ensemble import RandomForestClassifier
except ImportError:
    from hutan_mini import RandomForestClassifier

NAMA = {1: "Hutan alam", 2: "Hutan tanaman", 3: "Pertanian", 4: "Air", 5: "Terbuka", 6: "Deforestasi", 7: "Terbakar", 8: "Panen HTI"}
N_STRATUM = 20            # sampel per kelas peta
LUAS_PIKSEL = 0.01        # ha

st, gt, prj = U.baca_tif(os.path.join(U.HASIL, "fitur_deret.tif"))
acuan, _, _ = U.baca_tif(os.path.join(U.PAKET, "acuan", "Acuan_Dinamika.tif"))
latih = QgsVectorLayer(os.path.join(U.PAKET, "vektor", "Titik_Latih.gpkg"), "latih")
pl = [(U.xy_ke_piksel(gt, f.geometry().asPoint().x(), f.geometry().asPoint().y()), f["KELAS"]) for f in latih.getFeatures()]
tr = np.array([p[0] for p in pl])
yl = np.array([p[1] for p in pl])
pakai_latih = np.zeros(acuan.shape, bool)
pakai_latih[tr[:, 0], tr[:, 1]] = True

# peta B: hanya 7 fitur NDVI (band 1 sampai 7 pada tumpukan fitur)
idx = list(range(7))
mB = RandomForestClassifier(n_estimators=100, random_state=1).fit(st[idx][:, tr[:, 0], tr[:, 1]].T, yl)
petaB = mB.predict(st[idx].reshape(len(idx), -1).T).reshape(acuan.shape).astype("uint8")
jB = os.path.join(U.HASIL, "peta_ndvi_saja.tif")
U.tulis_tif(jB, petaB, gt, prj, nodata=0)
jA = os.path.join(U.HASIL, "peta_rf.tif")


def sampel_berstrata(peta, seed):
    """Sampel acak berstrata: N_STRATUM piksel acak dari tiap kelas peta, di luar titik latih.
    Padanan di antarmuka QGIS: poligonkan peta, lebur per kelas, lalu Random points in polygons (hasil tiap jalan tidak identik)."""
    rng = np.random.default_rng(seed)
    hasil = {}
    for k in np.unique(peta):
        if k == 0:
            continue
        kand = np.argwhere((peta == k) & ~pakai_latih)
        pilih = rng.choice(len(kand), size=min(N_STRATUM, len(kand)), replace=False)
        hasil[int(k)] = [((int(kand[p][0]), int(kand[p][1])), int(acuan[kand[p][0], kand[p][1]])) for p in pilih]
    return hasil

def analisis(peta, sampel, nama, cetak=True):
    kelas = sorted(sampel)
    pop = peta[~pakai_latih]                          # populasi sampel = piksel di luar titik latih
    A = pop.size * LUAS_PIKSEL
    W = {k: (pop == k).sum() / pop.size for k in kelas}
    n = {k: len(sampel[k]) for k in kelas}
    nij = np.zeros((len(kelas), len(kelas)))
    for a, k in enumerate(kelas):
        for _, r in sampel[k]:
            if r in kelas:
                nij[a, kelas.index(r)] += 1
    p = np.array([[W[k] * nij[a, b] / n[k] for b in range(len(kelas))] for a, k in enumerate(kelas)])
    oa = np.trace(p)
    ua = np.diag(p) / p.sum(1)
    pa = np.diag(p) / np.maximum(p.sum(0), 1e-12)
    # piksel titik latih sudah punya label pasti: luasnya ditambahkan apa adanya
    luas_peta = np.array([(peta == k).sum() * LUAS_PIKSEL for k in kelas])
    luas_adj = p.sum(0) * A + np.array([(acuan[pakai_latih] == k).sum() * LUAS_PIKSEL for k in kelas])
    se = np.array([np.sqrt(sum(W[ki] ** 2 * (nij[a, b] / n[ki]) * (1 - nij[a, b] / n[ki]) / (n[ki] - 1) for a, ki in enumerate(kelas))) for b in range(len(kelas))]) * A
    if cetak:
        print("\n=== %s: %d sampel; baris = kelas peta, kolom = kelas acuan ===" % (nama, nij.sum()))
        print("          " + " ".join("%5d" % k for k in kelas))
        for a, k in enumerate(kelas):
            print("kelas %d   " % k + " ".join("%5d" % v for v in nij[a]) + "   | n=%d" % n[k])
        print("Akurasi keseluruhan (berbobot luas): %.3f" % oa)
        print("%-14s %8s %8s %10s %12s %10s" % ("Kelas", "UA", "PA", "luas peta", "luas terkoreksi", "+/- 95%"))
        for a, k in enumerate(kelas):
            print("%-14s %8.2f %8.2f %9.2f ha %11.2f ha %8.2f ha" % (NAMA[k], ua[a], pa[a], luas_peta[a], luas_adj[a], 1.96 * se[a]))
    return oa, kelas, luas_adj, se


for nama, jalur in (("Peta A (12 fitur)", jA), ("Peta B (NDVI saja)", jB)):
    peta, _, _ = U.baca_tif(jalur)
    oa, kelas, luas_adj, se = analisis(peta, sampel_berstrata(peta, 11), nama)
    print("Luas acuan sebenarnya: " + ", ".join("%s %.2f" % (NAMA[k], (acuan == k).sum() * LUAS_PIKSEL) for k in (6, 7)))
    # kestabilan terhadap pengundian sampel
    oas = [analisis(peta, sampel_berstrata(peta, s), nama, cetak=False)[0] for s in range(21, 26)]
    print("Akurasi keseluruhan pada 5 undian sampel lain: " + ", ".join("%.3f" % v for v in oas))

Penjelasan baris demi baris. Statistiknya sengaja dijelaskan dengan kata-kata dulu, bukan rumus.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2, tanpa scikit-learn):

=== Peta A (12 fitur): 158 sampel; baris = kelas peta, kolom = kelas acuan ===
              1     2     3     4     5     6     7     8
kelas 1      20     0     0     0     0     0     0     0   | n=20
kelas 2       0    20     0     0     0     0     0     0   | n=20
kelas 3       0     0    20     0     0     0     0     0   | n=20
kelas 4       0     0     0    20     0     0     0     0   | n=20
kelas 5       0     0     0     0    18     0     0     0   | n=18
kelas 6       0     0     0     0     0    20     0     0   | n=20
kelas 7       1     0     0     0     0     0    19     0   | n=20
kelas 8       0     0     0     0     0     0     0    20   | n=20
Akurasi keseluruhan (berbobot luas): 0.995
Kelas                UA       PA  luas peta luas terkoreksi    +/- 95%
Hutan alam         1.00     0.99      7.99 ha        8.05 ha     0.12 ha
Hutan tanaman      1.00     1.00      0.96 ha        0.96 ha     0.00 ha
Pertanian          1.00     1.00      1.80 ha        1.80 ha     0.00 ha
Air                1.00     1.00      0.80 ha        0.80 ha     0.00 ha
Terbuka            1.00     1.00      0.43 ha        0.43 ha     0.00 ha
Deforestasi        1.00     1.00      1.30 ha        1.30 ha     0.00 ha
Terbakar           0.95     1.00      1.52 ha        1.46 ha     0.12 ha
Panen HTI          1.00     1.00      1.20 ha        1.20 ha     0.00 ha
Luas acuan sebenarnya: Deforestasi 1.30, Terbakar 1.51
Akurasi keseluruhan pada 5 undian sampel lain: 0.995, 1.000, 1.000, 1.000, 1.000

=== Peta B (NDVI saja): 158 sampel; baris = kelas peta, kolom = kelas acuan ===
              1     2     3     4     5     6     7     8
kelas 1      20     0     0     0     0     0     0     0   | n=20
kelas 2       0    20     0     0     0     0     0     0   | n=20
kelas 3       0     0    20     0     0     0     0     0   | n=20
kelas 4       0     0     0    20     0     0     0     0   | n=20
kelas 5       0     0     0     0    18     0     0     0   | n=18
kelas 6       0     0     0     0     0    20     0     0   | n=20
kelas 7      10     0     0     0     0     0    10     0   | n=20
kelas 8       0     0     0     0     0     0     0    20   | n=20
Akurasi keseluruhan (berbobot luas): 0.915
Kelas                UA       PA  luas peta luas terkoreksi    +/- 95%
Hutan alam         1.00     0.85      6.87 ha        8.07 ha     0.54 ha
Hutan tanaman      1.00     1.00      0.96 ha        0.96 ha     0.00 ha
Pertanian          1.00     1.00      1.80 ha        1.80 ha     0.00 ha
Air                1.00     1.00      0.80 ha        0.80 ha     0.00 ha
Terbuka            1.00     1.00      0.43 ha        0.43 ha     0.00 ha
Deforestasi        1.00     1.00      1.30 ha        1.30 ha     0.00 ha
Terbakar           0.50     1.00      2.64 ha        1.44 ha     0.54 ha
Panen HTI          1.00     1.00      1.20 ha        1.20 ha     0.00 ha
Luas acuan sebenarnya: Deforestasi 1.30, Terbakar 1.51
Akurasi keseluruhan pada 5 undian sampel lain: 0.940, 0.917, 0.874, 0.889, 0.925

Bandingkan dua peta. Peta A hampir sempurna: akurasi keseluruhan 0,995 dan luas terbakar terkoreksi 1,46 hektare, dekat dengan acuan 1,51 hektare. Peta B (hanya NDVI) lebih lemah: akurasi 0,915. Lihat baris kelas 7 pada matriks Peta B: dari 20 sampel yang peta sebut Terbakar, hanya 10 yang benar terbakar dan 10 adalah hutan alam. Akurasi penggunanya 0,50. Luas terbakar menurut peta B adalah 2,64 hektare, sekitar 1,7 kali luas acuan (1,51 hektare). Setelah dikoreksi, luasnya 1,44 hektare dengan rentang plus minus 0,54. Kesimpulan: luas dari peta mentah bisa jauh meleset, dan koreksi luas memperbaikinya. Lima undian ulang menunjukkan Peta B kurang stabil (0,874 sampai 0,940), sedangkan Peta A hampir tidak bergeser.

Bagian yang boleh diubah.

• Baris 17, N_STRATUM. Jumlah sampel per kelas. Saya mengubahnya menjadi 10: sampel Peta A turun menjadi 80 dan akurasi keseluruhan naik menjadi 1,000. Akurasi Peta B menjadi 0,863, dan lima undian ulangnya melebar (0,881 sampai 0,983). Sampel yang lebih sedikit membuat angka lebih goyah.
• Baris 84 dan 87, benih. Ganti benih 11 atau rentang range(21, 26) untuk undian lain.
• Baris 30, range(7). Daftar fitur untuk Peta B. Ganti untuk menguji kelompok fitur lain.

Jebakan.

• Sampel kecil menghasilkan angka nan. Galat baku membagi dengan n - 1. Saya mengubah N_STRATUM menjadi 1: layar tetap menampilkan tabel, tetapi kolom rentang berisi nan ha, tanpa galat. Skrip diam. Bila ada kelas dengan satu sampel, jangan percaya rentangnya.
• Titik latih dikeluarkan dari populasi. Itu sengaja, supaya tidak menguji dengan soal yang sudah dipelajari. Akibatnya, luas terkoreksi menambah luas titik latih "apa adanya" (baris 68). Bentuk ini tidak persis sama dengan rumus bakunya [CEK].
• Rumus statistik tidak dinamai di skrip. Pendekatannya mirip penaksir berstrata yang dibahas Olofsson dan kawan-kawan (2014, Remote Sensing of Environment 148: 42-57). Saya hanya memastikan artikel itu ada, dan belum mencocokkan rumus baris demi baris dengan skrip [CEK].
• Dua peta, satu pengujian. Peta A harus sudah dibuat Skrip 26.13, dan tumpukan fitur harus sudah dibuat Skrip 26.12.
• Nomor fitur dikunci. Baris 30 mengandaikan tujuh fitur pertama adalah fitur NDVI.
• Luas 0,01 hektare per piksel dikunci.

Latihan variasi.

1. Kurangi sampel per kelas menjadi 10.
2. Lihat apakah akurasi Peta B bertambah goyah.

Contoh jawaban, hanya baris 17:

N_STRATUM = 10            # sampel per kelas peta

Layar menampilkan Akurasi keseluruhan (berbobot luas): 0.863 untuk Peta B, dan lima undian ulang 0.881, 0.915, 0.915, 0.983, 0.915. Rentang undian ulang (0,881 sampai 0,983) lebih lebar daripada yang asli (0,874 sampai 0,940). Jadi, Peta B menjadi lebih goyah.

Cek paham #

1. Mengapa grafik NDVI hutan alam di Skrip 26.2 sering jatuh ke nol, padahal hutannya utuh?
2. Mengapa uji silang blok (Skrip 26.13) memberi akurasi jauh lebih rendah daripada uji silang acak?
3. Skrip 26.9 menyimpan tanggal sebagai bilangan seperti 20241115. Mengapa baris 54 harus menyebut tipe=gdal.GDT_Int32, dan apa yang terjadi bila Anda membuang bagian itu?

Jawaban.

1. Awan dan bayangannya ikut terhitung sebagai permukaan, sehingga NDVI jatuh pada tanggal berawan. Hutannya tidak hilang. Kelompok B membuang atau menambal piksel itu.
2. Pada uji acak, tetangga dekat bisa ada di kedua sisi, jadi mesin "mengintip" jawabannya. Pada uji blok, blok yang diuji benar-benar asing bagi mesin. Angka blok lebih jujur untuk data berlokasi.
3. Fungsi U.tulis_tif memberi tipe Byte (0 sampai 255) untuk semua larik non-desimal, jadi 20241115 terpotong menjadi 255. Dengan tipe=gdal.GDT_Int32 di baris 54 Skrip 26.9, berkas menyimpan tanggal utuh. Skrip tidak mengeluh dalam kedua kasus, jadi tipe berkas harus diperiksa sendiri.

Kesalahan umum #

• Menjalankan skrip tanpa skrip pendahulunya. Skrip 26.4 berhenti dengan pesan bila masukannya tidak ada, tetapi tidak semua skrip dijaga begitu. Perbaikan: ikuti tabel urutan di awal bab, dan cocokkan angka keluaran dengan contoh di bab.
• Lupa menghapus hasil lama sebelum mengulang. Alat Polygonize menambah poligon ke berkas yang ada. Skrip 26.8 dan 26.11 sudah menghapus .gpkg lama lebih dulu. Bila Anda menulis skrip sendiri dengan alat ini, tambahkan penghapusan yang sama.
• Mempercayai satu ambang atau satu akurasi. Ambang literatur belum tentu cocok untuk data Anda (Skrip 26.10), dan uji silang acak terlalu manis (Skrip 26.13). Perbaikan: uji beberapa ambang, dan pakai uji blok serta sampel berstrata.
• Membiarkan angka yang dikunci untuk data contoh. Ukuran 40 kali 40, luas 0,01 hektare per piksel, dan batas koordinat tertulis langsung di kode. Perbaikan: gunakan tabel di bawah ini sebagai daftar periksa saat Anda mengganti data.

Ringkasan dan latihan #

• Deret waktu adalah tumpukan citra bertanggal. Tiap piksel punya riwayat angka yang bisa dibaca.
• Pembersihan awan memakai tiga cara: buang (mask), gabung (median), dan tambal (interpolasi).
• Deteksi perubahan punya dua cara: dua waktu (cepat, tanpa tanggal) dan titik patah (menyebut tanggal).
• Kebakaran lebih jelas terlihat di NBR daripada di NDVI. Titik panas adalah petunjuk, bukan peta bekas bakar.
• Mesin belajar butuh fitur. Uji blok memberi angka yang lebih jujur, dan validasi dengan sampel berstrata memberi luas terkoreksi.
• Skrip yang jalan tanpa galat belum tentu benar. Cocokkan angkanya.

Latihan. Gabungkan tiga skrip. Jalankan Skrip 26.1, 26.4, dan 26.5 untuk periode "kemarau2025" (tanggal ke-17 sampai ke-20). Lalu hitung dNDVI antara komposit "sebelum" dan "kemarau2025" dengan pola Skrip 26.8. Bagian mana dari Skrip 26.8 yang harus Anda ganti, dan bagian mana yang boleh dibiarkan?

Tabel perbandingan: pekerjaan, skrip, dan perkakas Python #