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