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