# [SKRIP 3.2: Batch per fitur: potong raster tinggi tajuk (CHM) untuk tiap KPH]
# Penulis: Badar Mubarok Yogaswara
import os
import processing
from qgis.core import QgsVectorLayer, QgsRasterLayer

DATA = "D:/Latihan/paket-i4"
HASIL = "D:/Latihan/hasil_i4"
folder_keluar = HASIL + "/CHM_per_KPH"
os.makedirs(folder_keluar, exist_ok=True)

batas = QgsVectorLayer(DATA + "/Batas_KPH.gpkg", "Batas_KPH", "ogr")
dsm = QgsRasterLayer(DATA + "/DSM_Drone.tif", "DSM", "gdal")
dtm = QgsRasterLayer(DATA + "/DTM_Drone.tif", "DTM", "gdal")

# 1. CHM untuk seluruh area (satu kali saja, di luar pengulangan)
chm_jalur = processing.run("native:rastercalc", {
    "EXPRESSION": '"DSM@1" - "DTM@1"', "LAYERS": [dsm, dtm], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

# 2. Satu putaran untuk tiap KPH: pilih satu poligon, lalu potong CHM dengannya
for fitur in batas.getFeatures():
    nama = fitur["NAMA_KPH"]                                 # contoh: KPH Alpha
    aman = nama.replace(" ", "_")                            # untuk nama berkas: KPH_Alpha
    satu_kph = processing.run("native:extractbyexpression", {
        "INPUT": batas, "EXPRESSION": '"NAMA_KPH" = \'%s\'' % nama,
        "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    keluar = folder_keluar + "/CHM_" + aman + ".tif"
    processing.run("gdal:cliprasterbymasklayer", {
        "INPUT": chm_jalur, "MASK": satu_kph, "CROP_TO_CUTLINE": True,
        "NODATA": -9999, "OUTPUT": keluar})
    hasil = QgsRasterLayer(keluar, aman, "gdal")
    st = processing.run("native:rasterlayerstatistics", {"INPUT": hasil, "BAND": 1})
    print("%-10s %4d x %4d piksel | rata-rata %.2f m | tertinggi %.2f m" % (
        aman, hasil.width(), hasil.height(), st["MEAN"], st["MAX"]))
print("Berkas di folder:", sorted(os.listdir(folder_keluar)))