# -*- coding: utf-8 -*-
# [SKRIP I3.1] Menghitung akurasi dan presisi GNSS dari 100 pembacaan per mode terhadap titik referensi.
# Penulis: Badar Mubarok Yogaswara
# Dijalankan di Python Console QGIS (atau python-qgis.bat). Ubah JALUR bila folder data Anda berbeda.
import math
import os
from qgis.core import QgsVectorLayer

JALUR = r"D:/KPH_Contoh/paket-i3/Gnss_Uji.gpkg"

pembacaan = QgsVectorLayer(JALUR + "|layername=Gnss_Uji", "Gnss_Uji", "ogr")
referensi = QgsVectorLayer(JALUR + "|layername=Titik_Referensi", "Titik_Referensi", "ogr")
assert pembacaan.isValid() and referensi.isValid(), "Layer tidak terbaca: periksa JALUR"

titik_ref = next(referensi.getFeatures()).geometry().asPoint()

# kelompokkan koordinat per mode
per_mode = {}
for fitur in pembacaan.getFeatures():
    p = fitur.geometry().asPoint()
    per_mode.setdefault(fitur["Mode"], []).append((p.x(), p.y()))


def persentil(nilai, q):
    """Persentil dengan interpolasi linear (sama dengan numpy.percentile)."""
    u = sorted(nilai)
    k = (len(u) - 1) * q / 100.0
    a = int(math.floor(k))
    b = min(a + 1, len(u) - 1)
    return u[a] + (u[b] - u[a]) * (k - a)


print("%-15s %8s %8s %8s %8s %8s" % ("Mode", "Bias", "Presisi", "RMSE", "CEP50", "R95"))
for mode, pts in per_mode.items():
    jarak = [math.hypot(x - titik_ref.x(), y - titik_ref.y()) for x, y in pts]       # galat ke titik benar
    cx = sum(x for x, _ in pts) / len(pts)
    cy = sum(y for _, y in pts) / len(pts)
    bias = math.hypot(cx - titik_ref.x(), cy - titik_ref.y())                          # akurasi (geser rata-rata)
    presisi = math.sqrt(sum((x - cx) ** 2 + (y - cy) ** 2 for x, y in pts) / len(pts))  # sebaran di sekitar rata-rata
    rmse = math.sqrt(sum(j * j for j in jarak) / len(jarak))
    print("%-15s %8.3f %8.3f %8.3f %8.3f %8.3f" % (mode, bias, presisi, rmse, persentil(jarak, 50), persentil(jarak, 95)))