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