# -*- coding: utf-8 -*-
# SKRIP 1.1: Susun deret waktu NDVI dari 24 citra Sentinel-2 sintetis (jalankan di Python Console QGIS)
# Penulis: Badar Mubarok Yogaswara
# Hasil: NDVI per tanggal (hasil/ndvi/NDVI_YYYYMMDD.tif) dan satu tumpukan 24 band: hasil/NDVI_deret.tif
import os
import processing
from qgis.core import QgsRasterLayer, QgsProject
from osgeo import gdal
import m3_umum as U

os.makedirs(os.path.join(U.HASIL, "ndvi"), exist_ok=True)
daftar = U.daftar_citra()
keluaran = []
for tgl, pr, ps in daftar:
    nama = "S2_" + tgl.strftime("%Y%m%d") + "_refl"
    lyr = QgsRasterLayer(pr, nama)
    QgsProject.instance().addMapLayer(lyr, False)            # tanpa menampilkan di kanvas
    out = os.path.join(U.HASIL, "ndvi", "NDVI_%s.tif" % tgl.strftime("%Y%m%d"))
    ekspresi = '("%s@4" - "%s@3") / ("%s@4" + "%s@3")' % (nama, nama, nama, nama)   # band 4 = B08, band 3 = B04
    processing.run("native:rastercalc", {"LAYERS": [lyr], "EXPRESSION": ekspresi, "OUTPUT": out})
    keluaran.append(out)

# gabungkan menjadi satu berkas: tiap berkas menjadi satu band (SEPARATE)
stack = os.path.join(U.HASIL, "NDVI_deret.tif")
processing.run("gdal:merge", {"INPUT": keluaran, "SEPARATE": True, "DATA_TYPE": 5, "OUTPUT": stack})
# beri nama band sesuai tanggal
ds = gdal.Open(stack, gdal.GA_Update)
for i, (tgl, _, _) in enumerate(daftar):
    ds.GetRasterBand(i + 1).SetDescription(tgl.isoformat())
ds = None
lyr = QgsRasterLayer(stack, "NDVI_deret")
print("Tumpukan:", stack, "valid:", lyr.isValid(), "band:", lyr.bandCount())