BAB 25: Dari Lemari Berkas ke Gudang Data: PostGIS, Jaringan, Hidrologi, dan Kesesuaian Lahan dengan Skrip

Studi kasus: "Datanya ada di mana-mana, hasilnya beda-beda" #

Setiap analis di KPH Contoh menyimpan salinan data sendiri. Petak ada di satu berkas, jalan di berkas lain, dan angka luas dihitung dengan cara yang berbeda-beda. Kepala Seksi minta satu "gudang data" bersama. Ia juga minta jawaban cepat: rute patroli tercepat antar pos, luas daerah aliran sungai, peta curah hujan, dan peta kesesuaian lahan. Analis Muda menjawabnya dengan skrip, supaya setiap angka bisa dicek dan diulang.

Konsep: gudang data dalam tiga kalimat #

Berkas GeoPackage itu seperti lemari arsip pribadi: enak dibawa, tetapi tiap orang punya salinannya sendiri. PostGIS adalah gudang arsip pusat. Ia berupa program server yang melayani banyak orang sekaligus, dan Anda memintanya lewat surat bernama SQL (bahasa untuk bertanya kepada basis data).

Lima istilah baru yang perlu Anda kenal sekarang:

• Basis data spasial: tempat menyimpan tabel yang kolomnya boleh berisi bentuk (titik, garis, poligon). PostGIS menambahkan kemampuan itu ke basis data PostgreSQL.
• Skema: laci besar di dalam basis data. Semua tabel bab ini ada di laci bernama kph.
• SQL spasial: SQL yang mengerti bentuk, misalnya ST_Area (luas) dan ST_Intersects (bersinggungan).
• Variabel lingkungan: pengaturan yang disimpan di komputer, di luar skrip, yang boleh dibaca skrip. Cara membacanya os.environ (Bab 8).
• Tanpa jendela (headless): QGIS dijalankan dari terminal, tanpa membuka tampilan grafisnya.

Sebelum mulai: tiga persiapan #

Pertama, paket data. Siapkan folder paket-m2. Isinya delapan berkas GeoPackage (Petak.gpkg, Tanah.gpkg, Jalan.gpkg, Jalan_celah.gpkg, Fasilitas.gpkg, Kejadian.gpkg, Stasiun_Hujan.gpkg, Batas_Wilayah.gpkg), satu peta ketinggian DEM_10m.tif (200 x 200 sel, sel 10 meter, wilayah 2 x 2 km), dan folder hasil untuk keluaran. Semua data sintetis, bukan kondisi lapangan nyata.

Kedua, PostGIS (untuk skrip tertentu). Skrip 25.1, 25.2, 25.6, dan 25.10 butuh server PostgreSQL dengan PostGIS (dan pgRouting untuk 25.6), serta basis data bernama kph_contoh dengan skema kph. Skrip 25.3 tetap jalan tanpa PostGIS, hanya bagian pembandingnya dilewati. Skrip 25.4 sampai 25.9 tidak butuh PostGIS. Skrip 25.7 dan 25.8 butuh GRASS. Pada pemasangan QGIS 4.0.2 di komputer uji, GRASS sudah ikut terpasang. Uji untuk buku ini memakai QGIS 4.0.2, PostgreSQL 17, PostGIS 3.6, dan pgRouting 4.0. Bila server Anda berbeda versi, cocokkan hasilnya. [CEK]

Ketiga, koneksi lewat variabel lingkungan. Skrip tidak memuat alamat server atau kata sandi. Skrip membaca empat pengaturan dari komputer Anda: PGHOST (alamat server), PGPORT (nomor port), PGUSER (nama pengguna), dan PGDATABASE (nama basis data). Isi dulu di terminal, lalu jalankan skripnya. Contoh di Windows (terminal Command Prompt):

set PGHOST=localhost
set PGPORT=5432
set PGUSER=nama_pengguna_anda
set PGDATABASE=kph_contoh
"C:\Program Files\QGIS 4.0.2\bin\python-qgis.bat" m2_01b_impor.py D:\Latihan\paket-m2

Baris terakhir menjalankan skrip dengan Python milik QGIS (lihat Bab 3). Kata D:\Latihan\paket-m2 adalah argumen: folder paket data, yang dibaca skrip lewat sys.argv[1]. Di macOS dan Ubuntu, perintahnya berbeda. [CEK]

Berkas pembantu: _inisialisasi.py #

Sebagian besar skrip bab ini dimulai dengan import _inisialisasi. Berkas pendek ini menyalakan QGIS tanpa jendela. Bayangkan menyalakan dapur sebelum memasak: bila kompor tidak dinyalakan, resep sebagus apa pun tidak jalan. Di Konsol Python QGIS, dapurnya sudah menyala, jadi berkas ini tidak diperlukan.

Skrip utuh:

# -*- coding: utf-8 -*-
"""Berkas bantu M2: menyalakan QGIS tanpa jendela (headless) agar skrip bisa dijalankan dengan python-qgis.bat.
Penulis: Badar Mubarok Yogaswara. Di Python Console QGIS berkas ini TIDAK diperlukan."""
import os
import sys
from qgis.core import QgsApplication, Qgis

os.environ.setdefault("QT_QPA_PLATFORM", "offscreen")
qgs = QgsApplication([], False)
qgs.initQgis()
sys.path.append(os.path.join(QgsApplication.prefixPath(), "python", "plugins"))
import processing  # noqa: E402
from processing.core.Processing import Processing  # noqa: E402

Processing.initialize()

# GRASS: bila dijalankan tanpa jendela QGIS, beri tahu lokasi GRASS lewat variabel GISBASE
# (di jendela QGIS sendiri hal ini sudah otomatis dan bagian ini boleh dilewati).
_grass = os.path.join(os.path.dirname(QgsApplication.prefixPath()), "grass", "grass84")
if os.path.isdir(_grass) and "GISBASE" not in os.environ:
    os.environ["GISBASE"] = _grass

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Berkas ini tidak mencetak apa-apa. Buktinya bisa Anda lihat dengan dua baris pemeriksaan berikut, yang Anda simpan dalam berkas di folder yang sama:

import _inisialisasi
import os
from qgis.core import QgsApplication
print('alat GRASS ada:', QgsApplication.processingRegistry().algorithmById('grass:r.watershed') is not None)
print('GISBASE terisi:', bool(os.environ.get('GISBASE')))

alat GRASS ada: True
GISBASE terisi: True

Bagian yang boleh diubah. Baris 19, tulisan "grass84". Itu nama folder GRASS di instalasi QGIS 4.0.2 pada komputer penulis. Bila instalasi Anda memakai versi GRASS lain, ganti angkanya.

Jebakan.

• Versi GRASS berbeda. Bila folder grass84 tidak ada di komputer Anda, GISBASE tidak terisi. Alat grass: lalu gagal dengan QgsProcessingException: GRASS GIS folder is not configured. Please configure it before running GRASS GIS algorithms. (diuji dengan mengganti nama folder menjadi yang tidak ada).
• Tidak ada penutup. Berkas menyalakan QGIS tetapi tidak pernah mematikannya (tidak ada exitQgis()). Untuk skrip pendek, ini tidak masalah.
• Satu berkas, satu folder. Skrip lain menambahkan folder skripnya sendiri ke sys.path supaya import _inisialisasi ketemu. Pindahkan keduanya bersama-sama.

Skrip 25.1: Memasukkan berkas ke PostGIS dengan GDAL #

Berkas m2_01b_impor.py. Di dalam berkasnya tertulis "M2 Bab 1". Skrip ini memindahkan tujuh GeoPackage ke gudang. Pekerjaannya sama dengan perintah ogr2ogr, tetapi ditulis di Python.

Tujuan. Memasukkan tujuh berkas GeoPackage menjadi tujuh tabel di skema kph pada PostGIS.

Bahan. Paket paket-m2 (tujuh GeoPackage), server PostGIS dengan skema kph yang sudah dibuat (oleh m2_01a_skema.sql, lihat akhir bab), dan variabel PGHOST, PGPORT, PGUSER, PGDATABASE. Perintah: python-qgis.bat m2_01b_impor.py <folder paket-m2>.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 1: memasukkan berkas GeoPackage ke PostGIS memakai GDAL (setara ogr2ogr).
Penulis: Badar Mubarok Yogaswara
Pemakaian: python-qgis.bat m2_01b_impor.py <folder paket-m2>
Koneksi dibaca dari variabel lingkungan standar PostgreSQL: PGHOST, PGPORT, PGUSER, PGPASSWORD, PGDATABASE.
Kata sandi tidak pernah ditulis di skrip."""
import os
import sys
from osgeo import gdal

gdal.UseExceptions()
paket = sys.argv[1]
tujuan = "PG:dbname=%s" % os.environ.get("PGDATABASE", "kph_contoh")
# (berkas, layer sumber, nama tabel, jenis geometri)
daftar = [("Petak.gpkg", "Petak", "petak", "MULTIPOLYGON"),
          ("Tanah.gpkg", "Tanah", "tanah", "MULTIPOLYGON"),
          ("Jalan.gpkg", "Jalan", "jalan", "MULTILINESTRING"),
          ("Fasilitas.gpkg", "Fasilitas", "fasilitas", "POINT"),
          ("Kejadian.gpkg", "Kejadian", "kejadian", "POINT"),
          ("Stasiun_Hujan.gpkg", "Stasiun_Hujan", "stasiun_hujan", "POINT"),
          ("Batas_Wilayah.gpkg", "Batas_Wilayah", "batas_wilayah", "MULTIPOLYGON")]
for berkas, layer, tabel, gtype in daftar:
    opsi = gdal.VectorTranslateOptions(
        format="PostgreSQL", layerName=tabel, geometryType=gtype, accessMode="overwrite",
        layerCreationOptions=["SCHEMA=kph", "GEOMETRY_NAME=geom", "FID=gid", "SPATIAL_INDEX=GIST", "OVERWRITE=YES"],
        srcSRS="EPSG:32749", dstSRS="EPSG:32749", reproject=False)
    gdal.VectorTranslate(tujuan, os.path.join(paket, berkas), options=opsi, layers=[layer])
    print("dimasukkan:", tabel)

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Layar menampilkan tujuh baris:

dimasukkan: petak
dimasukkan: tanah
dimasukkan: jalan
dimasukkan: fasilitas
dimasukkan: kejadian
dimasukkan: stasiun_hujan
dimasukkan: batas_wilayah

Periksa di basis data. Jumlah baris tabel (dihitung dengan psql) adalah petak 25, tanah 8, jalan 37, fasilitas 6, kejadian 150, stasiun_hujan 30, dan batas_wilayah 1.

SELECT count(*) FROM kph.petak;

Hasilnya 25.

Bagian yang boleh diubah.

• Baris 13, nama basis data pengganti. Ganti "kph_contoh" dengan nama basis data Anda, atau isi PGDATABASE.
• Baris 15-21, daftar. Tambah atau hapus tuple untuk memasukkan berkas lain. Isi keempat bagian tuple.
• Baris 25, nama skema. Ganti SCHEMA=kph bila skema Anda bernama lain.
• Baris 25, pengaman ulang. OVERWRITE=YES membuat skrip aman dijalankan berulang. Hapus pilihan ini bila Anda ingin skrip berhenti ketika tabel sudah ada.
• Baris 26, sistem koordinat. Ganti kode EPSG:32749 bila data Anda berada di zona lain.

Jebakan.

• Lupa argumen folder. Tanpa <folder paket-m2>, baris 12 berhenti dengan IndexError: list index out of range.
• Skema belum ada. Bila skema kph belum dibuat, GDAL berhenti dengan RuntimeError: Terminating translation prematurely ... May be caused by: Schema "kph" does not exist. Jalankan m2_01a_skema.sql dulu.
• Dijalankan dua kali. Berkat OVERWRITE=YES di baris 25, jalan kedua berhasil: tabel lama diganti dengan isi berkas asal. Akibatnya, nilai yang sudah Anda ubah langsung di tabel akan hilang (diuji: nilai yang sengaja diubah kembali ke semula). Tanpa pilihan itu, jalan kedua gagal dengan RuntimeError ... Layer kph.petak already exists, CreateLayer failed. Use the layer creation option OVERWRITE=YES to replace it. Pilihan accessMode="overwrite" saja tidak cukup. Percobaannya ada di latihan variasi nomor 2.
• Nama basis data salah. Pada komputer penulis, server berbahasa Indonesia mengirim pesan dengan huruf yang bukan UTF-8, sehingga yang tampil UnicodeDecodeError: 'utf-8' codec can't decode byte 0xab .... Penyebab sebenarnya: basis data PGDATABASE tidak ada. Cek namanya dengan psql.
• Sandi. Skrip tidak mengisi kata sandi. Bila server memintanya, siapkan berkas pgpass.

Latihan variasi.

1. Masukkan hanya Petak.gpkg, dengan nama tabel petak_uji.
2. Buktikan sendiri gunanya OVERWRITE=YES: hapus pilihan itu, lalu jalankan skrip dua kali.

Contoh jawaban 1: ganti baris 15-21 dengan satu tuple. Tabel kph.petak_uji terbentuk dan SELECT count(*) FROM kph.petak_uji; menghasilkan 25.

daftar = [("Petak.gpkg", "Petak", "petak_uji", "MULTIPOLYGON")]

Contoh jawaban 2: pakai daftar semula, lalu ubah baris 25 menjadi seperti berikut (tanpa OVERWRITE=YES). Jalan pertama pada basis data kosong berhasil. Jalan kedua gagal dengan Layer kph.petak already exists, dan tabel dari jalan pertama tetap utuh. Kembalikan "OVERWRITE=YES" untuk memulihkan skrip.

        layerCreationOptions=["SCHEMA=kph", "GEOMETRY_NAME=geom", "FID=gid", "SPATIAL_INDEX=GIST"],

Skrip 25.2: Memuat tabel PostGIS di QGIS dengan skrip #

Berkas m2_01c_qgis_muat.py. Skrip ini membalik arah: kini QGIS yang membaca dari gudang. Ia juga menunjukkan cara memuat sebagian data dengan saringan SQL, dan mengekspor satu layer lewat alat Processing.

Tujuan. Memuat tabel kph.petak sebagai layer QGIS, membuat layer berisi petak jati saja, lalu mengekspor Fasilitas.gpkg ke tabel baru.

Bahan. Hasil Skrip 25.1 (tabel kph.petak sudah ada), Fasilitas.gpkg, empat variabel PG*. Jalankan dari terminal dengan folder paket sebagai argumen.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 1: memuat tabel PostGIS sebagai layer lewat PyQGIS, dan mengekspor layer ke PostGIS lewat Processing.
Penulis: Badar Mubarok Yogaswara. Koneksi dari variabel lingkungan PGHOST, PGPORT, PGUSER, PGDATABASE (tanpa kata sandi di skrip)."""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401  (hanya perlu bila dijalankan di luar jendela QGIS)
from qgis.core import QgsDataSourceUri, QgsVectorLayer, QgsProject
import processing

