# -*- coding: utf-8 -*-
# SKRIP 6.1: Bangun tumpukan fitur (12 band) dari deret waktu untuk machine learning
# Penulis: Badar Mubarok Yogaswara
# Fitur = angka ringkas yang menggambarkan "perilaku" tiap piksel sepanjang waktu. Hasil: hasil/fitur_deret.tif
import os
import numpy as np
import m3_umum as U

tgl, refl, scl, gt, prj = U.baca_deret()
hari = np.array([(t - tgl[0]).days for t in tgl], dtype="float64")
baik = U.bersih(scl)
ndvi = U.isi_waktu(U.nilai_ndvi(refl), baik, hari)           # NDVI tanpa lubang awan
nbr = U.isi_waktu(U.nilai_nbr(refl), baik, hari)
b11 = U.isi_waktu(refl[:, U.BAND["B11"]], baik, hari)
hc = hari - hari.mean()


def tren_tahun(x):
    return (hc[:, None, None] * (x - x.mean(0))).sum(0) / (hc ** 2).sum() * 365.0


def awal_akhir(x):
    return np.median(x[-4:], axis=0) - np.median(x[:4], axis=0)      # median 4 tanggal terakhir dikurangi 4 tanggal pertama


fitur = {
    "ndvi_rata": ndvi.mean(0), "ndvi_min": ndvi.min(0), "ndvi_maks": ndvi.max(0), "ndvi_simpang": ndvi.std(0),
    "ndvi_tren": tren_tahun(ndvi), "ndvi_turun_maks": np.diff(ndvi, axis=0).min(0), "ndvi_akhir_awal": awal_akhir(ndvi),
    "nbr_rata": nbr.mean(0), "nbr_min": nbr.min(0), "nbr_turun_maks": np.diff(nbr, axis=0).min(0),
    "nbr_akhir_awal": awal_akhir(nbr), "swir1_rata": b11.mean(0),
}
nama = list(fitur)
U.tulis_tif(os.path.join(U.HASIL, "fitur_deret.tif"), np.array([fitur[n] for n in nama]).astype("float32"), gt, prj, deskripsi=nama)
print("Fitur (%d band): %s" % (len(nama), ", ".join(nama)))