# -*- coding: utf-8 -*-
# I2 Bab 6: uji akurasi (matriks kesalahan, akurasi keseluruhan, kappa). Penulis: Badar Mubarok Yogaswara
# Jalankan setelah skrip bab 4 dan bab 5 (memakai peta hasilnya di folder hasil).
import os
import numpy as np
import processing
from qgis.core import QgsRasterLayer, QgsVectorLayer

PAKET_I2 = os.environ.get("I2_PAKET", r"C:/KPH_Contoh/paket-i2")
HASIL = r"C:/temp/hasil_i2"
NAMA = {1: "Hutan", 2: "Kebun", 3: "Sawah", 4: "Lahan Terbuka"}
PETA = [("Tak terbimbing, 6 kluster", "tak_terbimbing_k6.tif"),
        ("Tak terbimbing, 10 kluster", "tak_terbimbing_k10.tif"),
        ("Terbimbing, jarak minimum", "terbimbing_jarak_min.tif"),
        ("Terbimbing, kemiripan maksimum", "terbimbing_maks_lik.tif")]

titik = QgsVectorLayer(os.path.join(PAKET_I2, "Titik_Uji.gpkg") + "|layername=Titik_Uji", "Titik_Uji", "ogr")
acuan = np.array([f["Kode_Acuan"] for f in titik.getFeatures()])
print("Jumlah titik uji:", len(acuan), "| per kelas:", {NAMA[k]: int((acuan == k).sum()) for k in NAMA})


def matriks(acuan, hasil):
    """Baris = kelas acuan (lapangan), kolom = kelas hasil peta."""
    m = np.zeros((4, 4), int)
    for a, h in zip(acuan, hasil):
        m[a - 1, h - 1] += 1
    return m


def ringkas(m):
    n = m.sum()
    oa = np.trace(m) / n
    pe = (m.sum(axis=1) * m.sum(axis=0)).sum() / n ** 2        # kecocokan yang diharapkan karena kebetulan
    kappa = (oa - pe) / (1 - pe)
    with np.errstate(divide="ignore", invalid="ignore"):       # kolom kosong (kelas tidak pernah muncul di peta) menghasilkan nan
        pa = np.diag(m) / m.sum(axis=1)                        # akurasi pembuat: dari yang benar di lapangan, berapa yang terpetakan benar
        ua = np.diag(m) / m.sum(axis=0)                        # akurasi pengguna: dari yang dipetakan kelas X, berapa yang benar
    return oa, kappa, pa, ua


hasil_semua = {}
for judul, berkas in PETA:
    lyr = QgsRasterLayer(os.path.join(HASIL, berkas), "peta")
    sampel = processing.run("native:rastersampling", {"INPUT": titik, "RASTERCOPY": lyr, "COLUMN_PREFIX": "Peta_", "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    pasang = [(a, f["Peta_1"]) for a, f in zip(acuan, sampel.getFeatures())]
    hasil = np.array([int(h) for a, h in pasang if h is not None])
    acuan_ok = np.array([a for a, h in pasang if h is not None])
    if len(hasil) < len(acuan):
        print("PERINGATAN: %d titik uji di luar raster dilewati" % (len(acuan) - len(hasil)))
    assert hasil.min() >= 1, "Peta berisi kelas 0 (kluster tanpa nama)"
    hasil_semua[judul] = (matriks(acuan_ok, hasil), ringkas(matriks(acuan_ok, hasil)))

print("\nRingkasan:")
print("%-34s | akurasi keseluruhan | kappa" % "Peta")
for judul, (m, (oa, kappa, pa, ua)) in hasil_semua.items():
    print("%-34s | %5.0f%% (%d dari %d)    | %.2f" % (judul, 100 * oa, np.trace(m), m.sum(), kappa))

judul = "Terbimbing, kemiripan maksimum"
m, (oa, kappa, pa, ua) = hasil_semua[judul]
print("\nMatriks kesalahan: %s (baris = acuan, kolom = hasil peta)" % judul)
print("%-14s" % "" + "".join("%14s" % NAMA[k] for k in NAMA) + "   jumlah")
for i, k in enumerate(NAMA):
    print("%-14s" % NAMA[k] + "".join("%14d" % v for v in m[i]) + "%8d" % m[i].sum())
print("%-14s" % "jumlah" + "".join("%14d" % v for v in m.sum(axis=0)))
print("\nKelas          | akurasi pembuat | akurasi pengguna")
for i, k in enumerate(NAMA):
    print("%-14s | %6.0f%%          | %6.0f%%" % (NAMA[k], 100 * pa[i], 100 * ua[i]))

# Cara kedua (alat GRASS di Processing QGIS): r.kappa membandingkan dua raster, yaitu peta hasil dan raster acuan.
# Raster acuan dibuat dari titik uji: piksel yang tidak ada titiknya bernilai 0 (kosong) dan diabaikan.
acuan_tif = os.path.join(HASIL, "acuan_titik.tif")
processing.run("gdal:rasterize", {
    "INPUT": os.path.join(PAKET_I2, "Titik_Uji.gpkg"), "FIELD": "Kode_Acuan", "UNITS": 1, "WIDTH": 1.0, "HEIGHT": 1.0,
    "EXTENT": "312000,312300,9996000,9996300 [EPSG:32749]", "NODATA": 0, "DATA_TYPE": 0, "INIT": 0, "OUTPUT": acuan_tif})
laporan = os.path.join(HASIL, "kappa_maks_lik.txt")
processing.run("grass:r.kappa", {"classification": os.path.join(HASIL, "terbimbing_maks_lik.tif"), "reference": acuan_tif,
                                 "title": "Uji akurasi", "-h": False, "-w": False, "output": laporan})
with open(laporan) as f:
    baris = f.read().splitlines()
print("\nLaporan r.kappa (bagian ringkasan):")
for i, b in enumerate(baris):
    if b.startswith(("Kappa\t", "Obs Correct")):
        print(b)
        print(baris[i + 1])