paket = sys.argv[1]
uri = QgsDataSourceUri()
uri.setConnection(os.environ["PGHOST"], os.environ["PGPORT"], os.environ["PGDATABASE"], os.environ["PGUSER"], "")
uri.setDataSource("kph", "petak", "geom", "", "gid")          # skema, tabel, kolom geometri, filter, kolom kunci
lyr = QgsVectorLayer(uri.uri(False), "Petak (PostGIS)", "postgres")
print("sah:", lyr.isValid(), "| jumlah fitur:", lyr.featureCount(), "| CRS:", lyr.crs().authid(), "| geometri:", lyr.wkbType())
if not lyr.isValid():
    sys.exit("Layer kph.petak tidak sah: periksa PGHOST, PGPORT, PGUSER, PGDATABASE dan skema kph.")

# layer dengan filter SQL (hanya petak jati)
uri.setSql("jenis_tegakan = 'Jati'")
lyr2 = QgsVectorLayer(uri.uri(False), "Petak jati", "postgres")
print("petak jati:", lyr2.featureCount())
QgsProject.instance().addMapLayer(lyr)

# ekspor satu layer GeoPackage ke tabel baru lewat Processing (alat GDAL)
hasil = processing.run("gdal:importvectorintopostgisdatabasenewconnection", {
    "INPUT": os.path.join(paket, "Fasilitas.gpkg|layername=Fasilitas"), "GTYPE": 3,
    "HOST": os.environ["PGHOST"], "PORT": os.environ["PGPORT"], "USER": os.environ["PGUSER"], "DBNAME": os.environ["PGDATABASE"],
    "PASSWORD": "", "SCHEMA": "kph", "TABLE": "fasilitas_uji", "PK": "gid", "GEOCOLUMN": "geom", "DIM": 0,
    "OVERWRITE": True, "APPEND": False, "ADDFIELDS": False, "LAUNDER": False, "INDEX": False, "SKIPFAILURES": False,
    "MAKEVALID": False, "PROMOTETOMULTI": False, "PRECISION": True})
print("ekspor selesai:", hasil)

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Layar menampilkan:

sah: True | jumlah fitur: 25 | CRS: EPSG:32749 | geometri: 6
petak jati: 10
ekspor selesai: {}

Angka 6 pada baris pertama adalah kode jenis bentuk MultiPolygon (poligon majemuk). Kamus kosong {} di baris terakhir wajar, karena alat ekspor tidak punya keluaran. Buktinya, tabel kph.fasilitas_uji terisi 6 baris (dicek dengan psql). Pada uji ini, QGIS juga mencetak peringatan DeprecationWarning (fungsi lama). Peringatan itu boleh diabaikan.

Bagian yang boleh diubah.

• Baris 14, skema dan tabel. Ganti "kph", "petak" untuk memuat tabel lain.
• Baris 18, pesan pengaman. Ganti teks di dalam sys.exit(...) dengan kalimat yang paling membantu Anda.
• Baris 21, saringan. Ganti jenis_tegakan = 'Jati' dengan syarat SQL lain.
• Baris 28, berkas ekspor. Ganti berkas dan nama layer.
• Baris 30, nama tabel. Ganti fasilitas_uji agar tidak menimpa tabel lama.

Jebakan.

• Variabel belum diisi. Bila PGHOST kosong, baris 13 berhenti dengan KeyError: 'PGHOST'.
• Basis data tidak ada. Layer tidak sah, jadi layar menampilkan sah: False | jumlah fitur: 0. Lalu baris 17-18 menghentikan skrip dengan pesan Layer kph.petak tidak sah: periksa PGHOST, PGPORT, PGUSER, PGDATABASE dan skema kph. Tahap ekspor tidak sempat jalan. Baca baris sah: dan pesan itu, lalu periksa keempat variabel PG* serta skema kph. (Diuji dengan nama basis data yang tidak ada.)
• Layer kedua tidak dimuat. lyr2 hanya dihitung. Tambahkan addMapLayer(lyr2) bila ingin melihatnya di peta.
• Menimpa tabel. OVERWRITE: True mengganti tabel fasilitas_uji bila sudah ada.
• Hanya terminal. sys.argv[1] tidak ada di Konsol QGIS. Di Konsol, ganti baris 6-7 dan 11 dengan jalur manual. [CEK]

Latihan variasi.

1. Muat hanya petak mahoni.
2. Tambahkan layer petak jati ke peta.

Contoh jawaban 1: ganti baris 21 dengan satu baris berikut. Layar menampilkan petak jati: 5 (tulisan petak jati ikut tercetak apa adanya, sebab baris 23 tidak diganti).

uri.setSql("jenis_tegakan = 'Mahoni'")

Contoh jawaban 2: tambahkan satu baris setelah baris 24.

QgsProject.instance().addMapLayer(lyr2)

Skrip 25.3: Kueri yang sama di GeoPackage dan PostGIS #

Berkas m2_02c_silang_geopackage.py. "Silang" berarti mencocokkan dua cara. Skrip ini menjalankan tiga pertanyaan spasial di GeoPackage, lalu (bila PostGIS ada) di PostGIS, dan membandingkan jawabannya. Bila angkanya sama, Anda boleh percaya keduanya.

Tujuan. Menjawab tiga pertanyaan dengan SQL di GeoPackage, lalu membandingkan dengan jawaban PostGIS.

