# -*- 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))