BAB 7: Bab Jembatan ke Master: Proyek Mini Terpadu dan Ceklis Masuk M1

# [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()))

Kiriman: 4 diterima, 2 ditolak; titik gabungan = 18
NAMA_KPH   LUAS_HA  SEMPADAN_HA  HOTSPOT  TITIK_SURVEI  NDVI_RATA  CHM_RATA_M
KPH Alpha     4.00         1.18        2             6      0.604        2.60
KPH Beta      4.00         1.39        5             4      0.431        2.69
KPH Gamma     8.00         1.50        5             7      0.329        2.56
Proyek tersimpan? True
Isi proyek: ['Batas_KPH', 'CHM', 'Citra_KPH', 'Hotspot_Nasional', 'NDVI_KPH', 'Ringkasan_KPH', 'Sempadan_30m', 'Sungai', 'Titik_Survei_Gabungan']

# [SKRIP 7.2: Ringkasan per KPH sebagai tugas latar belakang (QgsTask)]
# Penulis: Badar Mubarok Yogaswara
# Pola: pekerjaan berat di fungsi tugas (latar belakang); proyek dan tabel di fungsi selesai (utama).
import csv
import os
import processing
from qgis.core import (QgsApplication, QgsTask, QgsMessageLog, Qgis, QgsProject,
                       QgsVectorLayer, QgsRasterLayer)

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


def jalankan(alat, **parameter):
    return processing.run(alat, parameter)["OUTPUT"]


def hitung_ringkasan(task, data, hasil):
    """Latar belakang. Hanya membaca berkas dan menjalankan alat; tidak menyentuh QgsProject."""
    batas = QgsVectorLayer(data + "/Batas_KPH.gpkg", "Batas_KPH", "ogr")
    sungai = QgsVectorLayer(data + "/Sungai.gpkg", "Sungai", "ogr")
    hotspot = QgsVectorLayer(data + "/Hotspot_Nasional.gpkg", "Hotspot_Nasional", "ogr")
    citra = QgsRasterLayer(data + "/Citra_KPH.tif", "Citra_KPH", "gdal")
    dsm = QgsRasterLayer(data + "/DSM_Drone.tif", "DSM", "gdal")
    dtm = QgsRasterLayer(data + "/DTM_Drone.tif", "DTM", "gdal")

    # tahap 1 (0-25%): sempadan sungai
    sem = jalankan("native:buffer", INPUT=sungai, DISTANCE=30, DISSOLVE=True, OUTPUT=SEMENTARA)
    iris = jalankan("native:intersection", INPUT=batas, OVERLAY=sem, INPUT_FIELDS=["NAMA_KPH"],
                    OVERLAY_FIELDS=[], OUTPUT=SEMENTARA)
    sem_ha = {}
    for f in iris.getFeatures():
        sem_ha[f["NAMA_KPH"]] = sem_ha.get(f["NAMA_KPH"], 0) + f.geometry().area() / 10000
    task.setProgress(25)
    if task.isCanceled():
        return None

    # tahap 2 (25-50%): titik panas
    hit = jalankan("native:countpointsinpolygon", POLYGONS=batas, POINTS=hotspot, FIELD="HOTSPOT",
                   OUTPUT=SEMENTARA)
    titik = {f["NAMA_KPH"]: f["HOTSPOT"] for f in hit.getFeatures()}
    task.setProgress(50)
    if task.isCanceled():
        return None

    # tahap 3 (50-100%): raster
    ndvi = jalankan("native:rastercalc", LAYERS=[citra], OUTPUT=hasil + "/NDVI_KPH.tif",
                    EXPRESSION='("Citra_KPH@4" - "Citra_KPH@1") / ("Citra_KPH@4" + "Citra_KPH@1")')
    chm = jalankan("native:rastercalc", LAYERS=[dsm, dtm], OUTPUT=hasil + "/CHM.tif",
                   EXPRESSION='"DSM@1" - "DTM@1"')
    task.setProgress(75)

    def rata(jalur, awalan):
        raster = QgsRasterLayer(jalur, "r", "gdal")
        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()}

    n, c = rata(ndvi, "n_"), rata(chm, "c_")
    task.setProgress(100)
    hasil_baris = []
    for f in batas.getFeatures():
        k = f["NAMA_KPH"]
        hasil_baris.append([k, round(f.geometry().area() / 10000, 2), round(sem_ha.get(k, 0), 2),
                            titik[k], round(n[k], 3), round(c[k], 2)])
    return hasil_baris


