BAB 21: Memanggil Alat dari Skrip: Satu Perintah, Satu Analisis, Satu Resep untuk Banyak Berkas

Studi kasus: "Sempadan sungai, tiga puluh ubin, dan sore yang pendek" #

Kepala Seksi minta dua hal sebelum akhir pekan. Pertama, luas jalur sempadan sungai 30 meter di tiap KPH, lengkap dengan jumlah titik panasnya. Kedua, NDVI untuk semua ubin foto drone yang baru masuk. Mengklik alat satu per satu akan menghabiskan sore Anda. Dengan skrip, Anda memanggil alat Processing lewat tulisan, merangkainya, mengulanginya untuk banyak berkas, dan membiarkan QGIS tetap hidup selagi bekerja.

Konsep: memanggil alat dengan formulir pesanan #

Di QGIS, Anda memakai alat lewat jendela dengan kolom isian. Di skrip, kolom isian itu ditulis sebagai kamus parameter. Bayangkan formulir pesanan di loket: Anda menyebut nama layanan, mengisi kolomnya, lalu petugas mengembalikan lembar jawaban. Lembar jawaban itu adalah kamus hasil.

Nama alat di loket itu disebut ID alat, misalnya native:buffer. Bagian sebelum titik dua adalah nama kelompok alat (native berarti alat bawaan QGIS), dan bagian sesudahnya nama alatnya. Kamus sudah dibahas di Bab 5, jadi yang baru di sini hanyalah isinya.

Lima istilah baru Anda:

• Alat Processing: satu fungsi analisis di QGIS, seperti Buffer, yang juga ada di Processing Toolbox.
• ID alat: nama unik alat untuk skrip, misalnya native:buffer.
• Kamus parameter: isian alat, berupa pasangan nama dan nilai.
• Kamus hasil: jawaban alat, diambil dengan nama keluarannya, misalnya ["OUTPUT"].
• Keluaran sementara (TEMPORARY_OUTPUT): hasil perantara yang tidak ditulis ke berkas pilihan Anda.

Dua cara menjalankan, dan satu tambahan untuk Processing #

Bab 14 menjelaskan dua cara menjalankan skrip PyQGIS: di dalam QGIS, atau mandiri dengan kerangka QgsApplication. Semua skrip bab ini berjalan dengan cara pertama tanpa tambahan, sebab Konsol Python QGIS sudah memuat qgis.core dan processing.

Untuk cara mandiri, ada tambahan khusus. Python biasa tidak menemukan perkakas processing begitu saja. Pada uji, tulisan import processing di skrip mandiri menghasilkan ModuleNotFoundError: No module named 'processing'. Kerangka berikut sudah diuji pada QGIS 4.0.2.

import sys
from qgis.core import QgsApplication

QgsApplication.setPrefixPath(r"C:/Program Files/QGIS 4.0.2/apps/qgis", True)
app = QgsApplication([], False)
app.initQgis()

sys.path.append(r"C:/Program Files/QGIS 4.0.2/apps/qgis/python/plugins")
from processing.core.Processing import Processing
Processing.initialize()
from qgis.analysis import QgsNativeAlgorithms
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# ... tempel isi skrip Anda di sini ...

app.exitQgis()

Kerangka ini hanya untuk cara mandiri. Bila bagian Processing.initialize() dan addProvider dibuang, skrip tetap bisa mengimpor processing, tetapi alatnya tidak dikenal. Pada uji, hasilnya QgsProcessingException: Error: Algorithm native:buffer not found.

Skrip 21.1: Memanggil satu alat, Buffer #

Tujuan. Mencari ID alat Buffer, mengintip nama-nama isiannya, menjalankannya pada layer sungai dengan jarak 30 meter, lalu membaca hasilnya.

Bahan. Paket data paket-i4, berkas Sungai.gpkg (dua garis sungai). Berkas b2_01_jalankan_alat.py dijalankan dari Editor Skrip QGIS, atau dengan kerangka mandiri di atas. Ganti DATA di baris 6 dengan lokasi paket-i4 Anda.

Skrip ini bekerja dalam enam tahap. Gambar berikut menunjukkan urutannya.

Skrip utuh:

# [SKRIP 2.1: Menjalankan satu alat Processing dari Python (Buffer)]
# Penulis: Badar Mubarok Yogaswara
import processing
from qgis.core import QgsProject, QgsVectorLayer, QgsApplication

DATA = "D:/Latihan/paket-i4"
HASIL = "D:/Latihan/hasil_i4"

sungai = QgsVectorLayer(DATA + "/Sungai.gpkg", "Sungai", "ogr")

# 0. Mencari ID alat dengan kata kunci
for a in sorted(QgsApplication.processingRegistry().algorithms(), key=lambda x: x.id()):
    if a.id().startswith("native:") and "buffer" in a.id():
        print("ID:", a.id(), "-", a.displayName())

# 1. Mengintip "formulir" alat: nama parameter dan artinya
alat = QgsApplication.processingRegistry().algorithmById("native:buffer")
print("Alat:", alat.displayName())
for p in alat.parameterDefinitions():
    print("  ", p.name(), "=", p.description())

# 2. Menjalankan alat. Parameter berupa kamus (dictionary): nama parameter -> nilai.
hasil = processing.run("native:buffer", {
    "INPUT": sungai,
    "DISTANCE": 30,                 # satuan meter karena CRS-nya UTM
    "SEGMENTS": 8,
    "END_CAP_STYLE": 0,             # 0 = bulat
    "JOIN_STYLE": 0,                # 0 = bulat
    "MITER_LIMIT": 2,
    "DISSOLVE": True,               # gabungkan semua penyangga jadi satu bentuk
    "OUTPUT": "TEMPORARY_OUTPUT",   # keluaran sementara, tidak menulis berkas
})

