# [SKRIP 7.1: Proyek mini terpadu I1-I3: dari kiriman lapangan sampai ringkasan per KPH]
# Penulis: Badar Mubarok Yogaswara
# Keluaran di folder HASIL: Proyek_Terpadu.qgz, Hasil_Terpadu.gpkg, ringkasan_kph.csv, log_kiriman.csv.
import csv
import glob
import os
import processing
from qgis.core import QgsProject, QgsVectorLayer, QgsRasterLayer, QgsCoordinateReferenceSystem

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

CRS = "EPSG:32749"
SEMENTARA = "TEMPORARY_OUTPUT"
proyek = QgsProject.instance()
proyek.removeAllMapLayers()
proyek.setCrs(QgsCoordinateReferenceSystem(CRS))
proyek.setTitle("KPH Contoh - proyek terpadu")


# ------------------------------------------------------------------ bagian 0: pembantu
def muat_vektor(jalur, nama):
    lyr = QgsVectorLayer(jalur, nama, "ogr")
    if not lyr.isValid():
        raise ValueError("vektor tidak bisa dibuka: " + jalur)
    return lyr


def muat_raster(jalur, nama):
    lyr = QgsRasterLayer(jalur, nama, "gdal")
    if not lyr.isValid():
        raise ValueError("raster tidak bisa dibuka: " + jalur)
    return lyr


def jalankan(alat, **parameter):
    """Menjalankan alat Processing lalu mengembalikan isi kunci OUTPUT."""
    return processing.run(alat, parameter)["OUTPUT"]


def hitung_kolom(layer, nama, tipe, rumus):
    """Menambah satu kolom dengan ekspresi. tipe: 0 = desimal, 1 = bilangan bulat."""
    return jalankan("native:fieldcalculator", INPUT=layer, FIELD_NAME=nama, FIELD_TYPE=tipe,
                    FIELD_LENGTH=10, FIELD_PRECISION=2 if tipe == 0 else 0, FORMULA=rumus,
                    OUTPUT=SEMENTARA)


batas = muat_vektor(DATA + "/Batas_KPH.gpkg", "Batas_KPH")
sungai = muat_vektor(DATA + "/Sungai.gpkg", "Sungai")
hotspot = muat_vektor(DATA + "/Hotspot_Nasional.gpkg", "Hotspot_Nasional")
citra = muat_raster(DATA + "/Citra_KPH.tif", "Citra_KPH")
dsm = muat_raster(DATA + "/DSM_Drone.tif", "DSM")
dtm = muat_raster(DATA + "/DTM_Drone.tif", "DTM")
for lyr in (batas, sungai, hotspot, citra):
    proyek.addMapLayer(lyr)

# ------------------------------------------------------------------ bagian I3: data lapangan
tujuan_crs = QgsCoordinateReferenceSystem(CRS)
diterima, log = [], []
for jalur in sorted(glob.glob(DATA + "/Kiriman_Lapangan/*.gpkg")):
    nama = os.path.basename(jalur)
    try:
        lyr = muat_vektor(jalur, nama)
        kurang = [f for f in ("Kondisi", "Tinggi_Phn", "Foto") if lyr.fields().indexOf(f) < 0]
        if kurang:
            raise ValueError("field hilang: " + ", ".join(kurang))
        if lyr.crs().authid() != CRS:
            lyr = jalankan("native:reprojectlayer", INPUT=lyr, TARGET_CRS=tujuan_crs, OUTPUT=SEMENTARA)
        diterima.append(lyr)
        log.append((nama, "DITERIMA"))
    except Exception as e:
        log.append((nama, "DITOLAK: %s" % e))

gabung = jalankan("native:mergevectorlayers", LAYERS=diterima, CRS=tujuan_crs, OUTPUT=SEMENTARA)
# Tiap berkas membawa nomor fid sendiri, jadi setelah digabung ada fid kembar. Buang kolom itu;
# GeoPackage akan memberi nomor baru yang unik saat disimpan.
survei = jalankan("native:deletecolumn", INPUT=gabung, COLUMN=["fid"], OUTPUT=SEMENTARA)
survei.setName("Titik_Survei_Gabungan")
proyek.addMapLayer(survei)
with open(HASIL + "/log_kiriman.csv", "w", newline="", encoding="utf-8") as f:
    csv.writer(f).writerows([("berkas", "status")] + log)
ditolak = sum(1 for _, s in log if s != "DITERIMA")
print("Kiriman:", len(log) - ditolak, "diterima,", ditolak, "ditolak;",
      "titik gabungan =", survei.featureCount())

