# -*- coding: utf-8 -*-
"""M2 Bab 6: kesesuaian lahan multi-kriteria (tumpang tindih berbobot) untuk "jenis pohon contoh" di KPH Contoh.
Kriteria dan ambang di bawah ini FIKTIF untuk latihan; bukan standar kesesuaian lahan resmi.
Penulis: Badar Mubarok Yogaswara. Pemakaian: python-qgis.bat m2_06a_kesesuaian.py <folder paket-m2>
Prasyarat: Bab 4 (hidrologi: Sungai.gpkg) dan Bab 5 (CH_tahunan.tif) sudah dijalankan; keluarannya ada di <paket-m2>/hasil."""
import os
import sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
from osgeo import gdal

gdal.UseExceptions()
paket = sys.argv[1]
out = os.path.join(paket, "hasil")
P = lambda n: os.path.join(out, n)
EXT = "312000,314000,9996000,9998000 [EPSG:32749]"
dem = os.path.join(paket, "DEM_10m.tif")
JARAK_SUNGAI = 30                                            # m: ambang kendala dari sungai (dipakai di rumus QGIS dan pemeriksaan NumPy)


def baca(path):
    ds = gdal.Open(path)
    return ds.GetRasterBand(1).ReadAsArray().astype("float64")


def reklas(inp, tabel, hasil, batas=1):
    flat = [v for baris in tabel for v in baris]
    processing.run("native:reclassifybytable", {"INPUT_RASTER": inp, "RASTER_BAND": 1, "TABLE": flat, "NO_DATA": 0, "RANGE_BOUNDARIES": batas,
                                                "NODATA_FOR_MISSING": True, "DATA_TYPE": 0, "OUTPUT": hasil})


# ---------------------------------------------------------------- 1. bobot dengan AHP (perbandingan berpasangan)
nama = ["lereng", "tanah", "curah hujan", "jarak jalan"]
A = np.array([[1, 2, 3, 4], [1 / 2, 1, 2, 3], [1 / 3, 1 / 2, 1, 2], [1 / 4, 1 / 3, 1 / 2, 1]])
vals, vecs = np.linalg.eig(A)
k = int(np.argmax(vals.real))
w = np.abs(vecs[:, k].real)
w = w / w.sum()
CI = (vals[k].real - len(A)) / (len(A) - 1)
CR = CI / 0.90                                               # RI untuk n = 4 adalah 0,90 (tabel Saaty)
print("bobot AHP:", dict(zip(nama, np.round(w, 3))), "| lambda_maks = %.3f | CI = %.4f | CR = %.4f" % (vals[k].real, CI, CR))

# ---------------------------------------------------------------- 2. kriteria
processing.run("gdal:slope", {"INPUT": dem, "BAND": 1, "SCALE": 1, "AS_PERCENT": True, "COMPUTE_EDGES": True, "ZEVENBERGEN": False,
                              "OUTPUT": P("Lereng_pct.tif")})
reklas(P("Lereng_pct.tif"), [[0, 8, 4], [8, 15, 3], [15, 25, 2], [25, 40, 1]], P("Skor_lereng.tif"))
reklas(P("CH_tahunan.tif"), [[0, 2400, 1], [2400, 2600, 2], [2600, 2800, 3], [2800, 99999, 4]], P("Skor_hujan.tif"))

tanah = processing.run("native:fieldcalculator", {
    "INPUT": os.path.join(paket, "Tanah.gpkg|layername=Tanah"), "FIELD_NAME": "skor", "FIELD_TYPE": 1, "FIELD_LENGTH": 3, "FIELD_PRECISION": 0,
    "FORMULA": "CASE WHEN \"jenis_tanah\"='Aluvial' THEN 4 WHEN \"jenis_tanah\"='Latosol' THEN 3 WHEN \"jenis_tanah\"='Podsolik' THEN 2 ELSE 1 END",
    "OUTPUT": "memory:"})["OUTPUT"]
