# Skrip 7.1: ringkasan per KPH (luas, NDVI rata-rata, luas lereng curam, luas vegetasi rapat). Penulis: Badar Mubarok Yogaswara
# Jalankan Skrip 2.1, 3.1, 5.1, dan 6.1 lebih dulu (butuh NDVI.tif, NDVI_kelas.tif, Lereng_kelas.tif di OUT).
import os
import processing
from qgis.core import QgsVectorLayer

DATA = os.environ.get("B3_DATA", r"C:/KPH_Contoh/paket-b3")
OUT = os.environ.get("B3_OUT", r"C:/kerja_b3")
BATAS = DATA + "/Batas_KPH.gpkg|layername=Batas_KPH"


def luas_kelas(raster, mask, kelas):
    """Luas (ha) piksel bernilai 'kelas' di dalam mask."""
    pot = processing.run("gdal:cliprasterbymasklayer", {
        "INPUT": raster, "MASK": mask, "NODATA": 255, "CROP_TO_CUTLINE": True,
        "KEEP_RESOLUTION": True, "OUTPUT": "TEMPORARY_OUTPUT"})
    u = processing.run("native:rasterlayeruniquevaluesreport", {
        "INPUT": pot["OUTPUT"], "BAND": 1, "OUTPUT_TABLE": "TEMPORARY_OUTPUT"})
    for f in u["OUTPUT_TABLE"].getFeatures():
        if int(f["value"]) == kelas:
            return f["count"] * 4 / 10000
    return 0.0


print("KPH | luas (ha) | NDVI rata-rata | curam (ha) | vegetasi rapat (ha)")
for nama in ("KPH Alpha", "KPH Beta", "KPH Gamma"):
    lyr = QgsVectorLayer(BATAS, "batas", "ogr")
    lyr.setSubsetString("\"NAMA_KPH\" = '%s'" % nama)
    luas = sum(f.geometry().area() for f in lyr.getFeatures()) / 10000
    z = processing.run("native:zonalstatisticsfb", {
        "INPUT": lyr, "INPUT_RASTER": OUT + "/NDVI.tif", "RASTER_BAND": 1,
        "COLUMN_PREFIX": "n_", "STATISTICS": [2], "OUTPUT": "TEMPORARY_OUTPUT"})
    ndvi = [f["n_mean"] for f in z["OUTPUT"].getFeatures()][0]
    curam = luas_kelas(OUT + "/Lereng_kelas.tif", lyr, 3)
    rapat = luas_kelas(OUT + "/NDVI_kelas.tif", lyr, 4)
    print(nama, "|", round(luas, 2), "|", round(ndvi, 3), "|", round(curam, 2), "|", round(rapat, 2))