Bahan. paket-m2 (tujuh GeoPackage) dan folder hasil yang sudah ada. PostGIS dengan hasil Skrip 25.1 bersifat pilihan. Perintah: python-qgis.bat m2_02c_silang_geopackage.py <folder paket-m2>.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 2: menjalankan SQL spasial yang sama pada GeoPackage (dialek SQLite/SpatiaLite lewat GDAL),
lalu membandingkan angkanya dengan hasil PostGIS. Penulis: Badar Mubarok Yogaswara.
Pemakaian: python-qgis.bat m2_02c_silang_geopackage.py <folder paket-m2>
Bila PostGIS tidak tersedia, bagian pembanding dilewati (hanya mencetak angka GeoPackage)."""
import os
import sys
from osgeo import gdal, ogr

gdal.UseExceptions()
ogr.UseExceptions()
paket = sys.argv[1]
gpkg = os.path.join(paket, "hasil", "kph_contoh.gpkg")

# 1) kumpulkan semua layer ke satu GeoPackage (sekali saja)
if not os.path.exists(gpkg):
    for berkas, layer in [("Petak", "Petak"), ("Tanah", "Tanah"), ("Jalan", "Jalan"), ("Fasilitas", "Fasilitas"),
                          ("Kejadian", "Kejadian"), ("Stasiun_Hujan", "Stasiun_Hujan"), ("Batas_Wilayah", "Batas_Wilayah")]:
        gdal.VectorTranslate(gpkg, os.path.join(paket, berkas + ".gpkg"),
                             options=gdal.VectorTranslateOptions(format="GPKG", layerName=layer, accessMode="append"))
ds = gdal.OpenEx(gpkg)


def tanya(sql):
    lyr = ds.ExecuteSQL(sql, dialect="SQLite")
    baris = [tuple(f.GetField(i) for i in range(f.GetFieldCount())) for f in lyr]
    ds.ReleaseResultSet(lyr)
    return baris


sqlite = {
    "total_ha": tanya("SELECT ROUND(SUM(ST_Area(geom)) / 10000.0, 2) FROM Petak")[0][0],
    "per_petak_teratas": tanya("SELECT p.kode, COUNT(k.fid) AS n FROM Petak p LEFT JOIN Kejadian k ON ST_Intersects(p.geom, k.geom) "
                               "GROUP BY p.kode ORDER BY n DESC, p.kode LIMIT 3"),
    "tebangan_dekat_jalan": tanya("SELECT COUNT(*) FROM Kejadian k WHERE k.jenis = 'Tebangan liar' AND "
                                  "EXISTS (SELECT 1 FROM Jalan j WHERE ST_Distance(k.geom, j.geom) <= 100)")[0][0],
}
print("GeoPackage/SpatiaLite:", sqlite)

try:
    import psycopg2
    con = psycopg2.connect(dbname=os.environ["PGDATABASE"])      # host, port, user dari PGHOST, PGPORT, PGUSER
    cur = con.cursor()
    cur.execute("SELECT ROUND((SUM(ST_Area(geom)) / 10000.0)::numeric, 2) FROM kph.petak")
    pg_total = float(cur.fetchone()[0])
    cur.execute("SELECT p.kode, COUNT(k.gid) n FROM kph.petak p LEFT JOIN kph.kejadian k ON ST_Intersects(p.geom, k.geom) "
                "GROUP BY p.kode ORDER BY n DESC, p.kode LIMIT 3")
    pg_top = [tuple(r) for r in cur.fetchall()]
    cur.execute("SELECT COUNT(*) FROM kph.kejadian k WHERE k.jenis = 'Tebangan liar' AND "
                "EXISTS (SELECT 1 FROM kph.jalan j WHERE ST_DWithin(k.geom, j.geom, 100))")
    pg_dekat = cur.fetchone()[0]
    print("PostGIS              :", {"total_ha": pg_total, "per_petak_teratas": pg_top, "tebangan_dekat_jalan": pg_dekat})
    cocok = (abs(pg_total - sqlite["total_ha"]) < 0.01 and pg_top == sqlite["per_petak_teratas"] and pg_dekat == sqlite["tebangan_dekat_jalan"])
    print("HASIL SAMA:", cocok)
except Exception as galat:                                          # PostGIS tidak ada atau tidak terhubung
    print("Pembanding PostGIS dilewati:", galat)

Penjelasan baris demi baris. Skrip ini punya tiga bagian: menyiapkan satu berkas gabungan, menjawab di GeoPackage, lalu menjawab di PostGIS. Nomor baris mengikuti berkas di atas.

Tiga pertanyaan dalam bahasa biasa: (1) berapa hektare luas seluruh petak? (2) tiga petak mana yang memuat kejadian terbanyak? (3) berapa kejadian "Tebangan liar" yang jaraknya 100 meter atau kurang dari jalan? Baris 33 memakai LEFT JOIN: semua petak ikut, walau tanpa kejadian. GROUP BY mengelompokkan per petak, ORDER BY ... LIMIT 3 mengambil tiga teratas.

Hasil. Dengan PostGIS menyala, layar menampilkan:

GeoPackage/SpatiaLite: {'total_ha': 400.0, 'per_petak_teratas': [('P-04', 13), ('P-09', 12), ('P-12', 12)], 'tebangan_dekat_jalan': 59}
PostGIS              : {'total_ha': 400.0, 'per_petak_teratas': [('P-04', 13), ('P-09', 12), ('P-12', 12)], 'tebangan_dekat_jalan': 59}
HASIL SAMA: True

Tanpa PostGIS (misalnya server mati), baris GeoPackage tetap tampil, lalu Pembanding PostGIS dilewati: diikuti pesan galat.

Bagian yang boleh diubah.

• Baris 32-36, isi kueri. Ganti dengan pertanyaan Anda, lalu ubah kueri PostGIS yang setara (baris 44-51).
• Baris 36 dan 50, ambang jarak. Angka 100 (meter) ada di dua tempat. Ganti keduanya bersama-sama.
• Baris 53, toleransi. Angka 0.01 menentukan seberapa kecil selisih luas yang dianggap sama.

Jebakan.

• Folder hasil tidak ada. Skrip tidak membuat folder hasil. Penggabungan di baris 19-20 berhenti dengan RuntimeError: GPKG driver failed to create ...hasil\kph_contoh.gpkg ... unable to open database file.
• Server tidak terjangkau. Pesan yang tampil: Pembanding PostGIS dilewati: connection to server at "localhost" ... failed: Connection refused. Perhatikan, skrip tetap berjalan.
• Penjaring galat terlalu lebar. Baris except Exception menangkap semuanya. Salah ketik di baris 41-54 pun tampil sebagai "dilewati". Baca pesan setelah titik dua.
• Pesan menyesatkan di komputer berbahasa Indonesia. Basis data yang tidak ada bisa menampilkan 'utf-8' codec can't decode byte 0xab (pesan server bukan UTF-8), bukan pesan "database tidak ada".
• Sambungan tidak ditutup. con.close() tidak ada. Untuk skrip pendek tidak masalah.

Latihan variasi. Ganti ambang jarak dari 100 menjadi 50 meter di kedua kueri. Contoh jawaban: ganti <= 100 menjadi <= 50 di baris 36, dan ST_DWithin(k.geom, j.geom, 100) menjadi 50 di baris 50. Layar menampilkan 'tebangan_dekat_jalan': 46 di baris GeoPackage dan PostGIS, dan HASIL SAMA: True.

Skrip 25.4: Hitungan yang sama tanpa basis data #

Berkas m2_02d_qgis_setara.py. Tanpa PostGIS pun Anda bisa menghitung jumlah kejadian per petak. Skrip ini memakai dua alat Processing QGIS, lalu membandingkan hasilnya. Skrip ini juga memperkenalkan lambda.

Tujuan. Menghitung kejadian per petak dengan dua cara di QGIS dan memeriksa bahwa hasilnya sama.

Bahan. Petak.gpkg dan Kejadian.gpkg dari paket-m2. Tidak butuh PostGIS. Tidak menulis berkas.

Konsep singkat: lambda. Lambda adalah fungsi mini yang ditulis dalam satu baris, tanpa nama resmi. Bayangkan secarik catatan berisi satu aturan: "ambil angka ini, lalu hitung begini". Bentuknya lambda masukan: hasil. Di skrip ini lambda dipakai sebagai aturan pengurutan.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 2: padanan Q2 (kejadian per petak) lewat alat Processing QGIS tanpa basis data:
(a) Count points in polygon, (b) Execute SQL (kueri pada layer maya, dialek SQLite; layer masukan bernama input1, input2, dan seterusnya).
Penulis: Badar Mubarok Yogaswara. Pemakaian: python-qgis.bat m2_02d_qgis_setara.py <folder paket-m2>"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
from qgis.core import QgsVectorLayer

paket = sys.argv[1]
petak = QgsVectorLayer(os.path.join(paket, "Petak.gpkg|layername=Petak"), "Petak", "ogr")
kej = QgsVectorLayer(os.path.join(paket, "Kejadian.gpkg|layername=Kejadian"), "Kejadian", "ogr")

a = processing.run("native:countpointsinpolygon", {"POLYGONS": petak, "POINTS": kej, "FIELD": "jumlah", "OUTPUT": "memory:"})["OUTPUT"]
hasil_a = sorted(((f["kode"], f["jumlah"]) for f in a.getFeatures()), key=lambda t: (-t[1], t[0]))[:3]
print("Count points in polygon:", hasil_a)

sql = ("SELECT p.kode AS kode, COUNT(k.id) AS jumlah FROM input1 p LEFT JOIN input2 k ON ST_Intersects(p.geometry, k.geometry) "
       "GROUP BY p.kode ORDER BY jumlah DESC, p.kode LIMIT 3")
b = processing.run("qgis:executesql", {"INPUT_DATASOURCES": [petak, kej], "INPUT_QUERY": sql, "INPUT_UID_FIELD": "",
                                       "INPUT_GEOMETRY_FIELD": "", "INPUT_GEOMETRY_TYPE": 1, "INPUT_GEOMETRY_CRS": None,
                                       "OUTPUT": "memory:"})["OUTPUT"]
hasil_b = [(f["kode"], f["jumlah"]) for f in b.getFeatures()]
print("Execute SQL           :", hasil_b)
print("SAMA:", hasil_a == hasil_b)

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Cara membaca baris 17. key=lambda t: (-t[1], t[0]) berarti: untuk tiap pasangan t, kunci urutnya adalah (-t[1], t[0]). t[1] adalah jumlah. Tanda minus membalik urutan, sehingga jumlah terbesar tampil lebih dulu. t[0] adalah kode petak, dipakai bila jumlahnya sama. Tanpa tanda minus, urutannya dari terkecil.

Hasil. Layar menampilkan:

Count points in polygon: [('P-04', 13.0), ('P-09', 12.0), ('P-12', 12.0)]
Execute SQL           : [('P-04', 13), ('P-09', 12), ('P-12', 12)]
SAMA: True

Angka cara pertama berbentuk 13.0 (desimal), cara kedua 13 (bulat). Python menganggap 13.0 == 13, jadi hasilnya tetap True. QGIS juga mencetak peringatan DeprecationWarning yang boleh diabaikan.

Bagian yang boleh diubah.

• Baris 16, nama kolom hitung. Ganti "jumlah" dengan nama lain, dan ganti juga di baris 17 dan 21.
• Baris 17 dan 21, angka 3. [:3] dan LIMIT 3 menentukan berapa petak teratas. Ganti keduanya bersama-sama.
• Baris 20-21, SQL. Ganti pertanyaannya. Tetap pakai input1, input2.

Jebakan.

• Urutan layer dan nama masukan. Nama input1 selalu layer pertama di daftar baris 22. Tukar urutannya, tukar pula maknanya.
• Pecah hasil seri. Dua petak sama-sama 12 kejadian. Tanpa p.kode di urutan SQL dan t[0] di lambda, urutan kedua cara bisa berbeda dan SAMA jadi False.
• Jenis bentuk hasil. Isian INPUT_GEOMETRY_TYPE bernilai 1, yang berarti tanpa bentuk. Daftar pilihan alat di QGIS 4.0.2: 0 otomatis, 1 tanpa bentuk, 2 titik, 3 garis, 4 poligon, lalu 5 sampai 7 untuk versi majemuk. Untuk hasil yang membawa bentuk, pilih nilai yang sesuai.

Latihan variasi. Tampilkan lima petak teratas. Contoh jawaban: ganti [:3] menjadi [:5] di baris 17 dan LIMIT 3 menjadi LIMIT 5 di baris 21. Layar menampilkan lima pasangan di tiap cara, lalu SAMA: True:

Count points in polygon: [('P-04', 13.0), ('P-09', 12.0), ('P-12', 12.0), ('P-08', 9.0), ('P-01', 7.0)]
Execute SQL           : [('P-04', 13), ('P-09', 12), ('P-12', 12), ('P-08', 9), ('P-01', 7)]
SAMA: True

Skrip 25.5: Rute terpendek dan tercepat di QGIS #

Berkas m2_03a_rute_qgis.py. Kepala Seksi ingin tahu jalur patroli dari Pos Jaga 1 ke Pos Jaga 4. Skrip ini mencarinya dengan alat bawaan QGIS, lalu menguji apa yang terjadi bila satu ruas jalan ditutup atau bila jalannya bercelah.

Tujuan. Menghitung rute terpendek (dalam meter) dan tercepat (dalam menit) antar dua pos, dalam lima skenario.

Bahan. Fasilitas.gpkg (titik "Pos Jaga 1" dan "Pos Jaga 4"), Jalan.gpkg (kolom status dan kecepatan_kmj), dan Jalan_celah.gpkg (jalan yang sengaja dibuat bercelah). Tidak butuh PostGIS. Tidak menulis berkas.

Konsep: jaringan dan biaya. Jaringan jalan itu seperti peta jalur kereta: ada simpul (persimpangan atau ujung jalan) dan ruas (jalan di antara simpul). Tiap ruas punya biaya (cost), yaitu panjangnya atau waktu tempuhnya. Rute terbaik adalah jalur dengan jumlah biaya terkecil. Dua istilah lain: strategi (apakah biaya yang diminimalkan panjang atau waktu) dan toleransi (jarak terjauh dua ujung jalan masih dianggap menyambung).

Skrip ini juga memakai tiga lambda. Gambar berikut menunjukkan bahwa lambda hanyalah cara singkat menulis fungsi.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 3: rute terpendek dan rute tercepat dengan alat bawaan QGIS (Processing, algoritma native). Area layanan: lihat m2_03c.
Penulis: Badar Mubarok Yogaswara. Pemakaian: python-qgis.bat m2_03a_rute_qgis.py <folder paket-m2>"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
from qgis.core import QgsVectorLayer, QgsPointXY, QgsProcessingFeedback

paket = sys.argv[1]
SRC = lambda b, l, f="": os.path.join(paket, b + ".gpkg|layername=" + l + (("|subset=" + f) if f else ""))
fas = QgsVectorLayer(SRC("Fasilitas", "Fasilitas"), "Fasilitas", "ogr")
titik = {f["nama"]: f.geometry().asPoint() for f in fas.getFeatures()}
awal, akhir = titik["Pos Jaga 1"], titik["Pos Jaga 4"]
fmt = lambda p: "%.3f,%.3f [EPSG:32749]" % (p.x(), p.y())


def rute(jaringan, strategi=0, toleransi=0.0, a=awal, b=akhir):
    lyr = QgsVectorLayer(jaringan, "jaringan", "ogr")
    hasil = processing.run("native:shortestpathpointtopoint", {
        "INPUT": lyr, "STRATEGY": strategi, "DIRECTION_FIELD": "", "VALUE_FORWARD": "", "VALUE_BACKWARD": "", "VALUE_BOTH": "",
        "DEFAULT_DIRECTION": 2, "SPEED_FIELD": "kecepatan_kmj", "DEFAULT_SPEED": 20, "TOLERANCE": toleransi,
        "START_POINT": fmt(a), "END_POINT": fmt(b), "OUTPUT": "memory:"})
    return hasil["TRAVEL_COST"], hasil["OUTPUT"]


jalan = SRC("Jalan", "Jalan")
jalan_buka = SRC("Jalan", "Jalan", "status<>'ditutup'")
celah = SRC("Jalan_celah", "Jalan_celah", "status<>'ditutup'")
tt = lambda s: "%.1f menit" % (s * 60)

c, o = rute(jalan)
print("1) terpendek, semua ruas      : %.1f m" % c)
c0, o0 = rute(jalan_buka)
print("2) terpendek, ruas tutup dibuang: %.1f m  (jumlah segmen keluaran: %d)" % (c0, o0.featureCount()))
cf, of = rute(jalan_buka, strategi=1)
print("3) tercepat (jam -> menit)    : %s" % tt(cf))
cg, og = rute(celah)
print("4) jalan ber-celah, toleransi 0: %.1f m" % cg)
cg5, og5 = rute(celah, toleransi=5.0)
print("5) jalan ber-celah, toleransi 5: %.1f m" % cg5)

print("Catatan: titik awal dan akhir diambil dari fasilitas yang berjarak 25 m di luar jalan; QGIS menempelkannya ke titik jalan terdekat.")

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas. Bagian awal sama dengan Skrip 25.4, jadi dijelaskan singkat.

Hasil. Layar menampilkan:

1) terpendek, semua ruas      : 2987.7 m
2) terpendek, ruas tutup dibuang: 3316.2 m  (jumlah segmen keluaran: 1)
3) tercepat (jam -> menit)    : 6.4 menit
4) jalan ber-celah, toleransi 0: 3351.3 m
5) jalan ber-celah, toleransi 5: 3316.2 m
Catatan: titik awal dan akhir diambil dari fasilitas yang berjarak 25 m di luar jalan; QGIS menempelkannya ke titik jalan terdekat.

Cara membacanya: menutup satu ruas menambah rute sekitar 328 meter (2987,7 ke 3316,2). Pada jalan bercelah, tanpa toleransi QGIS tidak melihat ujung yang bercelah itu menyambung, sehingga rutenya memutar (3351,3). Itu kesimpulan dari angkanya, celahnya sendiri tidak diperiksa langsung. [kemungkinan] Dengan toleransi 5 meter, celah dianggap tertutup dan hasilnya kembali sama dengan nomor 2.

Bagian yang boleh diubah.

• Baris 15, nama pos. Ganti dengan nama fasilitas lain: Pos Jaga 2, Pos Jaga 3, Kantor KPH, atau TPK.
• Baris 23, kecepatan bawaan. Ganti 20 (km per jam) untuk ruas tanpa isian kecepatan.
• Baris 29, saringan. Ganti status<>'ditutup' dengan syarat lain.
• Baris 37 dan 41, strategi dan toleransi. Ganti strategi=1 dan toleransi=5.0.

Jebakan.

• Nama pos salah. Bila nama tidak ada di kamus, baris 15 berhenti dengan KeyError: 'Pos Jaga 9'. Salin nama persis dari hasil baris 14.
• Satuan biaya. Untuk strategi tercepat, biayanya jam, bukan menit. Itu sebabnya ada lambda tt.
• Jalan bercelah. Tanpa toleransi, rute bisa memutar atau gagal. Selalu uji dengan toleransi kecil.
• Sistem koordinat tertulis tetap. EPSG:32749 di baris 16 harus cocok dengan data. Bila data Anda memakai sistem lain, ganti.
• Nilai bawaan dibaca sekali. a=awal di baris 19 dibaca saat fungsi dibuat. Mengganti awal sesudahnya tidak mengubah nilai bawaan.
• Lambda sulit dibaca pemula. Tiga lambda di sini boleh ditulis ulang sebagai def. Hasilnya sama.

Latihan variasi.

1. Cari rute dari Pos Jaga 1 ke TPK.
2. Tulis ulang tt sebagai fungsi def.

Contoh jawaban 1: ganti titik["Pos Jaga 4"] di baris 15 menjadi titik["TPK"]. Layar menampilkan lima baris dengan angka berbeda, antara lain 1) terpendek, semua ruas : 1521.4 m dan 3) tercepat (jam -> menit) : 5.1 menit.

Contoh jawaban 2: ganti baris 31 dengan dua baris. Hasil skrip tidak berubah (diuji dengan ketiga lambda ditulis ulang sebagai def).

def tt(s):
    return "%.1f menit" % (s * 60)

Skrip 25.6: Rute dan area layanan, QGIS lawan pgRouting #

Berkas m2_03c_silang_jaringan.py. Ini silang kedua. Rute yang sama dihitung dua kali: sekali di QGIS dan sekali di PostGIS dengan pgRouting, ekstensi PostgreSQL untuk rute. Bila angkanya sama, hitungan Anda layak dipercaya.

Tujuan. Membandingkan rute terpendek, rute tercepat, dan area layanan dari QGIS dengan pgRouting, memakai titik yang tepat di simpul.

Bahan. Tabel kph.jalan dan kph.fasilitas (hasil Skrip 25.1), pgRouting aktif, Jalan.gpkg, dan m2_03b_pgrouting.sql yang sudah dijalankan (membuat tabel kph.jalan_simpul dan kph.jalan_ruas). Perintah: python-qgis.bat m2_03c_silang_jaringan.py <folder paket-m2>.

Konsep singkat. pgRouting memakai simpul bernomor dan ruas yang menyebut simpul asal dan tujuannya. pgr_dijkstra mencari rute terpendek, dan pgr_drivingDistance mencari semua simpul yang terjangkau dalam batas jarak (area layanan). Skrip ini juga memakai %s, tanda tempat untuk nilai yang diisi oleh psycopg2 secara aman.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 3: membandingkan rute dan area layanan QGIS dengan pgRouting memakai titik awal dan akhir yang SAMA (tepat di simpul).
Penulis: Badar Mubarok Yogaswara. Syarat: m2_03b_pgrouting.sql sudah dijalankan; variabel PGHOST, PGPORT, PGUSER, PGDATABASE terisi.
Pemakaian: python-qgis.bat m2_03c_silang_jaringan.py <folder paket-m2>"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
import psycopg2
from qgis.core import QgsVectorLayer

paket = sys.argv[1]
con = psycopg2.connect(dbname=os.environ["PGDATABASE"])
cur = con.cursor()


def simpul(vid):
    cur.execute("SELECT ST_X(geom), ST_Y(geom) FROM kph.jalan_simpul WHERE id = %s", (vid,))
    x, y = cur.fetchone()
    return "%.3f,%.3f [EPSG:32749]" % (x, y)


jalan_buka = QgsVectorLayer(os.path.join(paket, "Jalan.gpkg|layername=Jalan|subset=status<>'ditutup'"), "j", "ogr")
a, b = simpul(5), simpul(24)                       # simpul terdekat Pos Jaga 1 dan Pos Jaga 4 (lihat keluaran 03b)


def qgis_rute(strategi):
    return processing.run("native:shortestpathpointtopoint", {
        "INPUT": jalan_buka, "STRATEGY": strategi, "DIRECTION_FIELD": "", "VALUE_FORWARD": "", "VALUE_BACKWARD": "", "VALUE_BOTH": "",
        "DEFAULT_DIRECTION": 2, "SPEED_FIELD": "kecepatan_kmj", "DEFAULT_SPEED": 20, "TOLERANCE": 0,
        "START_POINT": a, "END_POINT": b, "OUTPUT": "memory:"})["TRAVEL_COST"]


cur.execute("""SELECT MAX(agg_cost) FROM pgr_dijkstra('SELECT id, source, target, panjang_m AS cost FROM kph.jalan_ruas WHERE status <> ''ditutup''', 5, 24, false)""")
pg_m = cur.fetchone()[0]
cur.execute("""SELECT MAX(agg_cost) FROM pgr_dijkstra('SELECT id, source, target, waktu_dtk AS cost FROM kph.jalan_ruas WHERE status <> ''ditutup''', 5, 24, false)""")
pg_dtk = cur.fetchone()[0]
q_m, q_jam = qgis_rute(0), qgis_rute(1)
print("terpendek : QGIS %.2f m | pgRouting %.2f m | selisih %.4f" % (q_m, pg_m, abs(q_m - pg_m)))
print("tercepat  : QGIS %.3f menit | pgRouting %.3f menit" % (q_jam * 60, pg_dtk / 60))

# area layanan: simpul persimpangan yang terjangkau
cur.execute("SELECT ST_X(geom), ST_Y(geom), id FROM kph.jalan_simpul")
persimpangan = {(round(x, 1), round(y, 1)): i for x, y, i in cur.fetchall()}
cur.execute("SELECT ST_X(geom), ST_Y(geom) FROM kph.jalan_simpul WHERE id = 11")
kx, ky = cur.fetchone()
for batas in (500, 1000):
    h = processing.run("native:serviceareafrompoint", {
        "INPUT": jalan_buka, "STRATEGY": 0, "DIRECTION_FIELD": "", "VALUE_FORWARD": "", "VALUE_BACKWARD": "", "VALUE_BOTH": "",
        "DEFAULT_DIRECTION": 2, "SPEED_FIELD": "", "DEFAULT_SPEED": 50, "TOLERANCE": 0, "START_POINT": "%.3f,%.3f [EPSG:32749]" % (kx, ky),
        "TRAVEL_COST": batas, "INCLUDE_BOUNDS": True, "OUTPUT": "memory:"})["OUTPUT"]
    dalam = set()
    for f in h.getFeatures():
        if f["type"] == "within":
            g = f.geometry()
            for p in (g.asMultiPoint() if g.isMultipart() else [g.asPoint()]):
                k = (round(p.x(), 1), round(p.y(), 1))
                if k in persimpangan:
                    dalam.add(persimpangan[k])
    cur.execute("""SELECT COUNT(*) FROM pgr_drivingDistance('SELECT id, source, target, panjang_m AS cost FROM kph.jalan_ruas WHERE status <> ''ditutup''', 11, %s::float, false)""", (batas,))
    print("area layanan %4d m: persimpangan terjangkau QGIS %d | pgRouting %d" % (batas, len(dalam), cur.fetchone()[0]))

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Layar menampilkan:

terpendek : QGIS 3296.93 m | pgRouting 3296.93 m | selisih 0.0000
tercepat  : QGIS 6.373 menit | pgRouting 6.373 menit
area layanan  500 m: persimpangan terjangkau QGIS 3 | pgRouting 3
area layanan 1000 m: persimpangan terjangkau QGIS 10 | pgRouting 10

Dua sistem sepakat sampai dua desimal. Angka 3296,93 meter lebih kecil dari 3316,2 meter di Skrip 25.5, karena sekarang rute dihitung dari simpul ke simpul, bukan dari pos jaga yang berjarak dari jalan. [kemungkinan]

Bagian yang boleh diubah.

• Baris 25, nomor simpul. 5 dan 24 bisa diganti. Angka yang sama juga tertulis di baris 35 dan 37, jadi ubah di tiga tempat.
• Baris 46, nomor simpul kantor. 11 untuk titik pusat area layanan.
• Baris 48, batas jarak. Ganti (500, 1000) dengan batas lain, dalam meter.
• Baris 31 dan 50, kecepatan bawaan. Isi DEFAULT_SPEED bila ruas tidak punya data kecepatan.

Jebakan.

• 03b belum dijalankan. Tabel kph.jalan_simpul belum ada, dan baris 19 berhenti dengan psycopg2.errors.UndefinedTable: relasi « kph.jalan_simpul » tidak ada (pesan mengikuti bahasa server; diuji pada basis data tanpa tabel itu).
• Nomor simpul kode keras. Angka 5, 24, dan 11 berasal dari keluaran m2_03b_pgrouting.sql pada data ini. Bila data jalan berubah, nomor bisa bergeser. Baca keluaran 03b untuk mencarinya.
• Kutip ganda dua kali. ''ditutup'' di dalam SQL bersarang mudah salah ketik.
• Koneksi tidak ditutup. con.close() tidak ada. Untuk skrip pendek tidak masalah.
• Pencocokan titik dengan pembulatan. Titik dicocokkan lewat koordinat yang dibulatkan satu desimal. Bila data berubah, pencocokan bisa meleset. [kemungkinan]

Latihan variasi. Bandingkan rute dari Pos Jaga 1 ke Pos Jaga 2. Pos Jaga 2 berada di simpul 21 menurut keluaran 03b. Contoh jawaban: ganti simpul(24) menjadi simpul(21) di baris 25, dan 5, 24, false menjadi 5, 21, false di baris 35 dan 37. Layar menampilkan terpendek : QGIS 2739.99 m | pgRouting 2739.99 m | selisih 0.0000 dan tercepat : QGIS 7.338 menit | pgRouting 7.338 menit.

Skrip 25.7: Hidrologi dari peta ketinggian #

Berkas m2_04a_hidrologi.py. Dari satu peta ketinggian (DEM), skrip ini mencari ke mana air mengalir, di mana sungainya, dan daerah mana yang airnya bermuara ke satu titik. Hasilnya berguna untuk menilai petak mana yang dekat sungai.

Tujuan. Mengisi lubang pada DEM, menghitung arah dan akumulasi aliran, mengekstrak sungai, mencari titik keluar, dan membentuk daerah aliran sungai (DAS).

Bahan. DEM_10m.tif (200 x 200 sel, sel 10 meter) dan folder hasil. Butuh GRASS (lewat _inisialisasi). Tidak butuh PostGIS. Waktu jalan pada uji: sekitar 22 detik. Hasilnya disimpan di hasil.

Konsep: air mengalir ke bawah. Bayangkan DEM sebagai maket pegunungan dari tanah liat. Hujan yang jatuh di tiap titik mengalir ke tetangga yang paling rendah. Lima istilah baru:

• DEM (Digital Elevation Model): peta ketinggian berbentuk grid. Tiap sel berisi tinggi.
• Lubang (sink): sel yang lebih rendah daripada semua tetangganya. Air terjebak di sana. Biasanya itu kesalahan data, jadi diisi.
• Akumulasi aliran: untuk tiap sel, jumlah sel di hulu yang airnya lewat sini. Sel dengan angka besar adalah calon sungai.
• Ambang (threshold): batas akumulasi untuk menyebut sebuah sel "sungai". Di bab ini 500 sel, yaitu 500 x 100 m2 = 5 hektare.
• DAS dan titik keluar (outlet): DAS adalah seluruh daerah yang airnya bermuara ke satu titik keluar.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 4: hidrologi dari DEM: isi lubang (fill sinks), arah dan akumulasi aliran, sungai, dan DAS di satu titik keluar (outlet).
Penulis: Badar Mubarok Yogaswara. Pemakaian: python-qgis.bat m2_04a_hidrologi.py <folder paket-m2>
Keluaran disimpan di <paket-m2>/hasil. Ambang sungai: 500 sel (500 x 100 m2 = 5 ha)."""
import os
import sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
from osgeo import gdal