# 3. Hasilnya juga kamus. Kuncinya nama keluaran. Untuk vektor sementara, isinya sudah berupa layer.
sempadan = hasil["OUTPUT"]
print("Jenis keluaran:", type(sempadan).__name__)
print("Fitur:", sempadan.featureCount())
luas = sum(f.geometry().area() for f in sempadan.getFeatures()) / 10000
print("Luas sempadan sungai:", round(luas, 2), "ha")

# 4. Beri nama yang jelas lalu tampilkan di proyek
sempadan.setName("Sempadan_30m")
QgsProject.instance().addMapLayer(sempadan)
print("Layer di proyek:", [l.name() for l in QgsProject.instance().mapLayers().values()])

Penjelasan baris demi baris. Dua blok. Blok pertama mengenali alatnya, blok kedua menjalankannya.

Blok A: mengenali alat (baris 3-20).

Blok B: menjalankan dan membaca hasil (baris 22-44).

Angka pilihan untuk END_CAP_STYLE dan JOIN_STYLE bisa dilihat dengan processing.algorithmHelp("native:buffer"). Pada uji, bantuan itu mencatat END_CAP_STYLE: 0 Round, 1 Flat, 2 Square. Untuk JOIN_STYLE: 0 Round, 1 Miter, 2 Bevel.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

ID: native:buffer - Buffer
ID: native:bufferbym - Variable width buffer (by M value)
ID: native:multiringconstantbuffer - Multi-ring buffer (constant distance)
ID: native:singlesidedbuffer - Single sided buffer
ID: native:taperedbuffer - Tapered buffers
ID: native:wedgebuffers - Create wedge buffers
Alat: Buffer
   INPUT = Input layer
   DISTANCE = Distance
   SEGMENTS = Segments
   END_CAP_STYLE = End cap style
   JOIN_STYLE = Join style
   MITER_LIMIT = Miter limit
   DISSOLVE = Dissolve result
   SEPARATE_DISJOINT = Keep disjoint results separate
   OUTPUT = Buffered
Jenis keluaran: QgsVectorLayer
Fitur: 1
Luas sempadan sungai: 4.48 ha
Layer di proyek: ['Sempadan_30m']

Ada enam alat yang ID-nya memuat kata "buffer". Alat Buffer punya sembilan isian, dan skrip hanya mengisi delapan. Isian SEPARATE_DISJOINT tidak diisi, jadi alat memakai nilai bawaannya. Hasil: satu fitur seluas 4,48 hektare.

Bagian yang boleh diubah.

• Baris 6, lokasi data. Ganti sesuai komputer Anda.
• Baris 13, kata buffer. Ganti dengan kata lain, misalnya clip atau raster, untuk mencari alat lain.
• Baris 17, ID alat. Ganti dengan ID lain, lalu baca isiannya di blok yang mengintip formulir.
• Baris 25, jarak. Ganti angka 30. Satuan mengikuti CRS.
• Baris 30, DISSOLVE. False menghasilkan satu fitur per sungai.
• Baris 42, nama layer. Ganti sesuai kebutuhan.

Jebakan. Semua kesalahan di bawah sudah dicoba pada QGIS 4.0.2.

• Nama isian salah ketik tidak menimbulkan galat. "DISTANC": 30 menggantikan "DISTANCE": 30. Skrip jalan mulus, tetapi alat memakai nilai bawaannya sendiri. Luas menjadi 1.44 ha, bukan 4.48 ha. Tidak ada pesan apa pun. Salin nama isian dari hasil blok yang mengintip formulir.
• Layer masukan gagal dibuka tidak menimbulkan galat. Bila Sungai.gpkg diganti Sungay.gpkg, hasilnya Fitur: 0 dan Luas sempadan sungai: 0.0 ha. Alat bekerja pada layer kosong. Periksa sungai.isValid() lebih dulu (Bab 14).
• ID alat salah ketik menimbulkan galat. "native:buffr" menghasilkan QgsProcessingException: Error: Algorithm native:buffr not found.
• Keluaran ke berkas bukan layer. Bila "OUTPUT" diisi jalur berkas, hasil["OUTPUT"] berisi teks jalur itu. Baris 36 mencetak Jenis keluaran: str, lalu baris 37 berhenti dengan AttributeError: 'str' object has no attribute 'featureCount'.
• Satuan jarak mengikuti CRS. Sungai.gpkg diubah ke EPSG:4326 (derajat), lalu skrip dijalankan. Hasilnya Luas sempadan sungai: 0.28 ha, angka yang tidak bermakna, karena 30 dibaca sebagai 30 derajat. Pastikan CRS berupa meter sebelum memakai DISTANCE.

Latihan variasi.

1. Jalankan lagi dengan DISSOLVE bernilai False. Berapa fitur dan berapa luasnya?
2. Mengapa luasnya berbeda dari 4,48 hektare?

Contoh jawaban (sudah diuji). Untuk nomor 1, ubah baris 30 menjadi "DISSOLVE": False,. Layar menampilkan Fitur: 2 dan Luas sempadan sungai: 4.82 ha. Untuk nomor 2, kemungkinan besar bagian yang tumpang tindih di pertemuan dua sungai terhitung dua kali bila tidak digabung. Selisihnya 0,34 hektare [kemungkinan, belum diukur terpisah].

Skrip 21.2: Merangkai alat vektor #

Keluaran satu alat bisa langsung menjadi masukan alat berikutnya, seperti ban berjalan di pabrik. Skrip ini merangkai tiga alat, lalu menggabungkan hasilnya dalam satu tabel kecil.

Tujuan. Menghitung luas sempadan sungai 30 meter di tiap KPH, persennya terhadap luas KPH, dan jumlah titik panas per KPH.

