# -*- coding: utf-8 -*-
# I2 Bab 1: aljabar raster (NDVI, tinggi tajuk, mask). Penulis: Badar Mubarok Yogaswara
# Jalankan di Python Console QGIS. Ubah dua jalur di bawah sesuai komputer Anda.
import os
import processing
from qgis.core import QgsRasterLayer, QgsProcessingFeedback

PAKET = os.environ.get("I2_DATA", r"C:/KPH_Contoh/paket-data")     # folder paket-data
HASIL = r"C:/temp/hasil_i2"                                                      # folder hasil (dibuat bila belum ada)
os.makedirs(HASIL, exist_ok=True)

citra = QgsRasterLayer(os.path.join(PAKET, "Citra_Komposit.tif"), "Citra_Komposit")
dsm = QgsRasterLayer(os.path.join(PAKET, "DSM_Drone.tif"), "DSM_Drone")
dtm = QgsRasterLayer(os.path.join(PAKET, "DTM_Drone.tif"), "DTM_Drone")


def hitung(ekspresi, layers, keluar, ref):
    """Menjalankan Raster calculator dengan ukuran sel, luas, dan CRS mengikuti layer acuan."""
    return processing.run("native:rastercalc", {
        "EXPRESSION": ekspresi, "LAYERS": layers,
        "CELL_SIZE": ref.rasterUnitsPerPixelX(), "EXTENT": ref.extent(), "CRS": ref.crs(),
        "OUTPUT": os.path.join(HASIL, keluar)})["OUTPUT"]


def statistik(jalur, band=1):
    lyr = QgsRasterLayer(jalur, "x")
    s = lyr.dataProvider().bandStatistics(band)
    return s.minimumValue, s.maximumValue, s.mean, s.stdDev


# 1) NDVI = (NIR - R) / (NIR + R); band 1 = R, band 4 = NIR
ndvi = hitung('("Citra_Komposit@4" - "Citra_Komposit@1") / ("Citra_Komposit@4" + "Citra_Komposit@1")',
              [citra], "ndvi.tif", citra)
print("NDVI min, maks, rata-rata, simpangan baku:", [round(v, 3) for v in statistik(ndvi)])

# 2) NDVI > 0,3 menjadi 1 (vegetasi) dan selain itu 0
veg = hitung('"ndvi@1" > 0.3', [QgsRasterLayer(ndvi, "ndvi")], "veg.tif", citra)
print("Rata-rata veg (=bagian piksel bervegetasi):", round(statistik(veg)[2], 4))

# 3) tinggi tajuk (CHM) = DSM - DTM
chm = hitung('"DSM_Drone@1" - "DTM_Drone@1"', [dsm, dtm], "chm.tif", dsm)
print("CHM min, maks, rata-rata:", [round(v, 2) for v in statistik(chm)[:3]])

# 4) tajuk tinggi: tinggi >= 5 m
tinggi = hitung('"chm@1" >= 5', [QgsRasterLayer(chm, "chm")], "tajuk_tinggi.tif", dsm)
print("Bagian piksel dengan tinggi >= 5 m:", round(statistik(tinggi)[2], 4))