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