# ------------------------------------------------------------------ bagian I1: analisis vektor
sempadan = jalankan("native:buffer", INPUT=sungai, DISTANCE=30, DISSOLVE=True, OUTPUT=SEMENTARA)
sempadan.setName("Sempadan_30m")
proyek.addMapLayer(sempadan)
irisan = jalankan("native:intersection", INPUT=batas, OVERLAY=sempadan, INPUT_FIELDS=["NAMA_KPH"],
                  OVERLAY_FIELDS=[], OUTPUT=SEMENTARA)
sempadan_ha = {}
for f in irisan.getFeatures():
    sempadan_ha[f["NAMA_KPH"]] = sempadan_ha.get(f["NAMA_KPH"], 0) + f.geometry().area() / 10000

# Kolom baru lewat ekspresi: luas, titik panas, titik survei (agregasi dari layer lain)
cacah = "aggregate('%s', 'count', \"%s\", intersects($geometry, geometry(@parent)))"
ringkas = hitung_kolom(batas, "LUAS_HA", 0, "round($area / 10000, 2)")
ringkas = hitung_kolom(ringkas, "HOTSPOT", 1, cacah % ("Hotspot_Nasional", "ID_Titik"))
ringkas = hitung_kolom(ringkas, "TITIK_SURVEI", 1, cacah % ("Titik_Survei_Gabungan", "Foto"))
ringkas.setName("Ringkasan_KPH")
proyek.addMapLayer(ringkas)

# ------------------------------------------------------------------ bagian I2: analisis raster
ndvi_jalur = jalankan("native:rastercalc", LAYERS=[citra], OUTPUT=HASIL + "/NDVI_KPH.tif",
                      EXPRESSION='("Citra_KPH@4" - "Citra_KPH@1") / ("Citra_KPH@4" + "Citra_KPH@1")')
chm_jalur = jalankan("native:rastercalc", LAYERS=[dsm, dtm], OUTPUT=HASIL + "/CHM.tif",
                     EXPRESSION='"DSM@1" - "DTM@1"')
ndvi = muat_raster(ndvi_jalur, "NDVI_KPH")
chm = muat_raster(chm_jalur, "CHM")
proyek.addMapLayer(ndvi)
proyek.addMapLayer(chm)


def rata(raster, awalan):
    z = jalankan("native:zonalstatisticsfb", INPUT=batas, INPUT_RASTER=raster, RASTER_BAND=1,
                 COLUMN_PREFIX=awalan, STATISTICS=[2], OUTPUT=SEMENTARA)
    return {f["NAMA_KPH"]: f[awalan + "mean"] for f in z.getFeatures()}


ndvi_rata, chm_rata = rata(ndvi, "n_"), rata(chm, "c_")

# ------------------------------------------------------------------ bagian akhir: tabel dan proyek
judul = ["NAMA_KPH", "LUAS_HA", "SEMPADAN_HA", "HOTSPOT", "TITIK_SURVEI", "NDVI_RATA", "CHM_RATA_M"]
baris = []
for f in ringkas.getFeatures():
    n = f["NAMA_KPH"]
    baris.append([n, f["LUAS_HA"], round(sempadan_ha.get(n, 0), 2), f["HOTSPOT"], f["TITIK_SURVEI"],
                  round(ndvi_rata[n], 3), round(chm_rata[n], 2)])
with open(HASIL + "/ringkasan_kph.csv", "w", newline="", encoding="utf-8") as f:
    w = csv.writer(f)
    w.writerow(judul)
    w.writerows(baris)
print("%-10s %7s %12s %8s %13s %10s %11s" % tuple(judul))
for b in baris:
    print("%-10s %7.2f %12.2f %8d %13d %10.3f %11.2f" % tuple(b))

# Layer sementara hilang saat QGIS ditutup. Simpan yang penting ke satu GeoPackage,
# lalu ganti layer sementara di proyek dengan layer dari GeoPackage itu.
paket = HASIL + "/Hasil_Terpadu.gpkg"
processing.run("native:package", {"LAYERS": [survei, sempadan, ringkas], "OUTPUT": paket,
                                  "OVERWRITE": True})
for lyr in (survei, sempadan, ringkas):
    nama = lyr.name()
    proyek.removeMapLayer(lyr.id())
    proyek.addMapLayer(muat_vektor(paket + "|layername=" + nama, nama))
print("Proyek tersimpan?", proyek.write(HASIL + "/Proyek_Terpadu.qgz"))
print("Isi proyek:", sorted(l.name() for l in proyek.mapLayers().values()))