gdal.UseExceptions()
paket = sys.argv[1]
out = os.path.join(paket, "hasil")
os.makedirs(out, exist_ok=True)
dem = os.path.join(paket, "DEM_10m.tif")
P = lambda n: os.path.join(out, n)
BAKU = {"GRASS_REGION_PARAMETER": None, "GRASS_REGION_CELLSIZE_PARAMETER": 0}


def baca(path):
    ds = gdal.Open(path)
    return ds.GetRasterBand(1).ReadAsArray().astype("float64"), ds.GetGeoTransform()


def jumlah_lubang(a):
    """jumlah sel dalam (bukan tepi) yang lebih rendah daripada semua 8 tetangganya"""
    c = a[1:-1, 1:-1]
    lebih_tinggi = np.ones_like(c, dtype=bool)
    for dy in (-1, 0, 1):
        for dx in (-1, 0, 1):
            if dy == 0 and dx == 0:
                continue
            lebih_tinggi &= a[1 + dy:a.shape[0] - 1 + dy, 1 + dx:a.shape[1] - 1 + dx] > c
    return int(lebih_tinggi.sum())


z0, gt = baca(dem)
print("lubang (sink) pada DEM asli            :", jumlah_lubang(z0))

# 1) isi lubang (Wang & Liu), kemiringan minimum 0,1 derajat
processing.run("native:fillsinkswangliu", {"INPUT": dem, "BAND": 1, "MIN_SLOPE": 0.1, "OUTPUT_FILLED_DEM": P("DEM_terisi.tif"),
                                           "OUTPUT_FLOW_DIRECTIONS": P("Arah_WL.tif"), "OUTPUT_WATERSHED_BASINS": P("Basin_WL.tif")})
z1, _ = baca(P("DEM_terisi.tif"))
print("lubang sesudah diisi                   :", jumlah_lubang(z1))
print("sel yang naik / kenaikan terbesar (m)  : %d / %.2f" % (int((z1 > z0 + 1e-6).sum()), float((z1 - z0).max())))

# 2) akumulasi aliran, arah, dan cekungan dengan GRASS r.watershed (aliran tunggal D8)
processing.run("grass:r.watershed", dict(BAKU, elevation=P("DEM_terisi.tif"), threshold=500, convergence=5, memory=300, **{"-s": True, "-a": True},
               accumulation=P("Akumulasi.tif"), drainage=P("Arah.tif"), basin=P("Cekungan.tif")))