processing.run("gdal:rasterize", {"INPUT": tanah, "FIELD": "skor", "BURN": 0, "USE_Z": False, "UNITS": 0, "WIDTH": 200, "HEIGHT": 200, "EXTENT": EXT,
                                  "NODATA": 0, "DATA_TYPE": 0, "INIT": 0, "INVERT": False, "OUTPUT": P("Skor_tanah.tif")})


def jarak(vektor, tmp, hasil):
    processing.run("gdal:rasterize", {"INPUT": vektor, "FIELD": "", "BURN": 1, "USE_Z": False, "UNITS": 0, "WIDTH": 200, "HEIGHT": 200,
                                      "EXTENT": EXT, "NODATA": 0, "DATA_TYPE": 0, "INIT": 0, "INVERT": False, "OUTPUT": tmp})
    processing.run("gdal:proximity", {"INPUT": tmp, "BAND": 1, "VALUES": "1", "UNITS": 0, "MAX_DISTANCE": 0, "REPLACE": 0, "NODATA": -9999,
                                      "DATA_TYPE": 5, "OUTPUT": hasil})


jarak(os.path.join(paket, "Jalan.gpkg|layername=Jalan|subset=status<>'ditutup'"), P("_jalan_sel.tif"), P("Jarak_jalan.tif"))
jarak(P("Sungai.gpkg"), P("_sungai_sel.tif"), P("Jarak_sungai.tif"))
reklas(P("Jarak_jalan.tif"), [[0, 100, 4], [100, 250, 3], [250, 500, 2], [500, 1e9, 1]], P("Skor_jalan.tif"))

# ---------------------------------------------------------------- 3. kendala (larangan): lereng >= 40% atau dalam 30 m dari sungai
processing.run("gdal:rastercalculator", {"INPUT_A": P("Skor_lereng.tif"), "BAND_A": 1, "INPUT_B": P("Jarak_sungai.tif"), "BAND_B": 1,
                                         "FORMULA": "logical_and(A>0, B>%d)" % JARAK_SUNGAI, "NO_DATA": 255, "RTYPE": 0, "EXTENT_OPT": 0, "EXTRA": "--hideNoData", "OUTPUT": P("Kendala_boleh.tif")})

# ---------------------------------------------------------------- 4. tumpang tindih berbobot
rumus = "A*%.6f + B*%.6f + C*%.6f + D*%.6f" % tuple(w)
processing.run("gdal:rastercalculator", {"INPUT_A": P("Skor_lereng.tif"), "BAND_A": 1, "INPUT_B": P("Skor_tanah.tif"), "BAND_B": 1,
                                         "INPUT_C": P("Skor_hujan.tif"), "BAND_C": 1, "INPUT_D": P("Skor_jalan.tif"), "BAND_D": 1,
                                         "FORMULA": rumus, "NO_DATA": -9999, "RTYPE": 5, "EXTENT_OPT": 0, "EXTRA": "--hideNoData", "OUTPUT": P("Kesesuaian_skor_mentah.tif")})
processing.run("gdal:rastercalculator", {"INPUT_A": P("Kesesuaian_skor_mentah.tif"), "BAND_A": 1, "INPUT_B": P("Kendala_boleh.tif"), "BAND_B": 1,
                                         "FORMULA": "(B==1)*(1*(A>=3.2) + 2*logical_and(A>=2.6, A<3.2) + 3*logical_and(A>=2.0, A<2.6) + 4*(A<2.0))",
                                         "NO_DATA": 255, "RTYPE": 0, "EXTENT_OPT": 0, "EXTRA": "--hideNoData", "OUTPUT": P("Kesesuaian_kelas.tif")})