def selesai(exception, baris=None):
    """Program utama. Aman menyentuh proyek, tabel, dan tampilan."""
    if exception is not None:
        QgsMessageLog.logMessage("Ringkasan gagal: %s" % exception, "Ringkasan KPH", Qgis.Critical)
        return
    if baris is None:
        QgsMessageLog.logMessage("Ringkasan dibatalkan", "Ringkasan KPH", Qgis.Warning)
        return
    with open(HASIL + "/ringkasan_latar.csv", "w", newline="", encoding="utf-8") as f:
        w = csv.writer(f)
        w.writerow(["NAMA_KPH", "LUAS_HA", "SEMPADAN_HA", "HOTSPOT", "NDVI_RATA", "CHM_RATA_M"])
        w.writerows(baris)
    for nama in ("NDVI_KPH", "CHM"):
        QgsProject.instance().addMapLayer(QgsRasterLayer(HASIL + "/" + nama + ".tif", nama, "gdal"))
    for b in baris:
        print("%-10s luas %5.2f ha | sempadan %4.2f ha | titik panas %d | NDVI %.3f | CHM %.2f m"
              % tuple(b))
    QgsMessageLog.logMessage("Ringkasan selesai: %d KPH" % len(baris), "Ringkasan KPH", Qgis.Success)


tugas = QgsTask.fromFunction("Ringkasan per KPH", hitung_ringkasan, on_finished=selesai,
                             data=DATA, hasil=HASIL)
QgsApplication.taskManager().addTask(tugas)
print("Tugas berjalan di latar belakang; QGIS tetap bisa dipakai.")

Tugas berjalan di latar belakang; QGIS tetap bisa dipakai.
KPH Alpha  luas  4.00 ha | sempadan 1.18 ha | titik panas 2 | NDVI 0.604 | CHM 2.60 m
KPH Beta   luas  4.00 ha | sempadan 1.39 ha | titik panas 5 | NDVI 0.431 | CHM 2.69 m
KPH Gamma  luas  8.00 ha | sempadan 1.50 ha | titik panas 5 | NDVI 0.329 | CHM 2.56 m

# [SKRIP 7.3: Uji kesiapan otomatis: sepuluh keterampilan yang dipakai Masterclass Seri 1 (M1)]
# Penulis: Badar Mubarok Yogaswara
# Tiap pemeriksaan mencetak LULUS atau GAGAL beserta alasannya.
# Skor akhir muncul setelah tugas latar belakang (pemeriksaan 10) selesai.
import os
import re
import time
import processing
from qgis.core import (QgsApplication, QgsTask, QgsProject, QgsVectorLayer, QgsRasterLayer,
                       QgsExpression, QgsExpressionContext, QgsExpressionContextUtils,
                       QgsProcessingException)

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

proyek = QgsProject.instance()
proyek.removeAllMapLayers()
catatan = []                                     # daftar (nama, lulus?, alasan)
NDVI = '("Citra_KPH@4" - "Citra_KPH@1") / ("Citra_KPH@4" + "Citra_KPH@1")'
CACAH = ("aggregate('Hotspot_Nasional', 'count', \"ID_Titik\", "
         "intersects($geometry, geometry(@parent)))")


def periksa(nama, fungsi):
    try:
        alasan = fungsi()
        catatan.append((nama, True, alasan or ""))
    except Exception as e:           # pemeriksaan yang gagal dicatat; skrip tidak berhenti
        catatan.append((nama, False, "%s: %s" % (type(e).__name__, e)))
    n, ok, alasan = catatan[-1]
    print("%-46s %s  %s" % (n, "LULUS" if ok else "GAGAL", alasan))


def tegaskan(syarat, pesan):
    if not syarat:
        raise AssertionError(pesan)


batas = QgsVectorLayer(DATA + "/Batas_KPH.gpkg", "Batas_KPH", "ogr")
hotspot = QgsVectorLayer(DATA + "/Hotspot_Nasional.gpkg", "Hotspot_Nasional", "ogr")
citra = QgsRasterLayer(DATA + "/Citra_KPH.tif", "Citra_KPH", "gdal")


def konteks(fitur):
    c = QgsExpressionContext()
    c.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(batas))
    c.setFeature(fitur)
    return c