acc, _ = baca(P("Akumulasi.tif"))
print("akumulasi maksimum (sel) / luas (ha)   : %d / %.1f" % (abs(acc).max(), abs(acc).max() * 100 / 1e4))

# 3) sungai: sel dengan akumulasi >= 500 sel (5 ha)
processing.run("grass:r.stream.extract", dict(BAKU, elevation=P("DEM_terisi.tif"), accumulation=P("Akumulasi.tif"), threshold=500, mexp=0,
               stream_length=0, memory=300, stream_raster=P("Sungai.tif"), stream_vector=P("sungai_mentah.gpkg"), direction=P("Arah_stream.tif"),
               GRASS_OUTPUT_TYPE_PARAMETER=2, GRASS_VECTOR_EXPORT_NOCAT=False))
if os.path.exists(P("Sungai.gpkg")):
    os.remove(P("Sungai.gpkg"))
gdal.VectorTranslate(P("Sungai.gpkg"), P("sungai_mentah.gpkg"), options=gdal.VectorTranslateOptions(
    format="GPKG", layerName="Sungai", selectFields=["cat"], geometryType="LINESTRING"))
from osgeo import ogr
ds = ogr.Open(P("Sungai.gpkg"))
lyr = ds.GetLayer(0)
panjang = sum(f.GetGeometryRef().Length() for f in lyr)
print("sungai: %d segmen, total panjang %.2f km" % (lyr.GetFeatureCount(), panjang / 1000.0))
ds = None