kelas = baca(P("Kesesuaian_kelas.tif"))
label = {0: "Terlarang (kendala)", 1: "S1 Sangat sesuai", 2: "S2 Cukup sesuai", 3: "S3 Sesuai marginal", 4: "N Tidak sesuai"}
print("\nLuas per kelas (ha), 1 sel = 100 m2 = 0,01 ha:")
for kd in range(5):
    print("  %-22s %6d sel  %7.2f ha  %5.1f%%" % (label[kd], (kelas == kd).sum(), (kelas == kd).sum() * 0.01, 100.0 * (kelas == kd).mean()))
print("  jumlah: %.2f ha (sel NoData: %d)" % (sum((kelas == kd).sum() for kd in range(5)) * 0.01, (kelas == 255).sum()))

# ---------------------------------------------------------------- 5. pemeriksaan silang dengan NumPy murni
s = [baca(P("Skor_lereng.tif")), baca(P("Skor_tanah.tif")), baca(P("Skor_hujan.tif")), baca(P("Skor_jalan.tif"))]
skor = sum(wi * si for wi, si in zip(w, s))
boleh = (s[0] > 0) & (baca(P("Jarak_sungai.tif")) > JARAK_SUNGAI)
np_kelas = np.where(boleh, np.select([skor >= 3.2, skor >= 2.6, skor >= 2.0], [1, 2, 3], 4), 0)
print("\nPemeriksaan NumPy: sel kelas sama = %.4f%%" % (100.0 * (np_kelas == kelas).mean()))

# ---------------------------------------------------------------- 6. kepekaan: bobot sama rata
skor2 = sum(0.25 * si for si in s)
np2 = np.where(boleh, np.select([skor2 >= 3.2, skor2 >= 2.6, skor2 >= 2.0], [1, 2, 3], 4), 0)
print("Kepekaan: luas S1 bobot AHP = %.2f ha; bobot sama rata = %.2f ha; sel yang berubah kelas = %.1f%%" %
      ((np_kelas == 1).sum() * 0.01, (np2 == 1).sum() * 0.01, 100.0 * (np2 != np_kelas).mean()))

# ---------------------------------------------------------------- 7. per petak (histogram zona) dan poligon untuk basis data
zh = processing.run("native:zonalhistogram", {"INPUT_RASTER": P("Kesesuaian_kelas.tif"), "RASTER_BAND": 1,
                                              "INPUT_VECTOR": os.path.join(paket, "Petak.gpkg|layername=Petak"), "COLUMN_PREFIX": "k_",
                                              "OUTPUT": "memory:"})["OUTPUT"]
baris = []
for f in zh.getFeatures():
    baris.append((f["kode"], f["jenis_tegakan"], *(f["k_%d" % kd] * 0.01 for kd in (1, 2, 3, 4, 0))))
baris.sort(key=lambda r: -r[2])
print("\nLima petak dengan S1 terluas (ha): kode, jenis, S1, S2, S3, N, Terlarang")
for r in baris[:5]:
    print("  %s %-12s %6.2f %6.2f %6.2f %6.2f %6.2f" % r)
print("  total S1 dari 25 petak = %.2f ha" % sum(r[2] for r in baris))
mentah = processing.run("gdal:polygonize", {"INPUT": P("Kesesuaian_kelas.tif"), "BAND": 1, "FIELD": "kelas", "EIGHT_CONNECTEDNESS": False, "EXTRA": "",
                                            "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
if os.path.exists(P("Kesesuaian_kelas.gpkg")):
    os.remove(P("Kesesuaian_kelas.gpkg"))
processing.run("native:dissolve", {"INPUT": mentah, "FIELD": ["kelas"], "SEPARATE_DISJOINT": False, "OUTPUT": P("Kesesuaian_kelas.gpkg")})
for tmp in ("_jalan_sel.tif", "_sungai_sel.tif"):
    try:
        os.remove(P(tmp))
    except OSError:
        pass
print("\npoligon kelas ditulis:", os.path.exists(P("Kesesuaian_kelas.gpkg")))