Lewati ke isi
Profil penulisSeri Buku GIS Kehutanan dan Pertanian/ I4
Tampilkan bagian untuk:

BAB 2: Alat Processing dari Skrip: Satu Perintah, Satu Analisis

#Studi kasus: "Sempadan sungai tiap petak"

Kepala Seksi ingin tahu berapa hektare jalur sempadan sungai di tiap petak, dan berapa titik panas di dalamnya. Untuk latihan, lebar sempadan ditetapkan 30 meter. Di seri I1, Anda menjalankan Buffer lalu Intersection satu per satu dengan klik. Kini seluruh rangkaian itu Anda tulis sekali, lalu dijalankan dengan satu tombol.

#Konsep: alat Processing dalam tiga kalimat

Setiap alat di Processing Toolbox bisa dipanggil dari Python dengan processing.run. Anda memberi ID alat dan kamus parameter, lalu QGIS mengembalikan kamus hasil. Ini seperti mengisi formulir pesanan: menu memberi kolom kosong, sedangkan skrip mengisinya dengan tulisan.

Istilah baru bab ini:

  • ID alat: nama unik alat, misalnya native:buffer.
  • Parameter: isian alat, seperti masukan, jarak, dan keluaran.
  • Keluaran sementara (TEMPORARY_OUTPUT): hasil yang disimpan di memori atau folder sementara, bukan di berkas pilihan Anda.
  • Rantai alat: keluaran satu alat menjadi masukan alat berikutnya.
  • Statistik zonal: ringkasan nilai raster di dalam tiap poligon.
Ilustrasi 2.1: Satu panggilan, kamus masuk dan kamus keluar
Skema: kamus parameter masuk ke processing.run, kamus hasil keluar, dan hasilnya diambil dengan kunci OUTPUT

#Bagian A: QGIS

#Bagian A: processing.run di QGIS

#A.1 Menemukan ID alat dan parameternya

Pertanyaan pertama: apa nama alat dan parameternya? Skrip 2.1 menjawabnya dengan Python. Bagian pertamanya mencari ID alat dengan kata kunci. Bagian keduanya mencetak nama parameter alat Buffer. Skrip lalu menjalankan Buffer pada layer sungai.

Cara lain yang tersedia: processing.algorithmHelp("native:buffer") mencetak bantuan lengkap alat itu, termasuk nilai yang diterima tiap parameter. Panel Processing ► History juga mencatat tiap alat yang pernah Anda jalankan lewat dialog. [CEK: bentuk catatannya]

PYTHON (PYQGIS) 2.1
# [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()])

Hasil uji di QGIS 4.0.2:

HASIL UJI SKRIP 2.1
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']

#A.2 Keluaran sementara: layer atau teks?

Isian OUTPUT menentukan ke mana hasil dikirim. Isi dengan jalur berkas untuk menyimpan permanen. Isi dengan "TEMPORARY_OUTPUT" bila hasil hanya perantara. Hasil uji di QGIS 4.0.2 menunjukkan perbedaan ini:

Jenis keluaranIsi hasil["OUTPUT"]Cara memakainya
Vektor, TEMPORARY_OUTPUTLayer (QgsVectorLayer) bernama "Buffered"Langsung dipakai sebagai masukan alat lain
Raster, TEMPORARY_OUTPUTTeks jalur berkas sementaraBungkus dengan QgsRasterLayer(jalur, nama, "gdal")
Vektor, "memory:nama"LayerSama seperti vektor di atas
Jalur berkasTeks jalur berkas ituMuat sebagai layer bila perlu

Layer sementara tidak muncul di panel Layers sebelum Anda memanggil addMapLayer. Layer sementara juga hilang saat QGIS ditutup. Simpan hasil penting ke GeoPackage (Bab 7 mencontohkannya).

#A.3 Rantai alat vektor

Sempadan per petak butuh tiga alat. Buffer membuat jalur 30 m di tepi sungai. Intersection memotong jalur itu per petak. Count points in polygon menghitung titik panas. Keluaran alat pertama langsung menjadi masukan alat kedua, tanpa berkas perantara.

Ilustrasi 2.2: Rantai alat
Skema: sungai dibuffer 30 meter, diiris dengan petak KPH, lalu luas per KPH dihitung
PYTHON (PYQGIS) 2.2
# [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"]))

Hasil uji di QGIS 4.0.2:

HASIL UJI SKRIP 2.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

Angka dapat dicek tanpa komputer: KPH Alpha dan Beta masing-masing berluas 4 ha, dan KPH Gamma 8 ha. Jumlah titik panas 2, 5, dan 5 sama dengan 12 titik di Hotspot_Nasional. Total sempadan di tiga petak (1,18 + 1,39 + 1,50 = 4,07 ha) lebih kecil dari 4,48 ha pada Skrip 2.1, karena sebagian jalur buffer jatuh di luar batas area.