# 4) titik keluar (outlet): sel tepi dengan akumulasi terbesar, lalu DAS-nya
acc, gt = baca(P("Akumulasi.tif"))
tepi = np.zeros_like(acc, dtype=bool)
tepi[0, :] = tepi[-1, :] = tepi[:, 0] = tepi[:, -1] = True
r, c = np.unravel_index(np.argmax(np.where(tepi, acc, -1)), acc.shape)
ox, oy = gt[0] + (c + 0.5) * gt[1], gt[3] + (r + 0.5) * gt[5]
print("outlet: baris %d kolom %d, X=%.1f Y=%.1f, akumulasi %d sel" % (r, c, ox, oy, acc[r, c]))
processing.run("grass:r.water.outlet", dict(BAKU, input=P("Arah.tif"), coordinates="%.1f,%.1f [EPSG:32749]" % (ox, oy), output=P("DAS_outlet.tif")))
das, _ = baca(P("DAS_outlet.tif"))
n_das = int((das == 1).sum())
print("DAS di outlet: %d sel = %.2f ha (%.1f%% wilayah)" % (n_das, n_das * 100 / 1e4, 100.0 * n_das / das.size))
mentah = processing.run("gdal:polygonize", {"INPUT": P("DAS_outlet.tif"), "BAND": 1, "FIELD": "nilai", "EIGHT_CONNECTEDNESS": False,
                                            "EXTRA": "", "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
dalam = processing.run("native:extractbyexpression", {"INPUT": mentah, "EXPRESSION": '"nilai" = 1', "OUTPUT": "memory:"})["OUTPUT"]
if os.path.exists(P("DAS_utama.gpkg")):
    os.remove(P("DAS_utama.gpkg"))
processing.run("native:dissolve", {"INPUT": dalam, "FIELD": [], "SEPARATE_DISJOINT": False, "OUTPUT": P("DAS_utama.gpkg")})
print("DAS poligon ditulis:", os.path.exists(P("DAS_utama.gpkg")))

Skrip ini 88 baris. Penjelasannya dipecah menjadi enam blok.

Blok 1, baris 1-19: persiapan.

Blok 2, baris 22-36: dua fungsi pembantu.

Blok 3, baris 39-47: mengisi lubang.

Blok 4, baris 49-53: arah dan akumulasi aliran.

Blok 5, baris 55-68: sungai.

Blok 6, baris 71-88: titik keluar dan DAS.

Hasil. Layar menampilkan:

lubang (sink) pada DEM asli            : 167
lubang sesudah diisi                   : 0
sel yang naik / kenaikan terbesar (m)  : 1994 / 13.67
akumulasi maksimum (sel) / luas (ha)   : 25581 / 255.8
sungai: 30 segmen, total panjang 7.75 km
outlet: baris 199 kolom 49, X=312495.0 Y=9996005.0, akumulasi 25581 sel
DAS di outlet: 26211 sel = 262.11 ha (65.5% wilayah)
DAS poligon ditulis: True

Cara membaca: ada 167 lubang di DEM asli, semuanya terisi. Sungai terdiri atas 30 segmen sepanjang 7,75 km. DAS yang bermuara di titik keluar menutupi 65,5 persen wilayah (262 ha dari 400 ha). Folder hasil kini berisi DEM_terisi.tif, Akumulasi.tif, Sungai.gpkg, DAS_utama.gpkg, dan lainnya. QGIS juga mencetak sedikit teks proses GRASS yang boleh diabaikan.

Bagian yang boleh diubah.

• Baris 56, ambang sungai. threshold=500: makin kecil, sungai makin rapat. Ganti ke 200 untuk sungai lebih rapat.
• Baris 43, kemiringan minimum. MIN_SLOPE: 0.1. Ganti bila air di tanah datar tidak mengalir sebagaimana mestinya.
• Baris 50, ambang cekungan. Isian threshold=500 di sini terpisah dari ambang sungai di baris 56.
• Baris 53 dan 81, luas sel. Angka 100 (m2) benar hanya untuk sel 10 meter. Untuk DEM lain, ganti.
• Baris 17, DEM. Ganti DEM_10m.tif dengan DEM Anda.

Jebakan.

• GRASS tidak terdeteksi. Bila folder GRASS tidak dikenali, baris 50 berhenti dengan QgsProcessingException: GRASS GIS folder is not configured. Please configure it before running GRASS GIS algorithms. (diuji dengan mengganti nama folder GRASS). Cek baris 19 _inisialisasi.py.
• Luas sel tertulis tetap. Baris 53 dan 81 memakai 100. Dengan DEM berukuran sel lain, hektarenya salah tanpa peringatan.
• Titik keluar ditentukan otomatis. Skrip memilih sel tepi dengan akumulasi terbesar. Itu bukan titik pilihan Anda.
• Hasil lama dihapus. Baris 59-60 dan 85-86 menghapus berkas lama. Jangan taruh hasil penting dengan nama sama di folder hasil.
• Hasil bergantung versi. Angka di atas berlaku untuk QGIS 4.0.2 dengan GRASS 8.4. Versi lain bisa sedikit berbeda. [CEK]

Latihan variasi. Buat sungai lebih rapat dengan ambang 200 sel. Contoh jawaban: ganti threshold=500 menjadi threshold=200 hanya di baris 56 (baris r.stream.extract). Baris sungai: menampilkan sungai: 105 segmen, total panjang 15.47 km. DAS dan titik keluar tidak berubah, karena keduanya tidak memakai ambang sungai.

Skrip 25.8: Peta curah hujan dari 24 stasiun #

Berkas m2_05a_interpolasi.py. Kepala Seksi punya 30 stasiun hujan, tetapi hujan jatuh di mana-mana, bukan hanya di stasiun. Skrip ini menebak curah hujan tahunan di seluruh wilayah dari titik-titik stasiun, memakai sembilan metode, lalu memilih yang paling jitu. Skrip ini panjang (167 baris), jadi dipecah menjadi enam blok.

Tujuan. Membuat peta curah hujan tahunan dengan banyak metode interpolasi, menilai tiap metode dengan stasiun yang disembunyikan, dan menyimpan peta metode terbaik sebagai CH_tahunan.tif.

Bahan. Stasiun_Hujan.gpkg (30 stasiun; kolom elevasi_m, ch_tahunan_mm, dan peran yang bernilai latih atau uji) dan DEM_10m.tif. Butuh GRASS, NumPy, dan SciPy. Tidak butuh PostGIS. Waktu jalan pada uji: sekitar 23 detik.

Konsep: menebak di tempat yang tidak diukur. Bayangkan Anda hanya punya termometer di beberapa desa dan ingin menebak suhu di desa lain. Itulah interpolasi: menebak nilai di tempat tak terukur dari titik terukur. Lima istilah baru:

• IDW (inverse distance weighting): rata-rata berbobot, di mana stasiun yang lebih dekat "bersuara" lebih keras. Pangkat p menentukan seberapa cepat suara itu mengecil dengan jarak.
• TIN: jaring segitiga yang menyambung stasiun. Di dalam tiap segitiga, nilai dihitung dari tiga sudutnya. Di luar jaring, tidak ada nilai.
• Data latih dan data uji: 24 stasiun dipakai membuat peta (latih). Enam stasiun lain disembunyikan, lalu dipakai menguji tebakan (uji). Seperti ujian: soalnya belum pernah dilihat.
• RMSE: ukuran galat. Rata-rata selisih tebakan dan kenyataan, dengan selisih besar dihitung lebih berat. Makin kecil makin baik. MAE adalah rata-rata selisih mutlak, dan bias adalah rata-rata selisih bertanda (positif berarti tebakan cenderung kelebihan).
• Tinggalkan-satu (leave-one-out, LOO): uji dengan menyembunyikan satu stasiun bergantian, sebanyak stasiun yang ada.

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 5: interpolasi curah hujan tahunan dari 24 stasiun latih, dinilai dengan 6 stasiun uji.
Metode: IDW (p = 1, 2, 3), TIN linear, TIN Clough-Tocher, B-spline (GRASS), RST (GRASS), kriging biasa (NumPy), dan regresi + residu IDW (NumPy).
Penulis: Badar Mubarok Yogaswara. Pemakaian: python-qgis.bat m2_05a_interpolasi.py <folder paket-m2>
Keluaran: <paket-m2>/hasil/CH_tahunan.tif (metode terbaik) dan tabel galat (RMSE) di layar."""
import os
import sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
from osgeo import gdal, ogr
from scipy.optimize import curve_fit

gdal.UseExceptions()
ogr.UseExceptions()
paket = sys.argv[1]
out = os.path.join(paket, "hasil")
os.makedirs(out, exist_ok=True)
P = lambda n: os.path.join(out, n)
X0, Y0, L, RES = 312000.0, 9996000.0, 2000.0, 10.0
EXT = "%f,%f,%f,%f [EPSG:32749]" % (X0, X0 + L, Y0, Y0 + L)

# ---- data stasiun: latih dan uji
ds = ogr.Open(os.path.join(paket, "Stasiun_Hujan.gpkg"))
lyr = ds.GetLayer(0)
st = [(f.GetGeometryRef().GetX(), f.GetGeometryRef().GetY(), f["elevasi_m"], f["ch_tahunan_mm"], f["peran"]) for f in lyr]
latih = np.array([s[:4] for s in st if s[4] == "latih"])
uji = np.array([s[:4] for s in st if s[4] == "uji"])
print("stasiun latih: %d, uji: %d" % (len(latih), len(uji)))
processing.run("native:extractbyexpression", {"INPUT": os.path.join(paket, "Stasiun_Hujan.gpkg|layername=Stasiun_Hujan"),
                                              "EXPRESSION": "\"peran\" = 'latih'", "OUTPUT": P("stasiun_latih.gpkg")})
from qgis.core import QgsVectorLayer
fi = QgsVectorLayer(P("stasiun_latih.gpkg"), "latih", "ogr").fields().indexOf("ch_tahunan_mm")      # indeks kolom menurut QGIS (kolom fid ikut dihitung)
data_in = "%s::~::0::~::%d::~::0" % (P("stasiun_latih.gpkg"), fi)


def contoh(raster_path, titik):
    r = gdal.Open(raster_path)
    gt = r.GetGeoTransform()
    a = r.GetRasterBand(1).ReadAsArray()
    nd = r.GetRasterBand(1).GetNoDataValue()
    h = np.array([a[int((gt[3] - y) / -gt[5]), int((x - gt[0]) / gt[1])] for x, y in titik[:, :2]], dtype="float64")
    if nd is not None:
        h[h == nd] = np.nan                           # di luar segitiga TIN: tidak ada nilai
    return h


hasil = {}
for p in (1, 2, 3):
    processing.run("qgis:idwinterpolation", {"INTERPOLATION_DATA": data_in, "DISTANCE_COEFFICIENT": p, "EXTENT": EXT, "PIXEL_SIZE": RES,
                                             "OUTPUT": P("CH_idw_p%d.tif" % p)})
    hasil["IDW p=%d" % p] = P("CH_idw_p%d.tif" % p)
for nama, m in (("TIN linear", 0), ("TIN Clough-Tocher", 1)):
    fn = P("CH_tin_%d.tif" % m)
    processing.run("qgis:tininterpolation", {"INTERPOLATION_DATA": data_in, "METHOD": m, "EXTENT": EXT, "PIXEL_SIZE": RES, "OUTPUT": fn})
    hasil[nama] = fn
BAKU = {"GRASS_REGION_PARAMETER": EXT, "GRASS_REGION_CELLSIZE_PARAMETER": RES}
processing.run("grass:v.surf.bspline", dict(BAKU, input=P("stasiun_latih.gpkg"), column="ch_tahunan_mm", ew_step=500, ns_step=500, method=1,
               lambda_i=0.01, solver=0, maxit=10000, error=1e-6, memory=300, raster_output=P("CH_bspline.tif")))
hasil["B-spline bikubik"] = P("CH_bspline.tif")
processing.run("grass:v.surf.rst", dict(BAKU, input=P("stasiun_latih.gpkg"), zcolumn="ch_tahunan_mm", tension=40, smooth=0.5, segmax=40, npmin=150,
               elevation=P("CH_rst.tif")))
hasil["RST (spline tegangan)"] = P("CH_rst.tif")

# ---- baris dan kolom sel pusat grid, untuk metode NumPy
n = int(L / RES)
gx = X0 + (np.arange(n) + 0.5) * RES
gy = Y0 + L - (np.arange(n) + 0.5) * RES
GX, GY = np.meshgrid(gx, gy)
dem = gdal.Open(os.path.join(paket, "DEM_10m.tif")).ReadAsArray().astype("float64")


def idw_np(xy, v, tx, ty, p=2):
    d = np.hypot(tx[..., None] - xy[:, 0], ty[..., None] - xy[:, 1])
    d = np.maximum(d, 1e-9)
    w = 1.0 / d ** p
    return (w * v).sum(-1) / w.sum(-1)


# kriging biasa dengan semivariogram eksponensial yang disesuaikan ke data latih
xy = latih[:, :2]
v = latih[:, 3]
ii, jj = np.triu_indices(len(xy), 1)
h = np.hypot(*(xy[ii] - xy[jj]).T)
gam = 0.5 * (v[ii] - v[jj]) ** 2
tepi = np.linspace(0, h.max() * 0.6, 8)
hc, gc = [], []
for a, b in zip(tepi[:-1], tepi[1:]):
    m = (h >= a) & (h < b)
    if m.sum() >= 5:
        hc.append(h[m].mean())
        gc.append(gam[m].mean())
model = lambda hh, c0, c, a: c0 + c * (1 - np.exp(-3 * hh / a))
(c0, c, a), _ = curve_fit(model, hc, gc, p0=[gc[0], max(gc) - gc[0], hc[-1]], bounds=([0, 1, 100], [1e7, 1e7, 5000]))
print("semivariogram eksponensial: nugget=%.0f, sill parsial=%.0f, jangkauan=%.0f m" % (c0, c, a))


def krig(tx, ty):
    N = len(xy)
    D = np.hypot(xy[:, None, 0] - xy[None, :, 0], xy[:, None, 1] - xy[None, :, 1])
    A = np.ones((N + 1, N + 1))
    A[:N, :N] = model(D, 0, c, a) + np.where(D == 0, 0, c0)
    A[N, N] = 0
    Ainv = np.linalg.inv(A)
    shp = tx.shape
    d0 = np.hypot(tx.ravel()[:, None] - xy[:, 0], ty.ravel()[:, None] - xy[:, 1])
    B = np.ones((d0.shape[0], N + 1))
    B[:, :N] = model(d0, 0, c, a) + c0
    lam = B @ Ainv
    return (lam[:, :N] @ v).reshape(shp)


# regresi pada elevasi (dan koordinat X), lalu residunya dengan IDW p=2
def elev_di(x, y):
    return dem[np.clip(((Y0 + L - y) / RES).astype(int), 0, n - 1), np.clip(((x - X0) / RES).astype(int), 0, n - 1)]


Xr = np.column_stack([np.ones(len(latih)), latih[:, 2], latih[:, 0] - X0])
beta, *_ = np.linalg.lstsq(Xr, v, rcond=None)
res = v - Xr @ beta
print("regresi CH = %.1f + %.2f * elevasi + %.4f * (X - X0)   (R2 = %.3f)" % (beta[0], beta[1], beta[2], 1 - res.var() / v.var()))


def regresi(tx, ty):
    return beta[0] + beta[1] * elev_di(tx, ty) + beta[2] * (tx - X0) + idw_np(xy, res, tx, ty, 2)


np_metode = {"Kriging biasa (NumPy)": krig, "Regresi elevasi + residu IDW (NumPy)": regresi}
tabel = []
for nama, fn in hasil.items():
    pred = contoh(fn, uji)
    tabel.append((nama, pred))
for nama, fn in np_metode.items():
    pred = fn(uji[:, 0], uji[:, 1])
    tabel.append((nama, pred))
print("\nGalat pada 6 stasiun uji (mm/tahun), dihitung dari prediksi - pengamatan")
print("%-40s %8s %8s %8s %4s" % ("metode", "RMSE", "MAE", "bias", "n"))
skor = {}
for nama, pred in tabel:
    e = (pred - uji[:, 3])
    e = e[~np.isnan(e)]
    skor[nama] = float(np.sqrt((e ** 2).mean())) if len(e) == len(uji) else 1e9      # metode yang tak mencakup semua titik uji tidak dipilih
    print("%-40s %8.1f %8.1f %8.1f %4d" % (nama, np.sqrt((e ** 2).mean()), np.abs(e).mean(), e.mean(), len(e)))
print("terbaik menurut RMSE:", min(skor, key=skor.get))

# simpan peta terbaik dari metode NumPy bila memang terbaik; jika tidak, salin raster QGIS terbaik
terbaik = min(skor, key=skor.get)
if terbaik in np_metode:
    arr = np_metode[terbaik](GX, GY).astype("float32")
    d = gdal.GetDriverByName("GTiff").Create(P("CH_tahunan.tif"), n, n, 1, gdal.GDT_Float32, ["COMPRESS=DEFLATE"])
    d.SetGeoTransform((X0, RES, 0, Y0 + L, 0, -RES))
    d.SetProjection(gdal.Open(os.path.join(paket, "DEM_10m.tif")).GetProjection())
    d.GetRasterBand(1).WriteArray(arr)
    d = None
else:
    gdal.Translate(P("CH_tahunan.tif"), hasil[terbaik], creationOptions=["COMPRESS=DEFLATE"])
print("CH_tahunan.tif dibuat dari:", terbaik)

# validasi silang tinggalkan-satu (LOO) untuk pangkat IDW pada 24 stasiun latih
print("\nTinggalkan-satu (LOO) pada 24 stasiun latih, RMSE IDW menurut pangkat:")
for p in (0.5, 1, 2, 3, 4):
    e = []
    for k in range(len(xy)):
        m = np.arange(len(xy)) != k
        e.append(idw_np(xy[m], v[m], np.array([xy[k, 0]]), np.array([xy[k, 1]]), p)[0] - v[k])
    print("  p=%.1f  RMSE=%.1f" % (p, np.sqrt(np.mean(np.square(e)))))

Blok A, baris 1-35: persiapan.

Blok B, baris 38-46: fungsi pembaca contoh.

Blok C, baris 49-64: metode QGIS dan GRASS.

Blok D, baris 66-126: metode NumPy. Tiga metode berikutnya ditulis sendiri dengan NumPy.

Blok E, baris 129-145: menilai semua metode dengan data uji.

Blok F, baris 147-167: simpan peta terbaik dan uji tinggalkan-satu.

Hasil. Layar menampilkan (setelah beberapa peringatan GRASS yang boleh diabaikan):

stasiun latih: 24, uji: 6
semivariogram eksponensial: nugget=0, sill parsial=66870, jangkauan=5000 m
regresi CH = 1944.2 + 3.04 * elevasi + 0.1371 * (X - X0)   (R2 = 0.912)

Galat pada 6 stasiun uji (mm/tahun), dihitung dari prediksi - pengamatan
metode                                       RMSE      MAE     bias    n
IDW p=1                                     183.4    175.2    150.1    6
IDW p=2                                     126.5    119.7     97.4    6
IDW p=3                                      94.6     83.8     61.0    6
TIN linear                                   55.9     40.5     38.6    2
TIN Clough-Tocher                            98.3     74.6    -45.6    3
B-spline bikubik                            100.7     84.1      9.5    6
RST (spline tegangan)                        89.9     76.2     41.0    6
Kriging biasa (NumPy)                        83.7     68.4     33.5    6
Regresi elevasi + residu IDW (NumPy)         42.2     33.3      9.6    6
terbaik menurut RMSE: Regresi elevasi + residu IDW (NumPy)
CH_tahunan.tif dibuat dari: Regresi elevasi + residu IDW (NumPy)

Tinggalkan-satu (LOO) pada 24 stasiun latih, RMSE IDW menurut pangkat:
  p=0.5  RMSE=169.0
  p=1.0  RMSE=149.4
  p=2.0  RMSE=115.3
  p=3.0  RMSE=96.0
  p=4.0  RMSE=90.1

Cara membaca: kolom n menunjukkan berapa dari enam titik uji yang tercakup. TIN linear hanya mencakup 2 titik, sehingga RMSE-nya (55,9) tampak bagus tetapi tidak adil dibandingkan. Skrip menolaknya. Pemenangnya regresi dengan RMSE 42,2 mm per tahun. Regresinya berbunyi: tiap kenaikan elevasi 100 meter menambah sekitar 304 mm hujan. Pada uji tinggalkan-satu, RMSE IDW terus turun sampai pangkat 4. Peta CH_tahunan.tif berukuran 200 x 200 sel, bernilai antara sekitar 2170 dan 3053 mm per tahun. Semua ini data sintetis.

Bagian yang boleh diubah.

• Baris 21, cakupan dan ukuran sel. X0, Y0, L, RES. Ubah RES ke 20 untuk sel lebih kasar dan lebih cepat.
• Baris 50, pangkat IDW. Ganti (1, 2, 3). Judul di tabel ikut berubah lewat baris 53.
• Baris 59-60, B-spline. ew_step, ns_step, dan lambda_i mengatur kehalusan.
• Baris 62-63, RST. tension dan smooth.
• Baris 95, batas jangkauan. Angka 5000 pada bounds membatasi jangkauan semivariogram.
• Baris 162, pangkat uji LOO. Ganti daftar (0.5, 1, 2, 3, 4).

Jebakan.

• Indeks kolom pada data masukan. Baris 34-35 mencari nomor kolom lewat QGIS. Kolom fid ikut dihitung, jadi nomornya 4, bukan 3. Salah nomor berarti yang diinterpolasi kolom lain.
• TIN tidak mencakup semua titik. Di luar segitiga, TIN tidak punya nilai, sehingga n kurang dari 6. Skrip membuang metode itu dengan skor 1e9, jadi pemenangnya tidak bisa TIN.
• Jangkauan semivariogram menempel di batas. Hasil jangkauan=5000 m sama dengan batas atas di baris 95. Itu tanda data tidak cukup untuk menentukan jangkauan dengan tegas. [kemungkinan]
• Matriks kriging bisa singular. Bila dua stasiun berada di titik yang sama, np.linalg.inv gagal. Pada data ini tidak terjadi.
• Butuh GRASS. Tanpa GRASS yang terdeteksi, baris 59 gagal (lihat Skrip 25.7).
• Hasil lama ditimpa. Berkas bernama sama di folder hasil akan diganti.
• Angka 24 dan 6 hanya di teks. Hitungan tidak memakainya. Jumlah latih dan uji ditentukan kolom peran di data.

Latihan variasi. Uji pangkat IDW 2, 4, dan 6 pada tinggalkan-satu. Contoh jawaban: ganti baris 162 menjadi for p in (2, 4, 6):. Layar menampilkan tiga baris terakhir:

  p=2.0  RMSE=115.3
  p=4.0  RMSE=90.1
  p=6.0  RMSE=92.4

Jadi pada data ini, RMSE terendah ada di sekitar pangkat 4, dan naik lagi di pangkat 6.

Skrip 25.9: Peta kesesuaian lahan dengan tumpang tindih berbobot #

Berkas m2_06a_kesesuaian.py. Kepala Seksi bertanya: lahan mana yang paling cocok untuk "jenis pohon contoh"? Skrip ini menilai tiap sel peta dari empat kriteria, menimbangnya, menerapkan larangan, lalu merangkum luasnya per kelas dan per petak. Kriteria dan ambangnya fiktif untuk latihan, bukan standar resmi (kata penulis di berkas).

Tujuan. Membuat peta kelas kesesuaian (S1 sangat sesuai sampai N tidak sesuai, plus kelas terlarang), memeriksanya dengan NumPy, menguji kepekaan bobot, dan merangkum per petak.

Bahan. DEM_10m.tif, Tanah.gpkg (kolom jenis_tanah), Jalan.gpkg, Petak.gpkg, dan dua keluaran bab sebelumnya di folder hasil: Sungai.gpkg (Skrip 25.7) serta CH_tahunan.tif (Skrip 25.8). Jalankan 25.7 dan 25.8 dulu. Tidak butuh PostGIS. Waktu jalan pada uji: sekitar 23 detik.

Konsep: skor, bobot, larangan. Bayangkan juri lomba dengan empat penilaian. Tiap sel peta mendapat skor 1 sampai 4 untuk tiap kriteria. Tiap kriteria punya bobot yang menyatakan seberapa penting. Skor akhir adalah jumlah skor kali bobot. Sebagian sel punya kendala (larangan) dan langsung dinyatakan terlarang, berapa pun skornya. Bobotnya tidak ditebak: dihitung dengan AHP (Analytic Hierarchy Process), metode yang menimbang kriteria lewat perbandingan berpasangan ("lereng lebih penting dari tanah, seberapa?"). Lima istilah baru: skor, bobot, kendala, AHP, dan reklasifikasi (mengubah nilai ke kelas, misalnya lereng 0 sampai 8 persen menjadi skor 4).

Skrip utuh:

# -*- coding: utf-8 -*-
"""M2 Bab 6: kesesuaian lahan multi-kriteria (tumpang tindih berbobot) untuk "jenis pohon contoh" di KPH Contoh.
Kriteria dan ambang di bawah ini FIKTIF untuk latihan; bukan standar kesesuaian lahan resmi.
Penulis: Badar Mubarok Yogaswara. Pemakaian: python-qgis.bat m2_06a_kesesuaian.py <folder paket-m2>
Prasyarat: Bab 4 (hidrologi: Sungai.gpkg) dan Bab 5 (CH_tahunan.tif) sudah dijalankan; keluarannya ada di <paket-m2>/hasil."""
import os
import sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import _inisialisasi  # noqa: F401
import processing
from osgeo import gdal

gdal.UseExceptions()
paket = sys.argv[1]
out = os.path.join(paket, "hasil")
P = lambda n: os.path.join(out, n)
EXT = "312000,314000,9996000,9998000 [EPSG:32749]"
dem = os.path.join(paket, "DEM_10m.tif")
JARAK_SUNGAI = 30                                            # m: ambang kendala dari sungai (dipakai di rumus QGIS dan pemeriksaan NumPy)


def baca(path):
    ds = gdal.Open(path)
    return ds.GetRasterBand(1).ReadAsArray().astype("float64")


def reklas(inp, tabel, hasil, batas=1):
    flat = [v for baris in tabel for v in baris]
    processing.run("native:reclassifybytable", {"INPUT_RASTER": inp, "RASTER_BAND": 1, "TABLE": flat, "NO_DATA": 0, "RANGE_BOUNDARIES": batas,
                                                "NODATA_FOR_MISSING": True, "DATA_TYPE": 0, "OUTPUT": hasil})


# ---------------------------------------------------------------- 1. bobot dengan AHP (perbandingan berpasangan)
nama = ["lereng", "tanah", "curah hujan", "jarak jalan"]
A = np.array([[1, 2, 3, 4], [1 / 2, 1, 2, 3], [1 / 3, 1 / 2, 1, 2], [1 / 4, 1 / 3, 1 / 2, 1]])
vals, vecs = np.linalg.eig(A)
k = int(np.argmax(vals.real))
w = np.abs(vecs[:, k].real)
w = w / w.sum()
CI = (vals[k].real - len(A)) / (len(A) - 1)
CR = CI / 0.90                                               # RI untuk n = 4 adalah 0,90 (tabel Saaty)
print("bobot AHP:", dict(zip(nama, np.round(w, 3))), "| lambda_maks = %.3f | CI = %.4f | CR = %.4f" % (vals[k].real, CI, CR))

# ---------------------------------------------------------------- 2. kriteria
processing.run("gdal:slope", {"INPUT": dem, "BAND": 1, "SCALE": 1, "AS_PERCENT": True, "COMPUTE_EDGES": True, "ZEVENBERGEN": False,
                              "OUTPUT": P("Lereng_pct.tif")})
reklas(P("Lereng_pct.tif"), [[0, 8, 4], [8, 15, 3], [15, 25, 2], [25, 40, 1]], P("Skor_lereng.tif"))
reklas(P("CH_tahunan.tif"), [[0, 2400, 1], [2400, 2600, 2], [2600, 2800, 3], [2800, 99999, 4]], P("Skor_hujan.tif"))

tanah = processing.run("native:fieldcalculator", {
    "INPUT": os.path.join(paket, "Tanah.gpkg|layername=Tanah"), "FIELD_NAME": "skor", "FIELD_TYPE": 1, "FIELD_LENGTH": 3, "FIELD_PRECISION": 0,
    "FORMULA": "CASE WHEN \"jenis_tanah\"='Aluvial' THEN 4 WHEN \"jenis_tanah\"='Latosol' THEN 3 WHEN \"jenis_tanah\"='Podsolik' THEN 2 ELSE 1 END",
    "OUTPUT": "memory:"})["OUTPUT"]
processing.run("gdal:rasterize", {"INPUT": tanah, "FIELD": "skor", "BURN": 0, "USE_Z": False, "UNITS": 0, "WIDTH": 200, "HEIGHT": 200, "EXTENT": EXT,
                                  "NODATA": 0, "DATA_TYPE": 0, "INIT": 0, "INVERT": False, "OUTPUT": P("Skor_tanah.tif")})


def jarak(vektor, tmp, hasil):
    processing.run("gdal:rasterize", {"INPUT": vektor, "FIELD": "", "BURN": 1, "USE_Z": False, "UNITS": 0, "WIDTH": 200, "HEIGHT": 200,
                                      "EXTENT": EXT, "NODATA": 0, "DATA_TYPE": 0, "INIT": 0, "INVERT": False, "OUTPUT": tmp})
    processing.run("gdal:proximity", {"INPUT": tmp, "BAND": 1, "VALUES": "1", "UNITS": 0, "MAX_DISTANCE": 0, "REPLACE": 0, "NODATA": -9999,
                                      "DATA_TYPE": 5, "OUTPUT": hasil})


jarak(os.path.join(paket, "Jalan.gpkg|layername=Jalan|subset=status<>'ditutup'"), P("_jalan_sel.tif"), P("Jarak_jalan.tif"))
jarak(P("Sungai.gpkg"), P("_sungai_sel.tif"), P("Jarak_sungai.tif"))
reklas(P("Jarak_jalan.tif"), [[0, 100, 4], [100, 250, 3], [250, 500, 2], [500, 1e9, 1]], P("Skor_jalan.tif"))

# ---------------------------------------------------------------- 3. kendala (larangan): lereng >= 40% atau dalam 30 m dari sungai
processing.run("gdal:rastercalculator", {"INPUT_A": P("Skor_lereng.tif"), "BAND_A": 1, "INPUT_B": P("Jarak_sungai.tif"), "BAND_B": 1,
                                         "FORMULA": "logical_and(A>0, B>%d)" % JARAK_SUNGAI, "NO_DATA": 255, "RTYPE": 0, "EXTENT_OPT": 0, "EXTRA": "--hideNoData", "OUTPUT": P("Kendala_boleh.tif")})

# ---------------------------------------------------------------- 4. tumpang tindih berbobot
rumus = "A*%.6f + B*%.6f + C*%.6f + D*%.6f" % tuple(w)
processing.run("gdal:rastercalculator", {"INPUT_A": P("Skor_lereng.tif"), "BAND_A": 1, "INPUT_B": P("Skor_tanah.tif"), "BAND_B": 1,
                                         "INPUT_C": P("Skor_hujan.tif"), "BAND_C": 1, "INPUT_D": P("Skor_jalan.tif"), "BAND_D": 1,
                                         "FORMULA": rumus, "NO_DATA": -9999, "RTYPE": 5, "EXTENT_OPT": 0, "EXTRA": "--hideNoData", "OUTPUT": P("Kesesuaian_skor_mentah.tif")})
processing.run("gdal:rastercalculator", {"INPUT_A": P("Kesesuaian_skor_mentah.tif"), "BAND_A": 1, "INPUT_B": P("Kendala_boleh.tif"), "BAND_B": 1,
                                         "FORMULA": "(B==1)*(1*(A>=3.2) + 2*logical_and(A>=2.6, A<3.2) + 3*logical_and(A>=2.0, A<2.6) + 4*(A<2.0))",
                                         "NO_DATA": 255, "RTYPE": 0, "EXTENT_OPT": 0, "EXTRA": "--hideNoData", "OUTPUT": P("Kesesuaian_kelas.tif")})
kelas = baca(P("Kesesuaian_kelas.tif"))
label = {0: "Terlarang (kendala)", 1: "S1 Sangat sesuai", 2: "S2 Cukup sesuai", 3: "S3 Sesuai marginal", 4: "N Tidak sesuai"}
print("\nLuas per kelas (ha), 1 sel = 100 m2 = 0,01 ha:")
for kd in range(5):
    print("  %-22s %6d sel  %7.2f ha  %5.1f%%" % (label[kd], (kelas == kd).sum(), (kelas == kd).sum() * 0.01, 100.0 * (kelas == kd).mean()))
print("  jumlah: %.2f ha (sel NoData: %d)" % (sum((kelas == kd).sum() for kd in range(5)) * 0.01, (kelas == 255).sum()))

# ---------------------------------------------------------------- 5. pemeriksaan silang dengan NumPy murni
s = [baca(P("Skor_lereng.tif")), baca(P("Skor_tanah.tif")), baca(P("Skor_hujan.tif")), baca(P("Skor_jalan.tif"))]
skor = sum(wi * si for wi, si in zip(w, s))
boleh = (s[0] > 0) & (baca(P("Jarak_sungai.tif")) > JARAK_SUNGAI)
np_kelas = np.where(boleh, np.select([skor >= 3.2, skor >= 2.6, skor >= 2.0], [1, 2, 3], 4), 0)
print("\nPemeriksaan NumPy: sel kelas sama = %.4f%%" % (100.0 * (np_kelas == kelas).mean()))

# ---------------------------------------------------------------- 6. kepekaan: bobot sama rata
skor2 = sum(0.25 * si for si in s)
np2 = np.where(boleh, np.select([skor2 >= 3.2, skor2 >= 2.6, skor2 >= 2.0], [1, 2, 3], 4), 0)
print("Kepekaan: luas S1 bobot AHP = %.2f ha; bobot sama rata = %.2f ha; sel yang berubah kelas = %.1f%%" %
      ((np_kelas == 1).sum() * 0.01, (np2 == 1).sum() * 0.01, 100.0 * (np2 != np_kelas).mean()))

# ---------------------------------------------------------------- 7. per petak (histogram zona) dan poligon untuk basis data
zh = processing.run("native:zonalhistogram", {"INPUT_RASTER": P("Kesesuaian_kelas.tif"), "RASTER_BAND": 1,
                                              "INPUT_VECTOR": os.path.join(paket, "Petak.gpkg|layername=Petak"), "COLUMN_PREFIX": "k_",
                                              "OUTPUT": "memory:"})["OUTPUT"]
baris = []
for f in zh.getFeatures():
    baris.append((f["kode"], f["jenis_tegakan"], *(f["k_%d" % kd] * 0.01 for kd in (1, 2, 3, 4, 0))))
baris.sort(key=lambda r: -r[2])
print("\nLima petak dengan S1 terluas (ha): kode, jenis, S1, S2, S3, N, Terlarang")
for r in baris[:5]:
    print("  %s %-12s %6.2f %6.2f %6.2f %6.2f %6.2f" % r)
print("  total S1 dari 25 petak = %.2f ha" % sum(r[2] for r in baris))
mentah = processing.run("gdal:polygonize", {"INPUT": P("Kesesuaian_kelas.tif"), "BAND": 1, "FIELD": "kelas", "EIGHT_CONNECTEDNESS": False, "EXTRA": "",
                                            "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
if os.path.exists(P("Kesesuaian_kelas.gpkg")):
    os.remove(P("Kesesuaian_kelas.gpkg"))
processing.run("native:dissolve", {"INPUT": mentah, "FIELD": ["kelas"], "SEPARATE_DISJOINT": False, "OUTPUT": P("Kesesuaian_kelas.gpkg")})
for tmp in ("_jalan_sel.tif", "_sungai_sel.tif"):
    try:
        os.remove(P(tmp))
    except OSError:
        pass
print("\npoligon kelas ditulis:", os.path.exists(P("Kesesuaian_kelas.gpkg")))

Skrip ini 124 baris dan dipecah menjadi tujuh blok.

Blok 1, baris 1-31: persiapan dan fungsi reklasifikasi.

Blok 2, baris 34-43: bobot AHP.

Blok 3, baris 45-68: menyiapkan empat kriteria.

Blok 4 dan 5, baris 70-87: kendala, tumpang tindih, dan kelas.

Blok 6, baris 89-100: pemeriksaan silang dan kepekaan.

Blok 7, baris 102-124: per petak dan poligon.

Hasil. Layar menampilkan (peringatan GDAL tentang lebar kolom teks boleh diabaikan):

bobot AHP: {'lereng': np.float64(0.467), 'tanah': np.float64(0.277), 'curah hujan': np.float64(0.16), 'jarak jalan': np.float64(0.095)} | lambda_maks = 4.031 | CI = 0.0103 | CR = 0.0115

Luas per kelas (ha), 1 sel = 100 m2 = 0,01 ha:
  Terlarang (kendala)      5321 sel    53.21 ha   13.3%
  S1 Sangat sesuai         6733 sel    67.33 ha   16.8%
  S2 Cukup sesuai         12345 sel   123.45 ha   30.9%
  S3 Sesuai marginal      11627 sel   116.27 ha   29.1%
  N Tidak sesuai           3974 sel    39.74 ha    9.9%
  jumlah: 400.00 ha (sel NoData: 0)

Pemeriksaan NumPy: sel kelas sama = 100.0000%
Kepekaan: luas S1 bobot AHP = 67.33 ha; bobot sama rata = 81.76 ha; sel yang berubah kelas = 31.6%

Lima petak dengan S1 terluas (ha): kode, jenis, S1, S2, S3, N, Terlarang
  P-05 Jati          11.95   5.36   0.36   0.01   1.47
  P-03 Mahoni         9.02   5.27   1.63   0.18   3.91
  P-08 Jati           4.46   5.46   2.96   0.00   2.95
  P-17 Jati           4.26   4.78   3.55   1.47   5.84
  P-06 Akasia         4.24   4.95   1.83   0.01   2.53
  total S1 dari 25 petak = 67.33 ha

poligon kelas ditulis: True

Cara membaca: lereng paling berpengaruh (bobot 0,467). Penilaian konsisten (CR 0,0115, jauh di bawah 0,1). Kelas S1 seluas 67,33 ha (16,8 persen), dan 13,3 persen wilayah terlarang oleh kendala. Pemeriksaan NumPy cocok 100 persen. Bila bobot disama-ratakan, S1 melebar ke 81,76 ha: pilihan bobot sangat menentukan. Petak P-05 adalah yang paling banyak S1.

Bagian yang boleh diubah.

• Baris 36, matriks AHP. Ubah angka perbandingan. Jaga agar sisi bawah tetap kebalikan sisi atas.
• Baris 48-49 dan 68, tabel skor. Ubah batas rentang dan skor.
• Baris 20, jarak aman sungai. Ubah angka 30 (meter) di JARAK_SUNGAI. Cukup di satu tempat: rumus kendala (baris 72) dan pemeriksaan NumPy (baris 92) ikut berubah.
• Baris 80, batas kelas. Ubah 3.2, 2.6, dan 2.0.
• Baris 86, luas sel. 0.01 hektare benar untuk sel 10 meter.

Jebakan.

• Bergantung pada 25.7 dan 25.8. Bila CH_tahunan.tif belum ada, hasilnya QgsProcessingException: Unable to execute algorithm dengan keterangan Could not load source layer for INPUT_RASTER: ...CH_tahunan.tif not found.
• Satu aturan, satu nama. Syarat sungai 30 meter disimpan di JARAK_SUNGAI (baris 20), lalu dipakai di rumus QGIS (baris 72) dan pemeriksaan NumPy (baris 92). Karena itu keduanya selalu sama dan pemeriksaan tetap 100 persen. Jangan menulis angka 30 langsung di dua baris itu. Bila ada yang diganti di satu tempat saja, kedua hitungan berbeda dan angka pemeriksaan turun di bawah 100 persen. Itulah alasan angka dijadikan satu nama.
• Skor 0 punya dua arti. NO_DATA: 0 memberi skor 0 pada nilai di luar rentang, dan baris 72 memakai A>0 sebagai tanda "lereng terlalu curam". Bila Anda menambah skor sah 0, kendalanya bisa tercampur. [kemungkinan]
• Jenis tanah tak terdaftar mendapat skor 1. ELSE 1 di rumus tanah (baris 53) berlaku untuk nama tanah apa pun yang tidak disebut.
• CR memakai angka acuan tetap. 0.90 hanya benar untuk matriks 4 x 4. Matriks lain butuh angka acuan lain.
• Luas sel tertulis tetap. 0.01 ha di beberapa tempat. Salah bila sel bukan 10 meter.
• Bobot dicetak dengan tulisan aneh. Tulisan np.float64(...) itu bentuk cetak NumPy versi baru, bukan galat.
• Kriteria fiktif. Jangan memakai ambang ini untuk keputusan nyata.

Latihan variasi. Perlebar jarak aman dari sungai menjadi 50 meter. Contoh jawaban: ganti JARAK_SUNGAI = 30 menjadi JARAK_SUNGAI = 50 di baris 20, cukup di satu tempat. Layar menampilkan antara lain Terlarang (kendala) 8249 sel 82.49 ha 20.6% dan S1 Sangat sesuai 6162 sel 61.62 ha 15.4%, dan Pemeriksaan NumPy: sel kelas sama = 100.0000%.

Skrip 25.10: Cadangan dan pemulihan dengan Bash #

Berkas m2_07d_cadangan.sh. Ini bukan Python, melainkan skrip Bash: bahasa perintah di terminal Linux, macOS, dan Git Bash di Windows. Isinya empat perintah PostgreSQL yang dibungkus supaya bisa diulang. Cadangan yang belum pernah diuji pemulihannya bukan cadangan, jadi skrip ini juga mengujinya.

Tujuan. Menyalin skema kph ke berkas cadangan (pg_dump), memulihkannya ke basis data baru (pg_restore), lalu membandingkan jumlah baris tabel asal dan hasil pulih.

Bahan. Alat PostgreSQL (pg_dump, pg_restore, psql) yang ada di PATH (daftar folder pencarian perintah) atau lewat variabel PG_BIN. Basis data kph_contoh dengan skema kph dan tabel petak, kejadian, jalan, sungai, kesesuaian, dan lembar_keputusan (tiga terakhir berasal dari hasil Skrip 25.7, 25.9, dan m2_07a_lembar_keputusan.sql). Hak membuat basis data. Kata sandi lewat pgpass atau PGPASSWORD, tidak ditulis di skrip. Perintah:

bash m2_07d_cadangan.sh D:/Latihan/cadangan

Skrip utuh:

#!/usr/bin/env bash
# M2 Bab 7: cadangan (pg_dump) dan pemulihan (pg_restore) skema kph, lalu pemeriksaan jumlah baris.
# Penulis: Badar Mubarok Yogaswara. Koneksi lewat PGHOST, PGPORT, PGUSER (kata sandi: berkas pgpass atau PGPASSWORD, JANGAN ditulis di skrip).
# Pemakaian: bash m2_07d_cadangan.sh <folder tujuan berkas cadangan> ; PG_BIN bisa diisi folder bin PostgreSQL bila tidak ada di PATH.
set -euo pipefail
B="${PG_BIN:-}"; [ -n "$B" ] && B="$B/"
TUJUAN="${1:-.}"
BERKAS="$TUJUAN/kph_contoh.dump"
"${B}pg_dump" -Fc -n kph -f "$BERKAS" kph_contoh                                    # -Fc = format kustom (terkompres), -n = hanya skema kph
echo "ukuran cadangan: $(wc -c < "$BERKAS") byte"
"${B}psql" -d postgres -c "DROP DATABASE IF EXISTS kph_pulih" -c "CREATE DATABASE kph_pulih"
"${B}psql" -d kph_pulih -c "CREATE EXTENSION postgis" -c "CREATE EXTENSION pgrouting"   # ekstensi dibuat dulu di basis data tujuan
"${B}pg_restore" -d kph_pulih "$BERKAS"
for t in petak kejadian jalan sungai kesesuaian lembar_keputusan; do
  a=$("${B}psql" -d kph_contoh -Atc "SELECT count(*) FROM kph.$t"); b=$("${B}psql" -d kph_pulih -Atc "SELECT count(*) FROM kph.$t")
  echo "$t asal=$a pulih=$b $([ "$a" = "$b" ] && echo SAMA || echo BEDA)"
done

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas.

Hasil. Pada uji, layar menampilkan:

ukuran cadangan: 114417 byte
DROP DATABASE
CREATE DATABASE
CREATE EXTENSION
CREATE EXTENSION
petak asal=25 pulih=25 SAMA
kejadian asal=150 pulih=150 SAMA
jalan asal=37 pulih=37 SAMA
sungai asal=30 pulih=30 SAMA
kesesuaian asal=5 pulih=5 SAMA
lembar_keputusan asal=25 pulih=25 SAMA

Angka byte berbeda tiap basis data, jadi angka Anda mungkin lain. Bila basis data kph_pulih belum ada, psql juga mencetak pemberitahuan NOTICE bahwa basis data itu tidak ada dan dilewati. Itu normal. Enam tabel sama persis, jadi cadangan terbukti bisa dipulihkan.

Bagian yang boleh diubah.

• Baris 9 dan 15, nama basis data sumber. Ganti kph_contoh.
• Baris 11-13 dan 15, nama basis data pulih. Ganti kph_pulih. Ubah di semua tempat.
• Baris 14, daftar tabel. Tambah atau kurangi sesuai isi skema Anda.
• Argumen folder tujuan. Ganti folder di perintah menjalankan.

Jebakan.

• Alat PostgreSQL tidak ada di PATH. Bila pg_dump tidak ditemukan, skrip berhenti di baris 9 dengan pg_dump: command not found. Isi PG_BIN dengan folder bin PostgreSQL, atau tambahkan ke PATH. Karena set -e, tidak ada perintah penghapus yang sempat jalan.
• Basis data pulih terhapus. Baris 11 menghapus basis data dengan nama itu. Ini aksi destruktif.
• Tabel belum ada. Tabel sungai, kesesuaian, dan lembar_keputusan berasal dari langkah lain. Bila belum ada di basis data Anda, psql berhenti dengan galat. [kemungkinan]
• Peran dan hak tidak ikut. Cadangan hanya skema kph. Peran pengguna dan hak akses perlu disiapkan terpisah di server tujuan. [kemungkinan]
• Kata sandi. Skrip tidak membawanya. Siapkan pgpass.
• Bukan Python. Skrip ini tidak jalan dengan python. Pakai bash.

Latihan variasi. Tambahkan tabel stasiun_hujan ke daftar pemeriksaan. Contoh jawaban: ubah baris 14 menjadi:

for t in petak kejadian jalan sungai kesesuaian lembar_keputusan stasiun_hujan; do

Tabel kph.stasiun_hujan berisi 30 baris pada data ini, jadi baris tambahan yang tercetak adalah stasiun_hujan asal=30 pulih=30 SAMA (diuji pada salinan skrip).

Sembilan berkas SQL yang hanya disebut #

Folder skrip M2 juga berisi sembilan berkas SQL. Berkas-berkas ini bukan Python, jadi tidak dibahas baris demi baris di sini. Anda menjalankannya dengan psql -f nama_berkas.sql. Dua di antaranya (m2_03b_pgrouting.sql dan m2_07a_lembar_keputusan.sql) dijalankan pada uji bab ini. Tujuh lainnya hanya dibaca judulnya. [CEK]

Cek paham #

1. Mengapa skrip PostGIS di bab ini tidak memuat alamat server atau kata sandi, dan dari mana skrip tahu ke mana harus tersambung?
2. Apa arti HASIL SAMA: True pada Skrip 25.3, dan mengapa membandingkan dua sistem berguna?
3. Mengapa Skrip 25.9 harus dijalankan setelah Skrip 25.7 dan 25.8?

Jawaban.

1. Karena kata sandi dan alamat tidak boleh bocor lewat berkas skrip. Skrip membaca PGHOST, PGPORT, PGUSER, dan PGDATABASE dari variabel lingkungan, dan sandi dari berkas pgpass atau PGPASSWORD.
2. Tiga jawaban (total luas, tiga petak teratas, kejadian dekat jalan) dari GeoPackage persis sama dengan dari PostGIS. Bila dua jalan berbeda memberi angka sama, hitungan kemungkinan besar benar.
3. Karena skrip itu membaca Sungai.gpkg (hasil 25.7) dan CH_tahunan.tif (hasil 25.8). Tanpa keduanya, alat berhenti dengan galat "tidak ditemukan".

Kesalahan umum #

• Lupa argumen folder. Hampir semua skrip butuh folder paket. Perbaikan: ketik folder setelah nama skrip, atau ganti sys.argv[1] dengan jalur manual bila memakai Konsol QGIS.
• Mengira Skrip 25.1 tidak menghapus apa-apa. Skrip ini aman dijalankan ulang karena OVERWRITE=YES mengganti tabel lama. Akibatnya, perubahan yang Anda buat langsung di tabel hilang. Perbaikan: simpan perubahan di skrip atau berkas sumber, bukan di tabel hasil impor.
• Mengganti angka di satu tempat, padahal ada di beberapa tempat. Contohnya jarak 100 meter (Skrip 25.3) dan nomor simpul (Skrip 25.6). Perbaikan: cari semua tulisan angka itu di skrip. Skrip 25.9 sudah menyimpan jarak sungai dalam satu nama, JARAK_SUNGAI, sehingga cukup diganti sekali.
• Menganggap "dilewati" tidak penting. Pesan Pembanding PostGIS dilewati bisa berarti server mati atau salah ketik. Perbaikan: baca pesan sesudah titik dua.
• Memberi nama berkas sama dengan modul bawaan Python. Berkas bernama enum.py di folder kerja membuat Python QGIS gagal dengan ImportError: cannot import name 'NULL' from partially initialized module (terjadi pada uji bab ini). Perbaikan: beri nama lain.

Ringkasan dan latihan #

• PostGIS adalah gudang data bersama. Skrip bicara kepadanya lewat GDAL (Skrip 25.1), QGIS (25.2), dan psycopg2 (25.3, 25.6).
• Koneksi dibaca dari variabel lingkungan PG*. Kata sandi tidak ditulis di skrip.
• Menyilang dua cara (GeoPackage dan PostGIS, QGIS dan pgRouting, alat QGIS dan NumPy) adalah cara terbaik menguji hitungan sendiri.
• Lambda adalah fungsi mini satu baris. Boleh ditulis ulang sebagai def bila lebih mudah dibaca.
• Analisis lahan (hidrologi, interpolasi, kesesuaian) berantai: keluaran satu skrip menjadi masukan skrip berikutnya.
• Jangan lupa menguji cadangan dengan memulihkannya (Skrip 25.10).

Latihan. Pilih satu skrip di bab ini dan temukan satu angka yang tertulis tetap (misalnya 100 di Skrip 25.3 atau 500). Tulis di bagian mana saja angka itu muncul, dan apa yang terjadi bila Anda mengubahnya hanya di satu tempat.

Tabel perbandingan: alat apa untuk pekerjaan apa #