# -*- coding: utf-8 -*-
"""Mengubah DSM menjadi jaring segitiga 3D (.obj) dan memeriksa pengaturan 3D QGIS. Penulis: Badar Mubarok Yogaswara
Jalankan: python-qgis.bat dsm_ke_obj.py
"""
import os
os.environ.setdefault("QT_QPA_PLATFORM", "offscreen")
from osgeo import gdal
from qgis.core import (QgsApplication, QgsRasterLayer, QgsPointCloudLayer, QgsProject)
gdal.UseExceptions()
AQ = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.abspath(os.path.join(AQ, "..", "..", "..", "data-sintetis", "paket-data"))
OUT = os.path.join(AQ, "keluaran", "tiga_d")


def dsm_ke_obj(tif, obj, langkah=16, skala_z=2.0):
    """Ambil tiap `langkah` piksel; x,y relatif terhadap sudut barat-daya (meter); z dikali skala_z."""
    ds = gdal.Open(tif)
    gt = ds.GetGeoTransform()
    band = ds.GetRasterBand(1)
    data = band.ReadAsArray()
    nodata = band.GetNoDataValue()
    baris = list(range(0, ds.RasterYSize, langkah))
    kolom = list(range(0, ds.RasterXSize, langkah))
    idx = {}
    n = 0
    with open(obj, "w", encoding="utf-8") as fh:
        fh.write("# DSM sintetis KPH Contoh; skala vertikal x%.1f\n" % skala_z)
        for j, r in enumerate(baris):
            for i, c in enumerate(kolom):
                z = float(data[r, c])
                if nodata is not None and z == nodata:
                    z = 0.0
                x = i * langkah * gt[1]
                y = (len(baris) - 1 - j) * langkah * abs(gt[5])    # utara di atas: baris 0 = utara
                n += 1
                idx[(j, i)] = n
                fh.write("v %.3f %.3f %.3f\n" % (x, z * skala_z, y))
        muka = 0
        for j in range(len(baris) - 1):
            for i in range(len(kolom) - 1):
                a, b, c, d = idx[(j, i)], idx[(j, i + 1)], idx[(j + 1, i + 1)], idx[(j + 1, i)]
                fh.write("f %d %d %d\nf %d %d %d\n" % (a, d, c, a, c, b))
                muka += 2
    return n, muka, float(data.min()), float(data.max())


if __name__ == "__main__":
    app = QgsApplication([], False)
    app.initQgis()
    os.makedirs(OUT, exist_ok=True)
    v, f, zmin, zmax = dsm_ke_obj(os.path.join(DATA, "DSM_Drone.tif"), os.path.join(OUT, "dsm.obj"))
    print("OBJ: titik sudut", v, "segitiga", f, "| z min/maks DSM:", round(zmin, 2), round(zmax, 2))
    vs = [l.split()[1:] for l in open(os.path.join(OUT, "dsm.obj"), encoding="utf-8") if l.startswith("v ")]
    xs, ys = [float(a[0]) for a in vs], [float(a[1]) for a in vs]
    print("OBJ: x dari %.0f ke %.0f m | tinggi (sudah x2) dari %.2f ke %.2f" % (min(xs), max(xs), min(ys), max(ys)))
    r = QgsRasterLayer(os.path.join(DATA, "DSM_Drone.tif"), "DSM")
    print("DSM valid:", r.isValid(), "ukuran:", r.width(), "x", r.height(), "CRS:", r.crs().authid())
    pc = QgsPointCloudLayer(os.path.join(DATA, "Hutan_Contoh.laz"), "Hutan", "pdal")
    print("point cloud valid:", pc.isValid(), "jumlah titik:", pc.pointCount())
    from qgis._3d import Qgs3DMapSettings, QgsDemTerrainSettings
    s = Qgs3DMapSettings()
    s.setCrs(r.crs())
    s.setExtent(r.extent())
    t = QgsDemTerrainSettings()
    t.setLayer(r)
    t.setVerticalScale(2.0)
    s.setTerrainSettings(t)
    print("3D: terrain =", s.terrainSettings().type() if hasattr(s.terrainSettings(), "type") else "?",
          "| skala vertikal:", s.terrainSettings().verticalScale(), "| ekstensi:", s.extent().toString(1))
    app.exitQgis()