#A.4 Alat raster: NDVI, tinggi tajuk, dan statistik zonal

NDVI adalah angka kehijauan: (NIR − Merah) ÷ (NIR + Merah). CHM adalah tinggi tajuk, yaitu DSM dikurangi DTM. Keduanya dihitung dengan native:rastercalc. Dalam ekspresinya, band ditulis "NamaLayer@nomorBand". Karena itu nama layer pada QgsRasterLayer(...) penting.

Setelah itu native:zonalstatisticsfb menghitung rata-rata nilai raster di dalam tiap poligon KPH. Pilihan statistik berupa angka: 0 hitungan, 1 jumlah, 2 rata-rata, dan seterusnya (daftarnya ada di algorithmHelp).

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

Hasil uji di QGIS 4.0.2:

HASIL UJI SKRIP 2.3
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

Angka NDVI cocok dengan hitungan tangan dari data sintetis: Alpha memakai Merah 40 dan NIR 160, jadi (160 − 40) ÷ (160 + 40) = 0,60. Beta memakai Merah 60 dan NIR 150, jadi 0,43. Gamma terdiri atas dua separuh (sekitar 0,52 dan 0,13) sehingga rata-ratanya 0,33.

#Bagian B: ArcGIS Pro

#Bagian B: Alat geoprocessing lewat ArcPy di ArcGIS Pro

Di ArcPy, setiap alat adalah fungsi: arcpy.analysis.Buffer, arcpy.analysis.Intersect, dan arcpy.analysis.SpatialJoin. Parameter ditulis berurutan atau dengan nama. Hasil antara disimpan di ruang kerja memory, misalnya memory\sempadan. Ruang ini padanan TEMPORARY_OUTPUT: cepat, tetapi memakai RAM dan hilang saat Pro ditutup.

Beda penting dengan QGIS. Menghitung titik per petak memakai Spatial Join (kolom Join_Count otomatis terbentuk). Statistik zonal memakai ZonalStatisticsAsTable dari Spatial Analyst. NDVI memakai Float dan band ditulis sebagai jalur\Band_n. [CEK]

PYTHON (ARCPY, PYTHON 3)
# [SKRIP 2.1-Pro: Rantai alat: sempadan sungai, titik panas, dan NDVI per KPH dengan ArcPy]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Zonal Statistics as Table dan Float butuh ekstensi Spatial Analyst.
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")

# Keluaran sementara: ruang kerja "memory" (padanan TEMPORARY_OUTPUT di QGIS)
arcpy.analysis.Buffer("Sungai.shp", r"memory\sempadan", "30 Meters", dissolve_option="ALL")
arcpy.analysis.Intersect(["Batas_KPH.shp", r"memory\sempadan"], r"memory\irisan")
arcpy.analysis.SpatialJoin("Batas_KPH.shp", "Hotspot_Nasional.shp", r"memory\kph_titik",
                           "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="INTERSECT")   # menambah Join_Count

sempadan = {}
with arcpy.da.SearchCursor(r"memory\irisan", ["NAMA_KPH", "SHAPE@AREA"]) as k:
    for nama, luas in k:
        sempadan[nama] = sempadan.get(nama, 0) + luas / 10000
with arcpy.da.SearchCursor(r"memory\kph_titik", ["NAMA_KPH", "Join_Count", "SHAPE@AREA"]) as k:
    for nama, jml, luas in k:
        print("%s: luas %.2f ha, sempadan %.2f ha, titik panas %d"
              % (nama, luas / 10000, sempadan.get(nama, 0), jml))

# NDVI = (NIR - Merah) / (NIR + Merah); band 4 = NIR, band 1 = Merah. Nama band ditulis jalur\Band_n.
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)
ndvi.save(HASIL + r"\NDVI_KPH.tif")
ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", ndvi, r"memory\ndvi_zonal", "DATA", "MEAN")
with arcpy.da.SearchCursor(r"memory\ndvi_zonal", ["NAMA_KPH", "MEAN"]) as k:
    for nama, rata in k:
        print("%s: NDVI rata-rata %.3f" % (nama, rata))
arcpy.management.Delete(r"memory\sempadan")      # ruang memory memakai RAM: bersihkan setelah selesai
arcpy.CheckInExtension("Spatial")

#Bagian C: ArcMap 10.8

#Bagian C: Alat geoprocessing lewat ArcPy di ArcMap 10.8

Di ArcMap, nama alat memakai bentuk lama: arcpy.Buffer_analysis, arcpy.Intersect_analysis, dan arcpy.SpatialJoin_analysis. Ruang kerja sementaranya bernama in_memory, bukan memory. Pembagian bilangan memakai 10000.0 karena Python 2.7.

PYTHON (ARCPY, PYTHON 2.7)
# -*- coding: utf-8 -*-
# [SKRIP 2.1-ArcMap: Rantai alat: sempadan sungai, titik panas, dan NDVI per KPH (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Zonal Statistics as Table dan Float butuh ekstensi Spatial Analyst.
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")

