# Skrip 6.1: menghitung NDVI, mengelompokkannya, dan merata-ratakan per KPH. Penulis: Badar Mubarok Yogaswara
# Jalankan Skrip 2.1 lebih dulu (butuh Citra_4band.tif).
import os
import processing
from qgis.core import QgsRasterLayer, QgsProject
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

DATA = os.environ.get("B3_DATA", r"C:/KPH_Contoh/paket-b3")
OUT = os.environ.get("B3_OUT", r"C:/kerja_b3")
os.makedirs(OUT, exist_ok=True)

stack = QgsRasterLayer(OUT + "/Citra_4band.tif", "Citra_4band")
nir = QgsRasterCalculatorEntry()
nir.ref, nir.raster, nir.bandNumber = "nir@1", stack, 4      # band 4 = inframerah dekat
red = QgsRasterCalculatorEntry()
red.ref, red.raster, red.bandNumber = "red@1", stack, 3      # band 3 = merah

rumus = '("nir@1" - "red@1") / ("nir@1" + "red@1")'
calc = QgsRasterCalculator(rumus, OUT + "/NDVI.tif", "GTiff", stack.extent(), stack.width(), stack.height(),
                           [nir, red], QgsProject.instance().transformContext())
print("Kode hasil (0 = sukses):", calc.processCalculation())

st = processing.run("native:rasterlayerstatistics", {"INPUT": OUT + "/NDVI.tif", "BAND": 1})
print("NDVI min/maks/rata-rata:", round(st["MIN"], 3), "/", round(st["MAX"], 3), "/", round(st["MEAN"], 3))

# Kelompokkan: sampai 0 = air; 0-0,3 = lahan terbuka; 0,3-0,6 = vegetasi jarang; di atas 0,6 = vegetasi rapat
processing.run("native:reclassifybytable", {
    "INPUT_RASTER": OUT + "/NDVI.tif", "RASTER_BAND": 1,
    "TABLE": [-1, 0, 1,  0, 0.3, 2,  0.3, 0.6, 3,  0.6, 1.01, 4],
    "RANGE_BOUNDARIES": 0, "NODATA_FOR_MISSING": True, "NO_DATA": 255, "DATA_TYPE": 0,
    "OUTPUT": OUT + "/NDVI_kelas.tif"})
u = processing.run("native:rasterlayeruniquevaluesreport", {"INPUT": OUT + "/NDVI_kelas.tif", "BAND": 1,
                                                            "OUTPUT_TABLE": "TEMPORARY_OUTPUT"})
nama = {1: "Air", 2: "Lahan terbuka", 3: "Vegetasi jarang", 4: "Vegetasi rapat"}
for f in u["OUTPUT_TABLE"].getFeatures():
    k = int(f["value"])
    print(k, nama.get(k, "NoData"), f["count"], "piksel =", round(f["count"] * 4 / 10000, 2), "ha")

# Rata-rata NDVI per KPH
z = processing.run("native:zonalstatisticsfb", {
    "INPUT": DATA + "/Batas_KPH.gpkg|layername=Batas_KPH", "INPUT_RASTER": OUT + "/NDVI.tif", "RASTER_BAND": 1,
    "COLUMN_PREFIX": "ndvi_", "STATISTICS": [2], "OUTPUT": "TEMPORARY_OUTPUT"})
for f in z["OUTPUT"].getFeatures():
    print(f["NAMA_KPH"], "NDVI rata-rata:", round(f["ndvi_mean"], 3))