Bahan. Paket data paket-i4: Batas_KPH.gpkg, Sungai.gpkg, dan Hotspot_Nasional.gpkg. Berkas b2_02_rantai_vektor.py. Ganti DATA di baris 6.

Skrip utuh:

# [SKRIP 2.2: Rantai alat vektor: sempadan sungai dan titik panas per KPH]
# Penulis: Badar Mubarok Yogaswara
import processing
from qgis.core import QgsProject, QgsVectorLayer

DATA = "D:/Latihan/paket-i4"
HASIL = "D:/Latihan/hasil_i4"

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")

# Langkah 1: penyangga (buffer) 30 m di kiri-kanan sungai
sempadan = processing.run("native:buffer", {
    "INPUT": sungai, "DISTANCE": 30, "SEGMENTS": 8, "END_CAP_STYLE": 0, "JOIN_STYLE": 0,
    "MITER_LIMIT": 2, "DISSOLVE": True, "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

# Langkah 2: irisan (intersection). Hasil langkah 1 langsung menjadi masukan langkah 2.
irisan = processing.run("native:intersection", {
    "INPUT": batas, "OVERLAY": sempadan, "INPUT_FIELDS": ["NAMA_KPH"],
    "OVERLAY_FIELDS": [], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

# Langkah 3: hitung titik panas di tiap KPH
dengan_titik = processing.run("native:countpointsinpolygon", {
    "POLYGONS": batas, "POINTS": hotspot, "FIELD": "JML_HOTSPOT", "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

# Langkah 4: gabungkan angkanya dalam satu tabel kecil
luas_sempadan = {}
for f in irisan.getFeatures():
    nama = f["NAMA_KPH"]
    luas_sempadan[nama] = luas_sempadan.get(nama, 0) + f.geometry().area() / 10000

print("KPH         luas(ha)  sempadan(ha)  persen  titik_panas")
for f in dengan_titik.getFeatures():
    nama = f["NAMA_KPH"]
    luas = f.geometry().area() / 10000
    sem = luas_sempadan.get(nama, 0)
    print("%-11s %8.2f %13.2f %7.1f %12d" % (nama, luas, sem, 100 * sem / luas, f["JML_HOTSPOT"]))

Skrip ini mirip Skrip 21.1, jadi hanya bagian yang baru yang dijelaskan. Baris import, DATA, dan pemuatan layer (baris 3-11) sudah dibahas sebelumnya.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

KPH         luas(ha)  sempadan(ha)  persen  titik_panas
KPH Alpha       4.00          1.18    29.4            2
KPH Beta        4.00          1.39    34.6            5
KPH Gamma       8.00          1.50    18.7            5

Dari tabel ini Anda tahu: sempadan paling luas di KPH Gamma (1,50 hektare), tetapi persen terbesar ada di KPH Beta (34,6 persen).

Bagian yang boleh diubah.

• Baris 15, "DISTANCE": 30. Lebar sempadan.
• Baris 25, "FIELD": "JML_HOTSPOT". Nama kolom penghitung. Bila Anda mengubahnya, ubah juga f["JML_HOTSPOT"] di baris 38.
• Baris 20, INPUT_FIELDS. Tambahkan nama kolom lain yang ingin dibawa.
• Baris 33 dan 38, format tabel. Lebar kolom boleh disesuaikan, tetapi judul dan baris data harus selaras.

Jebakan. Semua dicoba pada QGIS 4.0.2.

• Nama kolom tidak sama. f["JML_HOTSPOT"] diganti f["JML_TITIK"]: judul tabel tercetak, lalu skrip berhenti dengan KeyError: 'JML_TITIK'.
• Layer sungai gagal dibuka. Dengan Sungay.gpkg, tidak ada galat. Kolom sempadan berisi 0.00 dan persen 0.0 untuk ketiga KPH. Angka tampak wajar, padahal kosong.
• Jarak nol. "DISTANCE": 0 menghentikan alat Intersection dengan QgsProcessingException: Could not write feature into OUTPUT: Could not add feature with geometry type GeometryCollection to layer of type MultiPolygon. Jarak nol menghasilkan bentuk yang tidak cocok untuk alat berikutnya.
• Pembagian dengan nol. Baris 38 membagi dengan luas. Data latihan tidak punya KPH berluas nol, jadi galat ZeroDivisionError tidak muncul di sini (Bab 10 membahas penangannya).

Latihan variasi. Ubah lebar sempadan menjadi 50 meter. Bagaimana tabelnya berubah?

Contoh jawaban (sudah diuji). Ganti baris 15 menjadi "INPUT": sungai, "DISTANCE": 50, "SEGMENTS": 8, "END_CAP_STYLE": 0, "JOIN_STYLE": 0,. Layar menampilkan:

KPH         luas(ha)  sempadan(ha)  persen  titik_panas
KPH Alpha       4.00          1.69    42.2            2
KPH Beta        4.00          2.32    58.0            5
KPH Gamma       8.00          2.54    31.8            5

Skrip 21.3: Alat raster, NDVI dan tinggi tajuk #

Raster adalah kisi petak (piksel), jadi alatnya sedikit berbeda dari vektor. Skrip ini memakai dua alat raster dan satu alat statistik. Hasilnya dua angka per KPH.

Dua istilah baru. NDVI adalah indeks kehijauan tanaman dari citra, dihitung (NIR - Merah) / (NIR + Merah). Nilai tinggi berarti vegetasi rapat. CHM adalah tinggi tajuk pohon, dihitung dari raster tinggi permukaan (DSM) dikurangi raster tinggi tanah (DTM). Statistik zonal meringkas nilai raster di dalam tiap poligon.

Tujuan. Menghitung NDVI dan tinggi tajuk, lalu mencetak rata-ratanya untuk tiap KPH.

Bahan. Paket data paket-i4: Batas_KPH.gpkg, Citra_KPH.tif (empat band, band 1 Merah dan band 4 NIR), DSM_Drone.tif, dan DTM_Drone.tif. Berkas b2_03_raster_alat.py. Ganti DATA di baris 6.

Hal yang perlu Anda tahu dulu: keluaran raster sementara berbeda dari vektor sementara.

Skrip utuh:

# [SKRIP 2.3: Alat raster dari Python: NDVI, tinggi tajuk (CHM), dan statistik per KPH]
# Penulis: Badar Mubarok Yogaswara
import processing
from qgis.core import QgsProject, QgsVectorLayer, QgsRasterLayer

DATA = "D:/Latihan/paket-i4"
HASIL = "D:/Latihan/hasil_i4"

batas = QgsVectorLayer(DATA + "/Batas_KPH.gpkg", "Batas_KPH", "ogr")
citra = QgsRasterLayer(DATA + "/Citra_KPH.tif", "Citra_KPH", "gdal")   # band 1=Merah, 4=NIR
dsm = QgsRasterLayer(DATA + "/DSM_Drone.tif", "DSM", "gdal")
dtm = QgsRasterLayer(DATA + "/DTM_Drone.tif", "DTM", "gdal")

# 1. NDVI = (NIR - Merah) / (NIR + Merah). Dalam ekspresi, band ditulis "NamaLayer@nomorBand".
jalur_ndvi = processing.run("native:rastercalc", {
    "EXPRESSION": '("Citra_KPH@4" - "Citra_KPH@1") / ("Citra_KPH@4" + "Citra_KPH@1")',
    "LAYERS": [citra], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
print("Keluaran raster sementara berupa:", type(jalur_ndvi).__name__)    # teks jalur, bukan layer
ndvi = QgsRasterLayer(jalur_ndvi, "NDVI_KPH", "gdal")                    # bungkus jadi layer
print("NDVI sah?", ndvi.isValid())

# 2. CHM = DSM - DTM
jalur_chm = processing.run("native:rastercalc", {
    "EXPRESSION": '"DSM@1" - "DTM@1"', "LAYERS": [dsm, dtm], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
chm = QgsRasterLayer(jalur_chm, "CHM", "gdal")

# 3. Statistik zonal: rata-rata nilai raster di dalam tiap poligon KPH
def rata_rata(raster, awalan):
    return processing.run("native:zonalstatisticsfb", {
        "INPUT": batas, "INPUT_RASTER": raster, "RASTER_BAND": 1,
        "COLUMN_PREFIX": awalan, "STATISTICS": [2],      # 2 = Mean
        "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

zn = rata_rata(ndvi, "ndvi_")
zc = rata_rata(chm, "chm_")
chm_per_kph = {f["NAMA_KPH"]: f["chm_mean"] for f in zc.getFeatures()}
print("KPH         NDVI rata-rata   CHM rata-rata (m)")
for f in zn.getFeatures():
    print("%-11s %12.3f %16.2f" % (f["NAMA_KPH"], f["ndvi_mean"], chm_per_kph[f["NAMA_KPH"]]))

QgsProject.instance().addMapLayer(ndvi)
QgsProject.instance().addMapLayer(chm)

Penjelasan baris demi baris. Hanya yang baru yang dijelaskan.

Pada uji, bantuan alat native:zonalstatisticsfb mencatat kode statistik: 0 Count, 1 Sum, 2 Mean, 3 Median, 4 St dev, 5 Minimum, 6 Maximum, dan seterusnya.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Keluaran raster sementara berupa: str
NDVI sah? True
KPH         NDVI rata-rata   CHM rata-rata (m)
KPH Alpha          0.604             2.60
KPH Beta           0.431             2.69
KPH Gamma          0.329             2.56

Dari hasil ini Anda tahu: vegetasi paling hijau ada di KPH Alpha (NDVI 0,604), sedangkan tajuk paling tinggi rata-rata di KPH Beta (2,69 meter).

Bagian yang boleh diubah.

• Baris 16, rumus NDVI. Nomor band harus sesuai urutan band citra Anda. Citra ini: band 1 Merah, band 4 NIR.
• Baris 31, "STATISTICS": [2]. Ganti kode statistik. Perhatikan: resep rata_rata dipakai untuk NDVI dan CHM sekaligus, jadi perubahannya berlaku untuk keduanya.
• Baris 34-36, awalan kolom. Bila Anda mengubah awalan, ubah juga nama kolom yang dibaca.
• Baris 24, rumus CHM. Ganti raster atau rumusnya.

Jebakan. Semua dicoba pada QGIS 4.0.2.

• Band tertukar, tanpa galat. Rumus NDVI ditulis terbalik: "Citra_KPH@1" di tempat "Citra_KPH@4" dan sebaliknya. Skrip jalan mulus, tetapi tabel menampilkan -0.604, -0.431, dan -0.329. Tanda minus adalah petunjuk bahwa band terbalik.
• Nama layer di rumus tidak cocok. "DSM@1" diganti "DSM_Drone@1", padahal layer bernama "DSM". Hasilnya QgsProcessingException: Error occurred while performing calculation.
• Nama kolom hasil salah. Menulis f["chm_rata"] menghasilkan KeyError: 'chm_rata'. Kolom hasil selalu awalan ditambah nama statistik (mean, max, dan seterusnya).
• Hasil raster adalah teks. Memanggil jalur_ndvi.width() gagal, karena str tidak punya metode itu. Bungkus dulu dengan QgsRasterLayer.

Latihan variasi. Ganti kolom CHM menjadi nilai tertinggi per KPH, bukan rata-rata.

Contoh jawaban (sudah diuji). Ubah tiga hal: pada baris 31 jadi "STATISTICS": [2, 6] (rata-rata dan tertinggi sekaligus), pada baris 36 jadi f["chm_max"], dan pada baris 37 ganti judul menjadi CHM tertinggi (m). Layar menampilkan:

KPH         NDVI rata-rata   CHM tertinggi (m)
KPH Alpha          0.604            13.94
KPH Beta           0.431            13.96
KPH Gamma          0.329            13.96

Konsep: batch, satu resep untuk banyak berkas #

Tiga skrip berikutnya membahas batch: pekerjaan yang sama diulang untuk banyak masukan. Bayangkan Anda memasak resep yang sama untuk tiga puluh porsi. Resepnya ditulis sekali, lalu perulangan for (Bab 6) memasaknya satu porsi demi satu porsi. Buku I4 Bab 3 menunjukkan cara bawaan QGIS (Execute as batch process). Cara itu cocok untuk beberapa berkas. Untuk puluhan berkas, skrip lebih ringkas.

Istilah baru bagian ini:

• Batch: pekerjaan sama, banyak masukan.
• Pola berkas (glob): cara mencari banyak berkas sekaligus, misalnya semua yang berakhiran .tif.
• Latar belakang: jalur kerja terpisah, supaya QGIS tidak membeku.

Skrip 21.4: Batch sederhana, NDVI untuk semua ubin #

Tujuan. Menghitung NDVI untuk semua ubin ortofoto di satu folder dengan satu perulangan for, menyimpan tiap hasil sebagai berkas, dan mencetak rata-ratanya.

Bahan. Paket data paket-i4, folder Ortofoto_Mentah (tiga ubin; pada data ini band 1 adalah NIR dan band 2 adalah Merah). Berkas b3_01_batch_ndvi.py. Ganti DATA di baris 8 dan HASIL di baris 9. Skrip membuat folder NDVI_ubin di dalam HASIL.

Skrip utuh:

# [SKRIP 3.1: Batch sederhana: NDVI untuk semua ubin ortofoto dalam satu folder]
# Penulis: Badar Mubarok Yogaswara
import glob
import os
import processing
from qgis.core import QgsRasterLayer

DATA = "D:/Latihan/paket-i4"
HASIL = "D:/Latihan/hasil_i4"
folder_masuk = DATA + "/Ortofoto_Mentah"
folder_keluar = HASIL + "/NDVI_ubin"
os.makedirs(folder_keluar, exist_ok=True)       # buat folder bila belum ada

# 1. Kumpulkan semua berkas .tif, diurutkan agar hasilnya selalu berurutan sama
daftar = sorted(glob.glob(folder_masuk + "/*.tif"))
print("Ditemukan", len(daftar), "ubin")

# 2. Ulangi pekerjaan yang sama untuk tiap berkas. Pada ubin ini band 1 = NIR dan band 2 = Merah.
berhasil = []
for jalur in daftar:
    nama = os.path.basename(jalur)                         # contoh: ubin_01.tif
    keluar = folder_keluar + "/NDVI_" + nama
    ubin = QgsRasterLayer(jalur, "ubin", "gdal")           # nama "ubin" dipakai di ekspresi
    processing.run("native:rastercalc", {
        "EXPRESSION": '("ubin@1" - "ubin@2") / ("ubin@1" + "ubin@2")',
        "LAYERS": [ubin], "OUTPUT": keluar})
    hasil = QgsRasterLayer(keluar, "hasil", "gdal")
    rata = processing.run("native:rasterlayerstatistics", {"INPUT": hasil, "BAND": 1})["MEAN"]
    print("  %s -> %s  | NDVI rata-rata %.3f" % (nama, os.path.basename(keluar), rata))
    berhasil.append(keluar)

print("Selesai:", len(berhasil), "dari", len(daftar), "berkas")

Penjelasan baris demi baris.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Ditemukan 3 ubin
  ubin_01.tif -> NDVI_ubin_01.tif  | NDVI rata-rata 0.471
  ubin_02.tif -> NDVI_ubin_02.tif  | NDVI rata-rata 0.471
  ubin_03.tif -> NDVI_ubin_03.tif  | NDVI rata-rata 0.467
Selesai: 3 dari 3 berkas

Di folder NDVI_ubin kini ada tiga berkas NDVI. QGIS dapat menambahkan berkas pendamping berakhiran .aux.xml. Berkas itu aman dibiarkan.

Bagian yang boleh diubah.

• Baris 10, folder masukan. Ganti /Ortofoto_Mentah dengan folder lain.
• Baris 15, pola "/*.tif". Ganti ke "/ubin_0*.tif" untuk memilih sebagian berkas.
• Baris 25, rumus. Urutan band harus sesuai data Anda. Perhatikan: di Citra_KPH.tif (Skrip 21.3), band 4 adalah NIR, sedangkan di ubin ini band 1.
• Baris 22, awalan "NDVI_". Ganti sesuai selera.

Jebakan. Semua dicoba pada QGIS 4.0.2.

• Folder salah, tanpa galat. Ortofoto_Mentah diganti Ortofoto_Mentahh. Layar menampilkan Ditemukan 0 ubin dan Selesai: 0 dari 0 berkas. Perulangan tidak berputar sama sekali. Baca baris pertama hasil.
• Satu berkas rusak menghentikan semuanya. Folder Ortofoto_Campur berisi empat ubin, dan ubin keempat rusak. Layar menampilkan Ditemukan 4 ubin, tiga baris hasil, lalu ERROR 1: Attempt to create -2147483648x-2147483648 dataset is illegal,sizes must be larger than zero. dan berakhir dengan QgsProcessingException: Error creating output file. Tiga ubin pertama sudah selesai, tetapi ringkasan akhir tidak pernah tercetak. Penanganannya (try dan except) dibahas di Bab 10.
• Band terbalik, tanpa galat. Rumus ditulis ("ubin@2" - "ubin@1") / .... Hasilnya -0.471, -0.471, dan -0.467. Tanda minus menandakan band terbalik.
• Band yang tidak ada. Menulis "ubin@4" pada ubin dua band menghasilkan QgsProcessingException: Error occurred while performing calculation.
• berhasil selalu sama panjang. Daftar berhasil selalu terisi setiap putaran, sehingga ringkasan "3 dari 3" belum berarti ada penyaringan. Ia baru berguna bila disertai penanganan galat.

Latihan variasi. Tambahkan satu baris di dalam perulangan yang mencetak ukuran piksel tiap hasil.

Contoh jawaban (sudah diuji). Sisipkan sebelum berhasil.append(keluar):

    print("     ukuran:", hasil.width(), "x", hasil.height(), "piksel")

Di bawah tiap baris NDVI, layar menampilkan ukuran: 100 x 100 piksel.

Skrip 21.5: Satu putaran untuk tiap fitur #

Perulangan tidak harus atas berkas. Skrip ini berputar atas fitur: untuk tiap KPH, ia memilih satu poligon lalu memotong raster tinggi tajuk (CHM) dengannya. Hasilnya satu raster per KPH.

Tujuan. Menghasilkan raster CHM terpisah untuk tiap KPH, lalu mencetak ukuran, rata-rata, dan nilai tertingginya.

Bahan. Paket data paket-i4: Batas_KPH.gpkg, DSM_Drone.tif, dan DTM_Drone.tif. Berkas b3_02_batch_per_kph.py. Ganti DATA di baris 7 dan HASIL di baris 8. Skrip membuat folder CHM_per_KPH.

Skrip utuh:

# [SKRIP 3.2: Batch per fitur: potong raster tinggi tajuk (CHM) untuk tiap KPH]
# Penulis: Badar Mubarok Yogaswara
import os
import processing
from qgis.core import QgsVectorLayer, QgsRasterLayer

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

batas = QgsVectorLayer(DATA + "/Batas_KPH.gpkg", "Batas_KPH", "ogr")
dsm = QgsRasterLayer(DATA + "/DSM_Drone.tif", "DSM", "gdal")
dtm = QgsRasterLayer(DATA + "/DTM_Drone.tif", "DTM", "gdal")

# 1. CHM untuk seluruh area (satu kali saja, di luar pengulangan)
chm_jalur = processing.run("native:rastercalc", {
    "EXPRESSION": '"DSM@1" - "DTM@1"', "LAYERS": [dsm, dtm], "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

# 2. Satu putaran untuk tiap KPH: pilih satu poligon, lalu potong CHM dengannya
for fitur in batas.getFeatures():
    nama = fitur["NAMA_KPH"]                                 # contoh: KPH Alpha
    aman = nama.replace(" ", "_")                            # untuk nama berkas: KPH_Alpha
    satu_kph = processing.run("native:extractbyexpression", {
        "INPUT": batas, "EXPRESSION": '"NAMA_KPH" = \'%s\'' % nama,
        "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
    keluar = folder_keluar + "/CHM_" + aman + ".tif"
    processing.run("gdal:cliprasterbymasklayer", {
        "INPUT": chm_jalur, "MASK": satu_kph, "CROP_TO_CUTLINE": True,
        "NODATA": -9999, "OUTPUT": keluar})
    hasil = QgsRasterLayer(keluar, aman, "gdal")
    st = processing.run("native:rasterlayerstatistics", {"INPUT": hasil, "BAND": 1})
    print("%-10s %4d x %4d piksel | rata-rata %.2f m | tertinggi %.2f m" % (
        aman, hasil.width(), hasil.height(), st["MEAN"], st["MAX"]))
print("Berkas di folder:", sorted(os.listdir(folder_keluar)))

Penjelasan baris demi baris. Hanya yang baru yang dijelaskan.

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

KPH_Alpha   400 x  400 piksel | rata-rata 2.60 m | tertinggi 13.94 m
KPH_Beta    400 x  400 piksel | rata-rata 2.69 m | tertinggi 13.96 m
KPH_Gamma   800 x  400 piksel | rata-rata 2.56 m | tertinggi 13.96 m
Berkas di folder: ['CHM_KPH_Alpha.tif', 'CHM_KPH_Alpha.tif.aux.xml', 'CHM_KPH_Beta.tif', 'CHM_KPH_Beta.tif.aux.xml', 'CHM_KPH_Gamma.tif']

Perhatikan daftar berkas di baris terakhir. Berkas .aux.xml dibuat oleh QGIS saat menghitung statistik, bukan oleh skrip secara sengaja. Isi daftar di komputer Anda bisa berbeda, misalnya CHM_KPH_Gamma.tif.aux.xml ikut muncul bila skrip dijalankan untuk kedua kalinya. KPH Gamma berukuran 800 kali 400 piksel karena areanya dua kali lebih luas.

Bagian yang boleh diubah.

• Baris 9-10, folder keluaran. Ganti nama folder.
• Baris 25, kolom NAMA_KPH. Ganti bila data Anda memakai kolom lain.
• Baris 29, CROP_TO_CUTLINE. True memotong ukuran hasil mengikuti KPH.
• Baris 30, NODATA. Nilai -9999 menandai piksel kosong. Nilai ini ikut dikenali saat statistik dihitung.
• Baris 17-18, rumus CHM. Ganti rumus atau raster, misalnya untuk NDVI.

Jebakan. Semua dicoba pada QGIS 4.0.2.

• Nama KPH yang memuat kutip tunggal. Pada uji, nama KPH D'Alpha disisipkan ke baris 25. Alat berhenti dengan Incomplete expression. You might not have finished the full expression. Kutip tunggal di tengah nama mengakhiri teks ekspresi lebih awal. Data latihan tidak memuat nama seperti itu.
• Nama kolom salah. fitur["NAMA_KPH"] diganti fitur["NAMA"]: KeyError: 'NAMA'.
• CROP_TO_CUTLINE bernilai False. Tidak ada galat. Ketiga hasil berukuran 800 x 800 piksel, karena ukuran tidak dipotong mengikuti KPH. Rata-rata Beta dan Gamma di layar tetap 2.69 dan 2.56.
• Berkas pendamping di daftar. Karena baris 35 mencetak semua isi folder, berkas .aux.xml ikut tampil.

Latihan variasi. Jalankan hanya untuk KPH Gamma. Bagian mana yang cukup ditambah?

Contoh jawaban (sudah diuji). Sisipkan dua baris setelah baris 22:

    if nama != "KPH Gamma":
        continue

continue berarti "lewati sisa putaran ini, lanjut ke fitur berikutnya". Layar menampilkan satu baris KPH_Gamma 800 x 400 piksel | rata-rata 2.56 m | tertinggi 13.96 m, disusul daftar isi folder (yang juga memuat berkas dari uji sebelumnya).

Konsep: latar belakang dan fungsi selesai #

Skrip 21.4 dan 21.5 membuat QGIS terasa macet selama bekerja. Bayangkan Anda pelayan restoran yang juga harus memasak. Selama Anda memasak, tamu tidak terlayani. Solusinya: serahkan masakan ke juru masak di dapur belakang, sementara Anda tetap melayani. QGIS menyediakan juru masak ini sebagai QgsTask, yaitu tugas latar belakang.

Pembagian kerjanya begini. Fungsi tugas berjalan di dapur belakang dan hanya boleh memasak (menghitung). Fungsi selesai berjalan di depan setelah masakan jadi, dan aman untuk menyentuh proyek, misalnya menambah layer. Fungsi selesai disebut juga callback, yaitu fungsi yang dipanggil kembali setelah tugas berakhir.

Tiga istilah baru:

• QgsTask: tugas latar belakang milik QGIS.
• Fungsi tugas: fungsi yang dikerjakan di latar belakang.
• Fungsi selesai: fungsi yang dipanggil di program utama setelah tugas berakhir.

Skrip 21.6: Batch di latar belakang dengan QgsTask #

Tujuan. Menjalankan batch NDVI di latar belakang, lengkap dengan bilah kemajuan dan kemungkinan pembatalan, lalu menambahkan hasilnya ke proyek.

Bahan. Paket data paket-i4, folder Ortofoto_Mentah. Berkas b3_03_batch_latar_belakang.py. Ganti DATA di baris 8 dan HASIL di baris 9. Skrip membuat folder NDVI_latar. Skrip ini paling pas dijalankan di Editor Skrip QGIS yang terbuka, karena tugas latar belakang butuh QGIS yang hidup.

Skrip utuh:

# [SKRIP 3.3: Batch di latar belakang dengan QgsTask supaya QGIS tidak macet]
# Penulis: Badar Mubarok Yogaswara
import glob
import os
import processing
from qgis.core import QgsApplication, QgsTask, QgsMessageLog, Qgis, QgsRasterLayer, QgsProject

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


def hitung_ndvi(task, folder_in, folder_out):
    """Berjalan di LATAR BELAKANG. Jangan menyentuh proyek atau antarmuka di sini."""
    daftar = sorted(glob.glob(folder_in + "/*.tif"))
    keluar = []
    for i, jalur in enumerate(daftar):
        if task.isCanceled():                       # pengguna menekan Cancel
            return None
        tujuan = folder_out + "/NDVI_" + os.path.basename(jalur)
        ubin = QgsRasterLayer(jalur, "ubin", "gdal")
        processing.run("native:rastercalc", {
            "EXPRESSION": '("ubin@1" - "ubin@2") / ("ubin@1" + "ubin@2")',
            "LAYERS": [ubin], "OUTPUT": tujuan})
        keluar.append(tujuan)
        task.setProgress(100 * (i + 1) / len(daftar))    # bilah kemajuan di pojok kanan bawah QGIS
    return keluar


def selesai(exception, hasil=None):
    """Berjalan di PROGRAM UTAMA setelah tugas berakhir. Aman untuk menyentuh proyek."""
    if exception is not None:
        QgsMessageLog.logMessage("Tugas gagal: %s" % exception, "Batch NDVI", Qgis.Critical)
    elif hasil is None:
        QgsMessageLog.logMessage("Tugas dibatalkan", "Batch NDVI", Qgis.Warning)
    else:
        for jalur in hasil:
            nama = os.path.basename(jalur)[:-4]
            QgsProject.instance().addMapLayer(QgsRasterLayer(jalur, nama, "gdal"))
        QgsMessageLog.logMessage("Selesai: %d berkas" % len(hasil), "Batch NDVI", Qgis.Success)
        print("Tugas selesai:", len(hasil), "berkas NDVI ditambahkan ke proyek")


tugas = QgsTask.fromFunction("Hitung NDVI semua ubin", hitung_ndvi, on_finished=selesai,
                             folder_in=folder_masuk, folder_out=folder_keluar)
QgsApplication.taskManager().addTask(tugas)
print("Tugas dikirim ke latar belakang. QGIS tetap bisa dipakai.")

Penjelasan baris demi baris. Skrip ini terdiri atas dua fungsi dan tiga baris pengiriman. Rumus NDVI di baris 25 sama dengan Skrip 21.4.

Blok A: fungsi tugas, yang bekerja di latar belakang (baris 15-29).

Blok B: fungsi selesai, yang bekerja di depan (baris 32-43).

Blok C: membuat dan mengirim tugas (baris 46-49).

Hasil. Layar menampilkan (berdasarkan uji pada QGIS 4.0.2):

Tugas dikirim ke latar belakang. QGIS tetap bisa dipakai.
Tugas selesai: 3 berkas NDVI ditambahkan ke proyek

Pesan pertama muncul langsung. Pesan kedua muncul beberapa saat kemudian, setelah tugas selesai. Di dalam QGIS, kemajuan tugas tampil di bilah status dan pesan lognya di panel Log Messages dengan tag Batch NDVI [CEK: tampilan visualnya, tidak dapat diuji tanpa antarmuka QGIS]. Tiga layer baru NDVI_ubin_01 sampai NDVI_ubin_03 masuk ke proyek.

Bagian yang boleh diubah.

• Baris 8-9, jalur. Sesuaikan dengan komputer Anda.
• Baris 17, daftar berkas. Ubah pola berkas atau batasi jumlahnya.
• Baris 25, rumus. Sesuaikan urutan band.
• Baris 35, 37, dan 42, tag log "Batch NDVI". Ganti dengan nama yang Anda suka.
• Baris 46, judul tugas. Judul ini tampil di daftar tugas QGIS.

Jebakan.

• Skrip tidak menunggu tugas. Pesan "Tugas dikirim" tampil sebelum NDVI selesai. Itu normal di dalam QGIS.
• Skrip mandiri harus menunggu sendiri. Pada uji, skrip ini dijalankan dengan kerangka mandiri, lalu dibiarkan habis begitu saja. Layar hanya menampilkan Tugas dikirim ..., kemudian proses berhenti mendadak (kode keluar 139, tanda proses mogok). Pesan "Tugas selesai" tidak pernah muncul. Perbaikannya: sebelum app.exitQgis(), tunggu tugas dengan perulangan berikut. Dengan perulangan ini, uji berjalan normal dan menampilkan kedua pesan.

import time

while QgsApplication.taskManager().countActiveTasks() > 0:
    app.processEvents()
    time.sleep(0.05)
for _ in range(20):
    app.processEvents()
    time.sleep(0.02)

• Fungsi tugas tidak boleh menyentuh proyek atau antarmuka. Komentar di baris 16 dan 33 dan Buku I4 mengingatkan bahwa menambah layer dari fungsi tugas dapat membuat QGIS macet atau tertutup sendiri. Itu sebabnya addMapLayer ada di fungsi selesai. Bagian ini tidak diuji [CEK].
• Berkas pendamping. Seperti Skrip 21.4, QGIS dapat membuat berkas .aux.xml di folder NDVI_latar.

Latihan variasi. Batasi tugas hanya untuk dua ubin pertama.

Contoh jawaban (sudah diuji). Ubah baris 17 menjadi:

    daftar = sorted(glob.glob(folder_in + "/*.tif"))[:2]

[:2] mengambil dua isi pertama dari daftar. Layar menampilkan Tugas selesai: 2 berkas NDVI ditambahkan ke proyek, dan folder NDVI_latar berisi NDVI_ubin_01.tif dan NDVI_ubin_02.tif.

Cek paham #

1. Apa isi hasil["OUTPUT"] bila OUTPUT berupa "TEMPORARY_OUTPUT" pada alat vektor? Bagaimana bila alatnya raster?
2. Mengapa CHM dihitung di luar perulangan pada Skrip 21.5?
3. Fungsi mana yang boleh menambahkan layer ke proyek pada Skrip 21.6, dan mengapa?

Jawaban.

1. Pada alat vektor, isinya layer (QgsVectorLayer) yang siap dipakai. Pada alat raster, isinya teks jalur berkas, yang harus dibungkus dengan QgsRasterLayer.
2. CHM seluruh area sama untuk semua KPH. Menghitungnya sekali cukup, dan menghitungnya di tiap putaran hanya membuang waktu.
3. Fungsi selesai. Ia berjalan di program utama, sedangkan fungsi tugas berjalan di latar belakang dan tidak boleh menyentuh proyek.

Kesalahan umum #

• Salah ketik nama isian atau kolom. Nama isian yang tak dikenal diabaikan diam-diam (Skrip 21.1), sedangkan nama kolom yang salah memicu KeyError. Perbaikan: salin nama dari parameterDefinitions() dan dari daftar kolom layer.
• Mengira layer yang gagal dibuka akan menimbulkan galat. Alat tetap jalan dengan layer kosong dan menghasilkan angka nol. Perbaikan: periksa isValid() dan baca baris pertama hasil.
• Urutan band terbalik. NDVI menjadi negatif tanpa pesan apa pun. Perbaikan: cek urutan band tiap dataset dan curigai tanda minus.
• Memakai satuan jarak yang salah. DISTANCE mengikuti CRS, jadi 30 pada CRS derajat berarti 30 derajat. Perbaikan: pastikan CRS metrik.
• Skrip mandiri yang tidak menunggu tugas latar belakang. Prosesnya berhenti mendadak sebelum tugas selesai. Perbaikan: tambahkan perulangan tunggu.

Ringkasan dan latihan #

• processing.run("ID_ALAT", {kamus parameter}) memanggil alat Processing. Hasilnya kamus, dan jawabannya diambil dengan kunci seperti ["OUTPUT"].
• Vektor sementara berupa layer. Raster sementara dan keluaran berupa berkas berupa teks jalur.
• Merangkai alat berarti memakai hasil satu alat sebagai masukan alat berikutnya.
• Batch berarti satu resep, banyak masukan: perulangan atas berkas (glob) atau atas fitur (getFeatures).
• QgsTask menjalankan pekerjaan berat di latar belakang. Fungsi tugas menghitung, fungsi selesai menyentuh proyek.

Latihan. Gabungkan Skrip 21.2 dan 21.5. Untuk tiap KPH, buat sempadan sungai 30 meter, lalu potong CHM dengan sempadan itu, dan cetak rata-rata tinggi tajuk di sempadan. Alat apa yang Anda butuhkan di setiap langkah, dan bagian mana yang dipinjam dari tiap skrip?

Tabel perbandingan: memanggil alat dengan klik dan dengan skrip #