# Keluaran sementara di ArcMap: ruang kerja in_memory (di Pro namanya memory)
arcpy.Buffer_analysis("Sungai.shp", r"in_memory\sempadan", "30 Meters", dissolve_option="ALL")
arcpy.Intersect_analysis(["Batas_KPH.shp", r"in_memory\sempadan"], r"in_memory\irisan")
arcpy.SpatialJoin_analysis("Batas_KPH.shp", "Hotspot_Nasional.shp", r"in_memory\kph_titik",
                           "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="INTERSECT")   # menambah Join_Count

sempadan = {}
with arcpy.da.SearchCursor(r"in_memory\irisan", ["NAMA_KPH", "SHAPE@AREA"]) as k:
    for nama, luas in k:
        sempadan[nama] = sempadan.get(nama, 0) + luas / 10000.0
with arcpy.da.SearchCursor(r"in_memory\kph_titik", ["NAMA_KPH", "Join_Count", "SHAPE@AREA"]) as k:
    for nama, jml, luas in k:
        print("%s: luas %.2f ha, sempadan %.2f ha, titik panas %d"
              % (nama, luas / 10000.0, sempadan.get(nama, 0), jml))

# NDVI = (NIR - Merah) / (NIR + Merah); band 4 = NIR, band 1 = Merah
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)
ndvi.save(HASIL + r"\NDVI_KPH.tif")
ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", ndvi, r"in_memory\ndvi_zonal", "DATA", "MEAN")
with arcpy.da.SearchCursor(r"in_memory\ndvi_zonal", ["NAMA_KPH", "MEAN"]) as k:
    for nama, rata in k:
        print("%s: NDVI rata-rata %.3f" % (nama, rata))
arcpy.Delete_management("in_memory")      # bersihkan seluruh ruang in_memory
arcpy.CheckInExtension("Spatial")

#Cek paham

  1. Alat vektor dijalankan dengan OUTPUT berisi "TEMPORARY_OUTPUT". Apa isi hasil["OUTPUT"]?
  2. Layer raster bernama Citra_KPH. Bagaimana menulis band keempat dalam ekspresi raster calculator?
  3. Mengapa DISTANCE: 30 pada layer berkoordinat derajat memberi hasil yang salah?

Jawaban:

  1. Sebuah layer vektor, yang bisa langsung dipakai alat berikutnya.
  2. "Citra_KPH@4", yaitu nama layer, tanda @, lalu nomor band.
  3. Satuan jarak mengikuti CRS layer. Pada derajat, angka 30 berarti 30 derajat. Ubah ke CRS metrik (UTM) lebih dulu.

#Kesalahan umum

  • Lupa `["OUTPUT"]`. processing.run mengembalikan kamus, bukan layer. Ambil isinya dengan kunci "OUTPUT".
  • Menganggap keluaran raster sementara sebagai layer. Isinya teks jalur. Bungkus dengan QgsRasterLayer.
  • Nama layer di ekspresi raster tidak cocok. Pada uji, nama yang salah menghasilkan galat perhitungan. Samakan nama di ekspresi dengan nama layer.
  • Hasil tidak muncul di peta. Layer sementara baru tampil setelah addMapLayer.

#Ringkasan dan latihan

Ringkasan: processing.run(ID, kamus) menjalankan alat apa pun dan mengembalikan kamus hasil. TEMPORARY_OUTPUT cocok untuk hasil antara. Vektor kembali sebagai layer, raster sebagai jalur teks. Alat bisa dirantai tanpa berkas perantara.

Latihan:

  1. Ubah lebar sempadan menjadi 50 m pada Skrip 2.2. Perhatikan perubahan luas sempadan tiap petak. (Total tiga petak seharusnya 6,55 ha.)
  2. Tambahkan statistik 5 (minimum) dan 6 (maksimum) pada STATISTICS di Skrip 2.3. Cetak nilai NDVI terkecil dan terbesar per petak.

#Tabel perbandingan: alat dari skrip

HalQGISArcGIS ProArcMap 10.8
Memanggil alatprocessing.run("native:buffer", {...})arcpy.analysis.Buffer(...)arcpy.Buffer_analysis(...)
Hasil antaraTEMPORARY_OUTPUTmemory\namain_memory\nama
Hitung titik per poligonnative:countpointsinpolygonSpatialJoin + Join_Count [CEK]SpatialJoin_Analysis + Join_Count [CEK]
Raster calculatornative:rastercalcAljabar peta di arcpy.sa [CEK]Aljabar peta di arcpy.sa [CEK]
Statistik zonalnative:zonalstatisticsfbZonalStatisticsAsTable [CEK]ZonalStatisticsAsTable [CEK]
EkstensiTidak perluSpatial Analyst [CEK]Spatial Analyst [CEK]