def k1():
    for lyr in (batas, hotspot, citra):
        tegaskan(lyr.isValid(), lyr.name() + " tidak sah")
        proyek.addMapLayer(lyr)
    tegaskan(len(proyek.mapLayers()) == 3, "jumlah layer tidak 3")


def k2():
    tegaskan(proyek.mapLayersByName("Batas_KPH")[0].featureCount() == 3, "Batas_KPH bukan 3 fitur")
    r = proyek.mapLayersByName("Citra_KPH")[0]
    tegaskan(r.bandCount() == 4 and r.rasterUnitsPerPixelX() == 2.0, "raster tidak sesuai")
    tegaskan(proyek.mapLayersByName("Tidak_Ada") == [], "nama tak ada harus memberi daftar kosong")
    return "band=%d, piksel=%.0f m" % (r.bandCount(), r.rasterUnitsPerPixelX())


def k3():
    v = processing.run("native:buffer", {"INPUT": batas, "DISTANCE": 10,
                                         "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    tegaskan(isinstance(v, QgsVectorLayer) and v.featureCount() == 3, "bukan layer vektor")
    r = processing.run("native:rastercalc", {"EXPRESSION": '"Citra_KPH@4" - "Citra_KPH@1"',
                                             "LAYERS": [citra], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    tegaskan(isinstance(r, str) and os.path.exists(r), "keluaran raster bukan jalur berkas")
    return "vektor -> layer, raster -> jalur teks"


def k4():
    jalur = processing.run("native:rastercalc", {"EXPRESSION": NDVI, "LAYERS": [citra],
                                                 "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    st = processing.run("native:rasterlayerstatistics",
                        {"INPUT": QgsRasterLayer(jalur, "ndvi", "gdal"), "BAND": 1})
    tegaskan(-1 <= st["MIN"] <= st["MAX"] <= 1, "NDVI di luar -1..1")
    return "NDVI %.2f sampai %.2f" % (st["MIN"], st["MAX"])


def k5():
    e = QgsExpression(CACAH)
    hasil = {f["NAMA_KPH"]: e.evaluate(konteks(f)) for f in batas.getFeatures()}
    benar = {f["NAMA_KPH"]: sum(1 for h in hotspot.getFeatures()
                                if h.geometry().intersects(f.geometry())) for f in batas.getFeatures()}
    tegaskan(len(hasil) == 3, "agregasi tidak menghasilkan 3 nilai: %s" % hasil)
    tegaskan(hasil == benar, "agregasi %s beda dari hitungan manual %s" % (hasil, benar))
    return str(hasil)


def k6():
    e = QgsExpression("'Titik api: ' || " + CACAH + " || ' titik'")
    beta = [x for x in batas.getFeatures() if x["NAMA_KPH"] == "KPH Beta"][0]
    teks = e.evaluate(konteks(beta))
    tegaskan(not e.hasEvalError() and teks == "Titik api: 5 titik", "teks label salah: %s" % teks)
    return teks


def k7():
    try:
        processing.run("native:buffer", {"INPUT": DATA + "/tidak_ada.gpkg", "DISTANCE": 5,
                                         "OUTPUT": "TEMPORARY_OUTPUT"})
    except QgsProcessingException:
        return "galat alat tertangkap"
    raise AssertionError("seharusnya gagal")


def k8():
    lyr = QgsVectorLayer(DATA + "/Kiriman_Lapangan/tim_e.gpkg", "rusak", "ogr")
    tegaskan(not lyr.isValid(), "berkas rusak dianggap sah")
    return "isValid() = False untuk berkas rusak"


def k9():
    berkas = HASIL + "/uji_kesiapan.qgz"
    tegaskan(proyek.write(berkas), "gagal menulis proyek")
    proyek.clear()
    tegaskan(len(proyek.mapLayers()) == 0, "proyek belum kosong")
    tegaskan(proyek.read(berkas) and len(proyek.mapLayers()) == 3, "proyek tidak terbaca utuh")
    return "tulis, kosongkan, baca ulang: 3 layer"


daftar = [("1. Muat vektor dan raster [Bab 1]", k1), ("2. Ambil layer dengan nama [Bab 1]", k2),
          ("3. processing.run sementara [Bab 2]", k3), ("4. Raster calculator NDVI [Bab 2]", k4),
          ("5. Agregasi dan @parent [Bab 5]", k5), ("6. Teks label dari ekspresi [Bab 5]", k6),
          ("7. Menangkap galat alat [Bab 4]", k7), ("8. Mendeteksi berkas rusak [Bab 4]", k8),
          ("9. Simpan dan buka proyek [Bab 1]", k9)]
for nama, fungsi in daftar:
    periksa(nama, fungsi)


# 10. Tugas latar belakang. Hasilnya tiba belakangan, jadi skor akhir dicetak di selesai().
def kerja(task):
    total = 0
    for i in range(1, 6):
        if task.isCanceled():
            return None
        time.sleep(0.05)
        total += i
        task.setProgress(i * 20)
    return total


def selesai(exception, nilai=None):
    ok = exception is None and nilai == 15
    catatan.append(("10. QgsTask sederhana [Bab 3]", ok, "hasil tugas = %s" % nilai))
    print("%-46s %s  hasil tugas = %s" % (catatan[-1][0], "LULUS" if ok else "GAGAL", nilai))
    lulus = sum(1 for _, o, _ in catatan if o)
    gagal = sorted({re.search(r"\[(Bab \d)\]", n).group(1) for n, o, _ in catatan if not o})
    akhir = "Siap masuk M1." if not gagal else "Ulangi: " + ", ".join(gagal) + "."
    print("SKOR: %d dari %d. %s" % (lulus, len(catatan), akhir))


tugas_uji = QgsTask.fromFunction("Uji tugas latar belakang", kerja, on_finished=selesai)
QgsApplication.taskManager().addTask(tugas_uji)

1. Muat vektor dan raster [Bab 1]              LULUS
2. Ambil layer dengan nama [Bab 1]             LULUS  band=4, piksel=2 m
3. processing.run sementara [Bab 2]            LULUS  vektor -> layer, raster -> jalur teks
4. Raster calculator NDVI [Bab 2]              LULUS  NDVI -0.02 sampai 0.84
5. Agregasi dan @parent [Bab 5]                LULUS  {'KPH Alpha': 2, 'KPH Beta': 5, 'KPH Gamma': 5}
6. Teks label dari ekspresi [Bab 5]            LULUS  Titik api: 5 titik
7. Menangkap galat alat [Bab 4]                LULUS  galat alat tertangkap
8. Mendeteksi berkas rusak [Bab 4]             LULUS  isValid() = False untuk berkas rusak
9. Simpan dan buka proyek [Bab 1]              LULUS  tulis, kosongkan, baca ulang: 3 layer
10. QgsTask sederhana [Bab 3]                  LULUS  hasil tugas = 15
SKOR: 10 dari 10. Siap masuk M1.

# [SKRIP 7.1-Pro: Proyek mini terpadu dengan ArcPy: ringkasan per KPH]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh Spatial Analyst. Data Shapefile dipakai agar sejajar dengan Bagian QGIS.
import csv
import arcpy
from arcpy.sa import Raster, Float, ZonalStatisticsAsTable

DATA = r"D:\Latihan\paket-i4"
HASIL = r"D:\Latihan\hasil_i4"
arcpy.env.workspace = DATA + r"\shp"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

try:
    # I1: sempadan 30 m dan titik panas
    arcpy.analysis.Buffer("Sungai.shp", r"memory\sem", "30 Meters", dissolve_option="ALL")
    arcpy.analysis.Intersect(["Batas_KPH.shp", r"memory\sem"], r"memory\iris")
    arcpy.analysis.SpatialJoin("Batas_KPH.shp", "Hotspot_Nasional.shp", r"memory\ringkas",
                               "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="INTERSECT")
    # I2: NDVI dan CHM, lalu rata-rata per KPH
    nir = Float(Raster(DATA + r"\Citra_KPH.tif\Band_4"))
    merah = Float(Raster(DATA + r"\Citra_KPH.tif\Band_1"))
    ndvi = (nir - merah) / (nir + merah)
    chm = Raster(DATA + r"\DSM_Drone.tif") - Raster(DATA + r"\DTM_Drone.tif")
    ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", ndvi, r"memory\z_ndvi", "DATA", "MEAN")
    ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", chm, r"memory\z_chm", "DATA", "MEAN")

    def tabel(jalur, kolom):
        with arcpy.da.SearchCursor(jalur, ["NAMA_KPH", kolom]) as k:
            return {n: v for n, v in k}

    sempadan = {}
    with arcpy.da.SearchCursor(r"memory\iris", ["NAMA_KPH", "SHAPE@AREA"]) as k:
        for n, luas in k:
            sempadan[n] = sempadan.get(n, 0) + luas / 10000
    n_rata, c_rata = tabel(r"memory\z_ndvi", "MEAN"), tabel(r"memory\z_chm", "MEAN")
    with open(HASIL + r"\ringkasan_kph_pro.csv", "w", newline="") as f:
        w = csv.writer(f)
        w.writerow(["NAMA_KPH", "LUAS_HA", "SEMPADAN_HA", "HOTSPOT", "NDVI_RATA", "CHM_RATA_M"])
        with arcpy.da.SearchCursor(r"memory\ringkas", ["NAMA_KPH", "SHAPE@AREA", "Join_Count"]) as k:
            for n, luas, jml in k:
                w.writerow([n, round(luas / 10000, 2), round(sempadan.get(n, 0), 2), jml,
                            round(n_rata[n], 3), round(c_rata[n], 2)])
    print("Ringkasan tersimpan.")
except arcpy.ExecuteError:
    print("Alat gagal:", arcpy.GetMessages(2))
finally:
    arcpy.CheckInExtension("Spatial")

# -*- coding: utf-8 -*-
# [SKRIP 7.1-ArcMap: Proyek mini terpadu dengan ArcPy: ringkasan per KPH (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh Spatial Analyst.
import csv
import arcpy
from arcpy.sa import Raster, Float, ZonalStatisticsAsTable

DATA = r"D:\Latihan\paket-i4"
HASIL = r"D:\Latihan\hasil_i4"
arcpy.env.workspace = DATA + r"\shp"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

try:
    # I1: sempadan 30 m dan titik panas
    arcpy.Buffer_analysis("Sungai.shp", r"in_memory\sem", "30 Meters", dissolve_option="ALL")
    arcpy.Intersect_analysis(["Batas_KPH.shp", r"in_memory\sem"], r"in_memory\iris")
    arcpy.SpatialJoin_analysis("Batas_KPH.shp", "Hotspot_Nasional.shp", r"in_memory\ringkas",
                               "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="INTERSECT")
    # I2: NDVI dan CHM, lalu rata-rata per KPH
    nir = Float(Raster(DATA + r"\Citra_KPH.tif\Band_4"))
    merah = Float(Raster(DATA + r"\Citra_KPH.tif\Band_1"))
    ndvi = (nir - merah) / (nir + merah)
    chm = Raster(DATA + r"\DSM_Drone.tif") - Raster(DATA + r"\DTM_Drone.tif")
    ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", ndvi, r"in_memory\z_ndvi", "DATA", "MEAN")
    ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", chm, r"in_memory\z_chm", "DATA", "MEAN")

    def tabel(jalur, kolom):
        with arcpy.da.SearchCursor(jalur, ["NAMA_KPH", kolom]) as k:
            return dict((n, v) for n, v in k)

    sempadan = {}
    with arcpy.da.SearchCursor(r"in_memory\iris", ["NAMA_KPH", "SHAPE@AREA"]) as k:
        for n, luas in k:
            sempadan[n] = sempadan.get(n, 0) + luas / 10000.0
    n_rata, c_rata = tabel(r"in_memory\z_ndvi", "MEAN"), tabel(r"in_memory\z_chm", "MEAN")
    with open(HASIL + r"\ringkasan_kph_arcmap.csv", "wb") as f:      # Python 2: mode "wb" untuk csv
        w = csv.writer(f)
        w.writerow(["NAMA_KPH", "LUAS_HA", "SEMPADAN_HA", "HOTSPOT", "NDVI_RATA", "CHM_RATA_M"])
        with arcpy.da.SearchCursor(r"in_memory\ringkas", ["NAMA_KPH", "SHAPE@AREA", "Join_Count"]) as k:
            for n, luas, jml in k:
                w.writerow([n, round(luas / 10000.0, 2), round(sempadan.get(n, 0), 2), jml,
                            round(n_rata[n], 3), round(c_rata[n], 2)])
    print("Ringkasan tersimpan.")
except arcpy.ExecuteError:
    print("Alat gagal: %s" % arcpy.GetMessages(2))
finally:
    arcpy.CheckInExtension("Spatial")