# [SKRIP 2.3: Alat raster dari Python: NDVI, tinggi tajuk (CHM), dan statistik per KPH]
# Penulis: Badar Mubarok Yogaswara
import processing
from qgis.core import QgsProject, QgsVectorLayer, QgsRasterLayer

DATA = "D:/Latihan/paket-i4"
HASIL = "D:/Latihan/hasil_i4"

batas = QgsVectorLayer(DATA + "/Batas_KPH.gpkg", "Batas_KPH", "ogr")
citra = QgsRasterLayer(DATA + "/Citra_KPH.tif", "Citra_KPH", "gdal")   # band 1=Merah, 4=NIR
dsm = QgsRasterLayer(DATA + "/DSM_Drone.tif", "DSM", "gdal")
dtm = QgsRasterLayer(DATA + "/DTM_Drone.tif", "DTM", "gdal")

# 1. NDVI = (NIR - Merah) / (NIR + Merah). Dalam ekspresi, band ditulis "NamaLayer@nomorBand".
jalur_ndvi = processing.run("native:rastercalc", {
    "EXPRESSION": '("Citra_KPH@4" - "Citra_KPH@1") / ("Citra_KPH@4" + "Citra_KPH@1")',
    "LAYERS": [citra], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
print("Keluaran raster sementara berupa:", type(jalur_ndvi).__name__)    # teks jalur, bukan layer
ndvi = QgsRasterLayer(jalur_ndvi, "NDVI_KPH", "gdal")                    # bungkus jadi layer
print("NDVI sah?", ndvi.isValid())

# 2. CHM = DSM - DTM
jalur_chm = processing.run("native:rastercalc", {
    "EXPRESSION": '"DSM@1" - "DTM@1"', "LAYERS": [dsm, dtm], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
chm = QgsRasterLayer(jalur_chm, "CHM", "gdal")

# 3. Statistik zonal: rata-rata nilai raster di dalam tiap poligon KPH
def rata_rata(raster, awalan):
    return processing.run("native:zonalstatisticsfb", {
        "INPUT": batas, "INPUT_RASTER": raster, "RASTER_BAND": 1,
        "COLUMN_PREFIX": awalan, "STATISTICS": [2],      # 2 = Mean
        "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

zn = rata_rata(ndvi, "ndvi_")
zc = rata_rata(chm, "chm_")
chm_per_kph = {f["NAMA_KPH"]: f["chm_mean"] for f in zc.getFeatures()}
print("KPH         NDVI rata-rata   CHM rata-rata (m)")
for f in zn.getFeatures():
    print("%-11s %12.3f %16.2f" % (f["NAMA_KPH"], f["ndvi_mean"], chm_per_kph[f["NAMA_KPH"]]))

QgsProject.instance().addMapLayer(ndvi)
QgsProject.instance().addMapLayer(chm)