BAB 31: Dua Gaya, Satu Pekerjaan: ArcPy di ArcGIS Pro dan ArcMap

Studi kasus: "Kantor cabang memakai ArcGIS" #

Kepala Seksi membawa kabar: kantor cabang tidak memakai QGIS. Sebagian komputer di sana memakai ArcGIS Pro, sebagian lagi masih ArcMap 10.8. Pekerjaannya sama persis: membaca batas KPH, menghitung luas, dan memasukkan layer ke peta. Yang berbeda hanya cara menuliskannya. Anda perlu tahu dua gaya itu supaya skrip bisa berpindah dari satu komputer ke komputer lain.

Konsep: ArcPy dalam tiga kalimat #

ArcPy adalah perkakas Python bawaan ArcGIS. Bayangkan ArcGIS sebagai kantor dengan banyak loket, dan ArcPy sebagai buku telepon tiap loket: satu panggilan, satu pekerjaan selesai. ArcGIS Pro dan ArcMap punya buku telepon yang hampir sama, tetapi sebagian nomornya ditulis berbeda.

Perbedaan terbesarnya ada di versi Python. ArcGIS Pro memakai Python 3. ArcMap 10.8 memakai Python 2.7, versi lama yang cara menulis beberapa hal berbeda. Karena itu setiap skrip di bab ini hadir dalam dua versi.

Lima istilah yang perlu Anda kenal sekarang:

• arcpy.env.workspace: folder kerja bawaan. Nama data yang Anda tulis dicari di folder ini.
• Describe: perintah untuk melihat sifat sebuah data, misalnya jenis bentuk dan sistem koordinatnya.
• Kursor (cursor): penunjuk yang membaca tabel baris demi baris, seperti jari yang menelusuri daftar.
• Proyek (di Pro) dan dokumen peta (di ArcMap): berkas yang menyimpan peta Anda.
• Layer: satu lapis data di peta. Istilah ini sama dengan di QGIS.

Skrip 31.1: Membaca data, atribut, dan proyek #

Tujuan. Melihat daftar data di sebuah folder, memeriksa sifat batas KPH, membaca nama dan luas tiap KPH, membuka sebuah raster, lalu memasukkan layer ke peta.

Bahan. Paket data paket-i4: folder shp (berisi Batas_KPH.shp) dan berkas Citra_KPH.tif. Bagian 4 skrip hanya berjalan dari jendela Python di dalam ArcGIS Pro, karena memakai proyek yang sedang terbuka. [CEK]

Versi ArcGIS Pro (Python 3) #

# [SKRIP 1.1-Pro: ArcPy dasar: membaca data, atribut, dan proyek]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi ArcGIS Pro Anda.
import arcpy

DATA = r"D:\Latihan\paket-i4"
arcpy.env.workspace = DATA + r"\shp"          # folder Shapefile; daftar data dibaca dari sini

# 1. Daftar data dan sifatnya
print(arcpy.ListFeatureClasses())              # contoh: ['Batas_KPH.shp', 'Hotspot_Nasional.shp', ...]
batas = "Batas_KPH.shp"
d = arcpy.Describe(batas)
print(d.shapeType, d.spatialReference.name, d.spatialReference.factoryCode)   # Polygon, ..., 32749
print("Jumlah fitur:", int(arcpy.management.GetCount(batas)[0]))
print("Field:", [f.name for f in arcpy.ListFields(batas)])

# 2. Membaca tiap baris dengan kursor (padanan getFeatures() di QGIS)
with arcpy.da.SearchCursor(batas, ["NAMA_KPH", "SHAPE@AREA"]) as kursor:
    for nama, luas in kursor:
        print(nama, round(luas / 10000, 2), "ha")          # luas dalam satuan peta (m2) dibagi 10.000

# 3. Raster
r = arcpy.Raster(DATA + r"\Citra_KPH.tif")
print("Raster:", r.width, "x", r.height, "piksel,", r.bandCount, "band; piksel", r.meanCellWidth, "m")

# 4. Proyek: tambah data lalu ambil layer dengan nama
aprx = arcpy.mp.ArcGISProject("CURRENT")      # proyek yang sedang terbuka
peta = aprx.listMaps()[0]
peta.addDataFromPath(DATA + r"\shp\Batas_KPH.shp")
layer = peta.listLayers("Batas_KPH")[0]       # listLayers() juga mengembalikan DAFTAR
print("Layer ditemukan:", layer.name)
aprx.save()

Alur skrip ini mirip skrip QGIS yang sudah Anda kenal: siapkan lokasi data, baca sifatnya, ulangi tiap baris, lalu masukkan ke peta.

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

Hasil yang diharapkan. Berikut perkiraan keluarannya. Angka-angka ini dihitung dari data dengan alat lain, bukan dijalankan di ArcGIS, jadi tulisan persisnya bisa sedikit berbeda. [CEK]

Jumlah fitur: 3
KPH Alpha 4.0 ha
KPH Beta 4.0 ha
KPH Gamma 8.0 ha
Raster: 200 x 200 piksel, 4 band; piksel 2.0 m
Layer ditemukan: Batas_KPH

Baris sistem koordinat memuat angka 32749, yaitu kode UTM zona 49 Selatan.

Versi ArcMap 10.8 (Python 2.7) #

# -*- coding: utf-8 -*-
# [SKRIP 1.1-ArcMap: ArcPy dasar: membaca data, atribut, dan dokumen peta (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi ArcMap Anda.
import arcpy

DATA = r"D:\Latihan\paket-i4"
arcpy.env.workspace = DATA + r"\shp"

# 1. Daftar data dan sifatnya
print(arcpy.ListFeatureClasses())
batas = "Batas_KPH.shp"
d = arcpy.Describe(batas)
print("%s %s %s" % (d.shapeType, d.spatialReference.name, d.spatialReference.factoryCode))
print("Jumlah fitur: %d" % int(arcpy.GetCount_management(batas).getOutput(0)))
print("Field: %s" % [f.name for f in arcpy.ListFields(batas)])

# 2. Membaca tiap baris dengan kursor
with arcpy.da.SearchCursor(batas, ["NAMA_KPH", "SHAPE@AREA"]) as kursor:
    for nama, luas in kursor:
        print("%s %.2f ha" % (nama, luas / 10000.0))        # 10000.0: hasil bagi desimal di Python 2

# 3. Raster
r = arcpy.Raster(DATA + r"\Citra_KPH.tif")
print("Raster: %d x %d piksel, %d band, piksel %s m" % (r.width, r.height, r.bandCount, r.meanCellWidth))

# 4. Dokumen peta: tambah layer lalu ambil dengan nama
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
arcpy.mapping.AddLayer(df, arcpy.mapping.Layer(DATA + r"\shp\Batas_KPH.shp"))
layer = arcpy.mapping.ListLayers(mxd, "Batas_KPH", df)[0]    # juga mengembalikan DAFTAR
print("Layer ditemukan: %s" % layer.name)
mxd.save()

Bagian yang boleh diubah. Nomor baris mengikuti versi ArcGIS Pro.

• Baris 6, lokasi data. Ganti D:\Latihan\paket-i4 dengan lokasi folder paket-i4 di komputer Anda. Pertahankan huruf r di depan kutip.
• Baris 11, nama data. Ganti Batas_KPH.shp dengan data lain di folder shp, misalnya Sungai.shp.
• Baris 18, nama kolom. Ganti NAMA_KPH dengan kolom lain. Daftar kolom ada di hasil baris 15.
• Baris 23, raster. Ganti Citra_KPH.tif dengan raster lain di folder yang sama.
• Baris 30, nama layer. Ganti Batas_KPH sesuai nama layer di peta Anda.

Jebakan.

• Satuan luas. SHAPE@AREA memakai satuan peta. Pada data UTM dalam meter, hasilnya meter persegi, dan pembagian 10.000 mengubahnya menjadi hektare. Pada data berkoordinat derajat, angkanya bukan hektare.
• CURRENT di luar ArcGIS. Bagian 4 butuh proyek (atau dokumen peta) yang sedang terbuka. Dari terminal biasa, bagian ini gagal. Jalankan dari jendela Python di dalam aplikasinya. [CEK]
• print dengan banyak isian di ArcMap. Gaya print(nama, luas) milik Python 3 bisa mencetak tanda kurung dan koma di Python 2.7. Pakai satu teks dengan %.
• Garis miring terbalik. Lokasi data tanpa huruf r di depan kutip bisa dibaca salah, karena \ dipakai Python untuk tanda khusus.

Latihan variasi. Ubah bagian 2 supaya hanya KPH yang lebih luas dari 5 hektare yang tercetak. Contoh jawaban untuk ArcGIS Pro:

with arcpy.da.SearchCursor(batas, ["NAMA_KPH", "SHAPE@AREA"]) as kursor:
    for nama, luas in kursor:
        if luas / 10000 > 5:
            print(nama, round(luas / 10000, 2), "ha")

Dari data contoh, hanya satu baris yang tercetak: KPH Gamma 8.0 ha. [CEK]

Skrip 31.2: Rantai alat, dari sempadan sungai sampai NDVI per KPH #

Berkas b2_arcpy_pro.py dan b2_arcpy_arcmap.py. Di dalam berkasnya tertulis "Skrip 2.1" mengikuti penomoran Buku I4. Padanan PyQGIS-nya ada di Bab 21.

Tujuan. Menjalankan beberapa alat ArcPy berurutan: membuat sempadan sungai 30 m, memotongnya dengan batas KPH, menghitung titik panas per KPH, lalu menghitung NDVI dan rata-ratanya untuk tiap KPH.

Bahan. Paket data paket-i4: folder shp (Sungai.shp, Batas_KPH.shp, Hotspot_Nasional.shp) dan Citra_KPH.tif (empat band: merah, hijau, biru, inframerah dekat). Siapkan juga folder D:\Latihan\hasil_i4: skrip tidak membuatnya. Ekstensi Spatial Analyst harus tersedia.

Enam istilah baru:

• Alat geoprocessing: satu pekerjaan siap pakai di ArcGIS, misalnya Buffer. Di skrip, alat ditulis sebagai perintah ArcPy.
• Sempadan (buffer): daerah selebar jarak tertentu di kiri-kanan sebuah garis.
• Ruang sementara (memory di Pro): tempat menyimpan hasil antara di memori komputer, tanpa menulis berkas.
• Kamus (dictionary): daftar pasangan "nama dan nilai". Dibahas di Bab 5.
• Statistik zonal: menghitung angka ringkasan, misalnya rata-rata, untuk nilai raster di dalam tiap poligon.
• Join_Count: kolom jumlah pasangan yang ditambahkan alat Spatial Join.

Skrip ini bekerja dalam enam tahap.

Versi ArcGIS Pro (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")

Penjelasan baris demi baris. Nomor baris mengikuti berkas di atas. Baris import, ekstensi, dan folder kerja sudah Anda kenal di Skrip 31.1 dan Bab 30.

Dokumentasi Esri mencatat bahwa ZonalStatisticsAsTable memerlukan Spatial Analyst atau Image Analyst, dengan urutan isian in_zone_data, zone_field, in_value_raster, out_table, ignore_nodata, statistics_type. Urutan itu cocok dengan baris 33. Dokumentasi Buffer memberi contoh jarak berformat "30 Meters", dan nilai "ALL" pada dissolve_option menggabungkan semua sempadan menjadi satu fitur. Dokumentasi SpatialJoin menyebut kolom Join_Count ditambahkan otomatis. [CEK]

Versi ArcMap 10.8 (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")

Hasil. Layar menampilkan enam baris, dan NDVI_KPH.tif muncul di folder hasil. Berikut perkiraannya. Luas dan sempadan dihitung dengan GDAL dan OGR dari data contoh, dan sejalan dengan README.txt paket data. Urutan baris mengikuti urutan tabel. [CEK]

KPH Alpha: luas 4.00 ha, sempadan 1.18 ha, titik panas 2
KPH Beta: luas 4.00 ha, sempadan 1.39 ha, titik panas 5
KPH Gamma: luas 8.00 ha, sempadan 1.50 ha, titik panas 5
KPH Alpha: NDVI rata-rata 0.604
KPH Beta: NDVI rata-rata 0.431
KPH Gamma: NDVI rata-rata 0.329

Skrip PyQGIS padanannya di Bab 21 menghasilkan angka yang sama. Logika Python skrip ini sudah diuji dengan modul ArcPy tiruan buatan sendiri yang membaca paket data. Perilaku alat ArcGIS yang sebenarnya belum teruji. [CEK]

Bagian yang boleh diubah. Nomor baris mengikuti versi Pro.

• Baris 7-8, lokasi data dan hasil. Ganti sesuai komputer Anda. Pertahankan huruf r.
• Baris 14, jarak. Ganti "30 Meters" untuk jarak lain. Satuan ikut ditulis dalam teksnya.
• Baris 29-30, nomor band. Ganti Band_4 dan Band_1 bila urutan band citra Anda berbeda.
• Baris 33, jenis statistik. Ganti "MEAN" untuk statistik lain yang didukung alat. Daftarnya ada di dokumentasi alat. [CEK]

Jebakan.

• Folder hasil tidak dibuat. Bila D:\Latihan\hasil_i4 belum ada, baris 32 kemungkinan besar gagal. Skrip 31.3 menunjukkan cara membuat folder lebih dulu. [CEK]
• Ruang sementara menumpuk. Di Pro, hanya memory\sempadan yang dihapus. Objek irisan, kph_titik, dan ndvi_zonal bertahan sampai sesi berakhir. [CEK]
• Satuan luas. Pembagian 10.000 benar hanya bila data berukuran meter (UTM). Pada data berkoordinat derajat, angkanya bukan hektare.
• Band tertukar. Bila band 1 dan 4 tertukar, NDVI berubah tanda tanpa peringatan.
• Nama band dalam jalur. Penulisan Citra_KPH.tif\Band_4 mengandaikan ArcGIS mengenali band lewat jalur seperti itu. [CEK]
• Angka antaralat sedikit berbeda. Sempadan bergantung pada jumlah segmen yang dipakai di ujung bulat. Alat yang berbeda bisa selisih sekitar 0,01 hektare. Itu bukan galat.

Latihan variasi. Ubah sempadan menjadi 50 m. Contoh jawaban, ganti baris 14 dengan:

arcpy.analysis.Buffer("Sungai.shp", r"memory\sempadan", "50 Meters", dissolve_option="ALL")

Dari data contoh, baris pertama sampai ketiga menjadi sempadan kira-kira 1.69, 2.32, dan 2.54 hektare. Baris NDVI tidak berubah. Perkiraan ini dihitung dengan OGR. [CEK]

Skrip 31.3: Batch NDVI semua ubin dan potongan CHM per KPH #

Berkas b3_arcpy_pro.py dan b3_arcpy_arcmap.py. Di dalam berkasnya tertulis "Skrip 3.1". Padanan PyQGIS-nya ada di Bab 21.

Tujuan. Menghitung NDVI untuk semua ubin foto drone dalam satu folder, lalu memotong peta tinggi tajuk (CHM) satu per satu untuk tiap KPH.

Bahan. Paket data paket-i4: folder Ortofoto_Mentah (tiga ubin; band 1 inframerah dekat, band 2 merah), DSM_Drone.tif, DTM_Drone.tif, dan shp\Batas_KPH.shp. Berbeda dari Skrip 31.2, skrip ini membuat folder hasil sendiri. Ekstensi Spatial Analyst harus tersedia. CHM (canopy height model) adalah tinggi tajuk pohon: permukaan tertinggi (DSM) dikurangi permukaan tanah (DTM).

Versi ArcGIS Pro (Python 3) #

# [SKRIP 3.1-Pro: Batch: NDVI semua ubin dan potongan CHM per KPH dengan ArcPy]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh ekstensi Spatial Analyst.
import os
import arcpy
from arcpy.sa import Raster, Float, ExtractByMask

DATA = r"D:\Latihan\paket-i4"
HASIL = r"D:\Latihan\hasil_i4"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
os.makedirs(HASIL + r"\NDVI_ubin", exist_ok=True)

# A. Semua ubin dalam satu folder (band 1 = NIR, band 2 = Merah)
arcpy.env.workspace = DATA + r"\Ortofoto_Mentah"
daftar = arcpy.ListRasters("*", "TIF")
for i, nama in enumerate(daftar, start=1):
    b1 = Float(Raster(os.path.join(DATA, "Ortofoto_Mentah", nama, "Band_1")))
    b2 = Float(Raster(os.path.join(DATA, "Ortofoto_Mentah", nama, "Band_2")))
    ((b1 - b2) / (b1 + b2)).save(os.path.join(HASIL, "NDVI_ubin", "NDVI_" + nama))
    print("%d/%d selesai: %s" % (i, len(daftar), nama))

# B. Satu putaran per KPH: pilih satu poligon, potong CHM dengannya
chm = Raster(DATA + r"\DSM_Drone.tif") - Raster(DATA + r"\DTM_Drone.tif")
batas = DATA + r"\shp\Batas_KPH.shp"
with arcpy.da.SearchCursor(batas, ["NAMA_KPH"]) as k:
    nama_kph = [b[0] for b in k]
for nama in nama_kph:
    arcpy.analysis.Select(batas, r"memory\satu_kph", "NAMA_KPH = '%s'" % nama)
    keluar = os.path.join(HASIL, "CHM_" + nama.replace(" ", "_") + ".tif")
    ExtractByMask(chm, r"memory\satu_kph").save(keluar)
    print("CHM tersimpan untuk", nama)
arcpy.CheckInExtension("Spatial")

Penjelasan baris demi baris. Hanya bagian yang baru dijelaskan lengkap. Baris yang mirip Skrip 30.3 dan 31.2 dirangkum.

Versi ArcMap 10.8 (Python 2.7) #

# -*- coding: utf-8 -*-
# [SKRIP 3.1-ArcMap: Batch: NDVI semua ubin dan potongan CHM per KPH (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh ekstensi Spatial Analyst.
import os
import arcpy
from arcpy.sa import Raster, Float, ExtractByMask

DATA = r"D:\Latihan\paket-i4"
HASIL = r"D:\Latihan\hasil_i4"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
if not os.path.exists(HASIL + r"\NDVI_ubin"):          # Python 2.7 tidak punya exist_ok
    os.makedirs(HASIL + r"\NDVI_ubin")

# A. Semua ubin dalam satu folder (band 1 = NIR, band 2 = Merah)
arcpy.env.workspace = DATA + r"\Ortofoto_Mentah"
daftar = arcpy.ListRasters("*", "TIF")
for i, nama in enumerate(daftar, 1):
    b1 = Float(Raster(os.path.join(DATA, "Ortofoto_Mentah", nama, "Band_1")))
    b2 = Float(Raster(os.path.join(DATA, "Ortofoto_Mentah", nama, "Band_2")))
    ((b1 - b2) / (b1 + b2)).save(os.path.join(HASIL, "NDVI_ubin", "NDVI_" + nama))
    print("%d/%d selesai: %s" % (i, len(daftar), nama))

# B. Satu putaran per KPH: pilih satu poligon, potong CHM dengannya
chm = Raster(DATA + r"\DSM_Drone.tif") - Raster(DATA + r"\DTM_Drone.tif")
batas = DATA + r"\shp\Batas_KPH.shp"
with arcpy.da.SearchCursor(batas, ["NAMA_KPH"]) as k:
    nama_kph = [b[0] for b in k]
for nama in nama_kph:
    arcpy.Select_analysis(batas, r"in_memory\satu_kph", "\"NAMA_KPH\" = '%s'" % nama)
    keluar = os.path.join(HASIL, "CHM_" + nama.replace(" ", "_") + ".tif")
    ExtractByMask(chm, r"in_memory\satu_kph").save(keluar)
    print("CHM tersimpan untuk %s" % nama)
arcpy.CheckInExtension("Spatial")

Hasil. Layar menampilkan enam baris. Folder hasil berisi tiga berkas NDVI di NDVI_ubin dan tiga berkas CHM. Berikut perkiraannya. Nama ubin dan KPH mengikuti paket data. [CEK]

1/3 selesai: ubin_01.tif
2/3 selesai: ubin_02.tif
3/3 selesai: ubin_03.tif
CHM tersimpan untuk KPH Alpha
CHM tersimpan untuk KPH Beta
CHM tersimpan untuk KPH Gamma

Berkas yang terbentuk: NDVI_ubin\NDVI_ubin_01.tif sampai NDVI_ubin_03.tif, dan CHM_KPH_Alpha.tif, CHM_KPH_Beta.tif, CHM_KPH_Gamma.tif langsung di folder hasil. Dengan GDAL dari data contoh, rata-rata NDVI ketiga ubin kira-kira 0,471, 0,471, dan 0,467. Rata-rata tinggi tajuk per KPH kira-kira 2,60 m, 2,69 m, dan 2,56 m, sama dengan hasil skrip PyQGIS di Bab 21. [CEK] Logika Python skrip sudah diuji dengan modul ArcPy tiruan buatan sendiri, dan cetakannya sama dengan perkiraan di atas. Perilaku ArcGIS yang sebenarnya belum teruji.

Bagian yang boleh diubah. Nomor baris mengikuti versi Pro.

• Baris 8-9, lokasi. Ganti sesuai komputer Anda.
• Baris 16, pola dan format. Ganti "*" atau "TIF" untuk memilih ubin lain.
• Baris 18-19, nomor band. Tukar Band_1 dan Band_2 bila urutan band berbeda.
• Baris 29, klausa pilih. Ganti NAMA_KPH dengan kolom pembeda lain.
• Baris 30, awalan nama keluaran. Ganti "CHM_" untuk awalan lain.

Jebakan.

• Band terbalik. Bila band 1 sebenarnya merah, NDVI negatif dan tidak ada peringatan.
• Kutip dalam klausa pilih. Pro memakai NAMA_KPH = '%s' tanpa kutip ganda, ArcMap memakai kutip ganda di nama kolom. Dokumentasi Esri untuk Select memberi contoh dengan nama kolom diapit kutip ganda. Mana yang perlu untuk Shapefile bergantung pada versi. [CEK]
• Nama KPH yang memuat tanda kutip tunggal. Teks 'KPH D'Arc' merusak klausa pilih, karena kutip tunggal di tengah nama dianggap penutup. [CEK]
• Objek sementara dipakai ulang. memory\satu_kph ditulis ulang tiap putaran. Itu aman karena overwriteOutput = True.
• Lokasi CHM berbeda dari versi QGIS. Di sini CHM disimpan langsung di folder hasil. Versi PyQGIS di Bab 21 memakai subfolder CHM_per_KPH.

Latihan variasi. Potong CHM hanya untuk KPH Beta. Contoh jawaban, tambahkan satu baris sebelum baris 28:

nama_kph = [n for n in nama_kph if n == "KPH Beta"]

Daftar disaring sehingga hanya berisi KPH Beta. Layar menampilkan tiga baris kemajuan ubin, lalu CHM tersimpan untuk KPH Beta. Jawaban ini sudah diuji dengan modul ArcPy tiruan. [CEK]

Skrip 31.4: Menangkap galat dengan ArcPy #

Berkas b4_arcpy_pro.py dan b4_arcpy_arcmap.py. Di dalam berkasnya tertulis "Skrip 4.1". Padanan PyQGIS-nya ada di Bab 22.

Tujuan. Memeriksa tiap kiriman Shapefile dari tim lapangan: apakah datanya ada, apakah kolom wajib lengkap, dan apakah sistem koordinatnya sudah EPSG 32749. Yang berbeda diproyeksikan, dan semua hasil dicatat dalam satu daftar log.

Bahan. Paket data paket-i4, folder Kiriman_Lapangan\shp berisi enam kiriman, tim_a.shp sampai tim_f.shp. Satu kiriman (tim_d) memakai sistem koordinat derajat (EPSG 4326), satu (tim_e) berkas rusak, dan satu (tim_f) tanpa kolom Kondisi. Skrip ini hanya menguji, tidak menyimpan hasil proyeksi. Tidak perlu Spatial Analyst.

Lima istilah baru, semuanya dibahas di Bab 7 dan 10:

• def: membuat fungsi, yaitu potongan kerja yang diberi nama.
• raise: sengaja memunculkan galat dengan pesan yang jelas.
• try dan except: "coba kerjakan ini; bila galat jenis tertentu, kerjakan itu".
• arcpy.ExecuteError: jenis galat yang dimunculkan ArcPy ketika sebuah alat gagal.
• glob: perkakas untuk mencari berkas dengan pola nama, misalnya semua *.shp.

Versi ArcGIS Pro (Python 3) #

# [SKRIP 4.1-Pro: Penanganan galat dengan ArcPy: periksa masukan lalu tangkap ExecuteError]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Kiriman lapangan berbentuk Shapefile di folder Kiriman_Lapangan\shp.
import glob
import os
import sys
import arcpy

DATA = r"D:\Latihan\paket-i4"
FIELD_WAJIB = ["Kondisi", "Tinggi_Phn", "Foto"]
log = []


def periksa(jalur):
    """Berhenti dengan pesan jelas bila masukan tidak layak. Mengembalikan kode EPSG data."""
    if not arcpy.Exists(jalur):
        raise ValueError("data tidak ada atau tidak bisa dibuka")
    nama_field = [f.name for f in arcpy.ListFields(jalur)]
    kurang = [f for f in FIELD_WAJIB if f not in nama_field]
    if kurang:
        raise ValueError("field hilang: " + ", ".join(kurang))
    return arcpy.Describe(jalur).spatialReference.factoryCode


for jalur in sorted(glob.glob(DATA + r"\Kiriman_Lapangan\shp\*.shp")):
    nama = os.path.basename(jalur)
    try:
        epsg = periksa(jalur)
        if epsg != 32749:
            arcpy.management.Project(jalur, "memory\\" + nama[:-4], arcpy.SpatialReference(32749))
            log.append((nama, "DITERIMA", "CRS %s diubah ke 32749" % epsg))
        else:
            log.append((nama, "DITERIMA", "sesuai"))
    except ValueError as e:                      # galat milik skrip kita
        log.append((nama, "DITOLAK", str(e)))
    except arcpy.ExecuteError:                   # galat dari alat geoprocessing
        log.append((nama, "GAGAL ALAT", arcpy.GetMessages(2)))
    except Exception:                            # galat Python lain
        log.append((nama, "TAK TERDUGA", str(sys.exc_info()[1])))
for baris in log:
    print("%-14s %-12s %s" % baris)

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

Dokumentasi Esri menyebut: alat yang gagal memunculkan arcpy.ExecuteError, GetMessages(2) mengambil pesan berderajat galat, dan blok try dengan beberapa except dipakai untuk membedakan galat alat dari galat Python lain. Ketiga hal itu cocok dengan baris 36-39. [CEK]

Versi ArcMap 10.8 (Python 2.7) #

# -*- coding: utf-8 -*-
# [SKRIP 4.1-ArcMap: Penanganan galat dengan ArcPy (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Kiriman lapangan diandaikan berbentuk Shapefile di folder Kiriman_Lapangan\shp.
import glob
import os
import sys
import arcpy

DATA = r"D:\Latihan\paket-i4"
FIELD_WAJIB = ["Kondisi", "Tinggi_Phn", "Foto"]
log = []


def periksa(jalur):
    """Berhenti dengan pesan jelas bila masukan tidak layak. Mengembalikan kode EPSG data."""
    if not arcpy.Exists(jalur):
        raise ValueError("data tidak ada atau tidak bisa dibuka")
    nama_field = [f.name for f in arcpy.ListFields(jalur)]
    kurang = [f for f in FIELD_WAJIB if f not in nama_field]
    if kurang:
        raise ValueError("field hilang: " + ", ".join(kurang))
    return arcpy.Describe(jalur).spatialReference.factoryCode


for jalur in sorted(glob.glob(DATA + r"\Kiriman_Lapangan\shp\*.shp")):
    nama = os.path.basename(jalur)
    try:
        epsg = periksa(jalur)
        if epsg != 32749:
            arcpy.Project_management(jalur, "in_memory\\" + nama[:-4], arcpy.SpatialReference(32749))
            log.append((nama, "DITERIMA", "CRS %s diubah ke 32749" % epsg))
        else:
            log.append((nama, "DITERIMA", "sesuai"))
    except ValueError as e:                      # galat milik skrip kita
        log.append((nama, "DITOLAK", str(e)))
    except arcpy.ExecuteError:                   # galat dari alat geoprocessing
        log.append((nama, "GAGAL ALAT", arcpy.GetMessages(2)))
    except Exception:                            # galat Python lain
        log.append((nama, "TAK TERDUGA", str(sys.exc_info()[1])))
for baris in log:
    print("%-14s %-12s %s" % baris)

Hasil. Layar menampilkan satu baris untuk tiap kiriman. Berikut perkiraannya. Dari enam kiriman, empat diterima dan dua ditolak, sama dengan hasil skrip PyQGIS di Bab 22. [CEK]

tim_a.shp      DITERIMA     sesuai
tim_b.shp      DITERIMA     sesuai
tim_c.shp      DITERIMA     sesuai
tim_d.shp      DITERIMA     CRS 4326 diubah ke 32749
tim_e.shp      DITOLAK      data tidak ada atau tidak bisa dibuka
tim_f.shp      DITOLAK      field hilang: Kondisi

Perkiraan ini disusun dengan membaca sistem koordinat dan kolom tiap Shapefile lewat OGR, lalu menjalankan logika Python skrip dengan modul ArcPy tiruan. Dua hal bisa berbeda di ArcGIS sungguhan. Pertama, kode sistem koordinat yang dibaca dari berkas .prj. Kedua, apakah arcpy.Exists menganggap tim_e.shp (berkas rusak) ada atau tidak. Bila Exists mengatakan ada, kiriman itu tertangkap di cabang lain. [CEK]

Bagian yang boleh diubah. Nomor baris mengikuti versi Pro.

• Baris 9, lokasi data. Ganti sesuai komputer Anda.
• Baris 10, kolom wajib. Tambah atau kurangi nama kolom di daftar FIELD_WAJIB.
• Baris 25, folder dan pola. Ganti Kiriman_Lapangan\shp dan *.shp untuk folder atau format lain.
• Baris 29-31, kode tujuan. Ganti 32749 (tiga tempat, termasuk teks di baris 31) bila sistem koordinat tujuan Anda lain.

Jebakan.

• Urutan except penting. Python memeriksa dari atas. Bila except Exception: diletakkan paling atas, ia menangkap semuanya, dan cabang ValueError serta ExecuteError tidak pernah dipakai. Semua kiriman akan berlabel TAK TERDUGA. Perilaku ini sudah diuji di Python 3.12. Tulis penangkap yang paling khusus di atas, yang paling umum di bawah.
• Hasil proyeksi tidak disimpan. Keluarannya ke memory\ dan tidak dipakai lagi. Skrip hanya menguji kelayakan. Berbeda dari versi PyQGIS di Bab 22, tidak ada penggabungan data dan tidak ada pemeriksaan titik di luar KPH.
• Data tanpa sistem koordinat. Bila factoryCode bernilai 0, baris 29 memicu Project yang gagal. Kasusnya tertangkap sebagai GAGAL ALAT. [CEK]
• Pesan galat tidak diperiksa lebih jauh. Cabang terakhir hanya mencatat teksnya. Untuk menyelidiki, cetak seluruh isi galat.
• Kolom wajib dibaca dari Shapefile. Nama kolom Shapefile maksimal 10 karakter. Tinggi_Phn pas 10 karakter, jadi aman. Nama yang lebih panjang terpotong. [CEK]

Latihan variasi. Tambahkan kolom wajib Catatan, yang tidak ada di kiriman mana pun. Contoh jawaban, ganti baris 10 dengan:

FIELD_WAJIB = ["Kondisi", "Tinggi_Phn", "Foto", "Catatan"]

Dari data contoh, keempat kiriman yang terbaca (tim_a sampai tim_d) berubah menjadi DITOLAK dengan keterangan field hilang: Catatan. Kiriman tim_f ditolak dengan keterangan field hilang: Kondisi, Catatan. Kiriman tim_e tetap ditolak karena datanya tidak terbuka. Jawaban ini sudah diuji dengan modul ArcPy tiruan. [CEK]

Skrip 31.5: Proyek mini terpadu, ringkasan per KPH ke CSV #

Berkas b7_arcpy_pro.py dan b7_arcpy_arcmap.py. Di dalam berkasnya tertulis "Skrip 7.1". Padanan PyQGIS-nya ada di Bab 23.

Tujuan. Menyatukan alat-alat sebelumnya menjadi satu proyek kecil: menghitung luas, sempadan sungai, titik panas, NDVI rata-rata, dan tinggi tajuk rata-rata per KPH, lalu menulisnya ke berkas CSV.

Bahan. Paket data paket-i4: folder shp, Citra_KPH.tif, DSM_Drone.tif, dan DTM_Drone.tif. Siapkan folder D:\Latihan\hasil_i4 lebih dulu: skrip tidak membuatnya. Ekstensi Spatial Analyst harus tersedia. Skrip ini memakai Shapefile, supaya sejajar dengan bagian QGIS, menurut komentar di baris 3.

Versi ArcGIS Pro (Python 3) #

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

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

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

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

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

Penjelasan baris demi baris. Bagian alat dan raster sudah Anda kenal di Skrip 31.2. Yang baru adalah try/finally, fungsi di dalam try, dan penulisan CSV.

Dokumentasi Esri mencatat pola try/except/finally untuk galat alat, dan memberi contoh arcpy.AddError(arcpy.GetMessages(2)) pada cabang except. Skrip ini mencetak pesannya, bukan memakai AddError. [CEK]

Versi ArcMap 10.8 (Python 2.7) #

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

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

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

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

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

Hasil. Layar menampilkan satu baris, dan berkas CSV terbentuk di folder hasil. Berikut perkiraannya. Angka dihitung dengan GDAL dan OGR dari data contoh, dan sama dengan hasil skrip PyQGIS di Bab 23. Isi CSV tampak seperti ini, dengan angka ditulis apa adanya oleh csv.writer tanpa nol di belakang. [CEK]

Ringkasan tersimpan.

NAMA_KPH,LUAS_HA,SEMPADAN_HA,HOTSPOT,NDVI_RATA,CHM_RATA_M
KPH Alpha,4.0,1.18,2,0.604,2.6
KPH Beta,4.0,1.39,5,0.431,2.69
KPH Gamma,8.0,1.5,5,0.329,2.56

Ringkasan ArcPy punya enam kolom. Ringkasan di Bab 23 punya satu kolom lebih banyak, yaitu TITIK_SURVEI. Logika Python skrip sudah diuji dengan modul ArcPy tiruan buatan sendiri, dan isi CSV di atas adalah keluarannya. Perilaku alat ArcGIS yang sebenarnya belum teruji.

Bagian yang boleh diubah. Nomor baris mengikuti versi Pro.

• Baris 8-9, lokasi. Ganti sesuai komputer Anda.
• Baris 16, jarak sempadan. Ganti "30 Meters".
• Baris 37, nama berkas CSV. Ganti ringkasan_kph_pro.csv.
• Baris 39, judul kolom. Ubah daftar judul, lalu ubah isi baris 42-43 agar urutannya sama.

Jebakan.

• Folder hasil tidak dibuat. Bila D:\Latihan\hasil_i4 belum ada, penulisan CSV gagal. [CEK]
• Nama KPH tidak ada di tabel zonal. Baris 43 memakai n_rata[n] tanpa .get. Nama yang tidak ada memicu galat KeyError, misalnya KeyError: 'KPH Delta'. Galat ini tidak ditangkap oleh except arcpy.ExecuteError, jadi skrip berhenti dengan pesan galat Python. Pesan KeyError sudah diuji pada kamus biasa di Python 3.
• Penyandian CSV. Berkas ditulis tanpa encoding. Di Windows, Python memakai penyandian bawaan komputer. Nama KPH di contoh hanya huruf ASCII, jadi aman. [CEK]
• newline="" di Python 3. Dokumentasi modul csv menganjurkan pembukaan dengan newline="" agar tidak muncul baris kosong ganda di Windows. [CEK]
• Pembagian dengan nol pada latihan. Bila luas sebuah KPH 0, hitungan persen gagal.

Latihan variasi. Tambahkan kolom SEMPADAN_PERSEN, yaitu persen luas sempadan terhadap luas KPH. Contoh jawaban: ganti baris 39 sampai 43 dengan:

        w.writerow(["NAMA_KPH", "LUAS_HA", "SEMPADAN_HA", "HOTSPOT", "NDVI_RATA", "CHM_RATA_M", "SEMPADAN_PERSEN"])
        with arcpy.da.SearchCursor(r"memory\ringkas", ["NAMA_KPH", "SHAPE@AREA", "Join_Count"]) as k:
            for n, luas, jml in k:
                w.writerow([n, round(luas / 10000, 2), round(sempadan.get(n, 0), 2), jml,
                            round(n_rata[n], 3), round(c_rata[n], 2),
                            round(100 * sempadan.get(n, 0) / (luas / 10000), 1)])

Dari data contoh, kolom baru bernilai kira-kira 29.4, 34.7, dan 18.7. Skrip PyQGIS di Bab 21 menampilkan 34.6 untuk KPH Beta. Selisih 0,1 berasal dari perbedaan jumlah segmen pada sempadan di dua alat. Jawaban ini sudah diuji dengan modul ArcPy tiruan. [CEK]

Cek paham #

1. Mengapa skrip ArcMap mencetak dengan satu teks hasil rakitan %, bukan beberapa isian dipisah koma?
2. Apa fungsi arcpy.env.workspace?
3. Mengapa bagian 4 Skrip 31.1 gagal bila dijalankan dari terminal biasa?
4. Apa beda memory\sempadan dan in_memory\sempadan?
5. Mengapa except arcpy.ExecuteError harus ditulis di atas except Exception?
6. Apa fungsi blok finally di Skrip 31.5?

Jawaban.

1. Karena print pada Python 2.7 berperilaku lain bila diberi beberapa isian. Satu teks hasil rakitan aman di kedua versi.
2. Menetapkan folder kerja. Nama data yang ditulis tanpa lokasi dicari di folder itu.
3. Karena CURRENT berarti proyek atau dokumen peta yang sedang terbuka. Di terminal biasa tidak ada proyek yang terbuka.
4. Keduanya ruang sementara di memori komputer. memory dipakai di ArcGIS Pro, sedangkan in_memory adalah nama gaya lama yang dipakai di ArcMap.
5. Python memeriksa penangkap dari atas. Exception mencakup semua jenis galat, jadi bila ditulis lebih dulu, penangkap yang lebih khusus tidak pernah dipakai.
6. Menjalankan bagian itu selalu, berhasil atau galat. Di skrip itu, ia mengembalikan ekstensi Spatial Analyst.

Kesalahan umum #

• Lupa huruf r pada lokasi data. Perbaikan: tulis r"D:\Latihan\paket-i4".
• Memakai gaya Python 3 di ArcMap. Perbaikan: ikuti versi ArcMap di bab ini, terutama untuk print, exist_ok, dan modul arcpy.mapping.
• Menganggap luas selalu hektare. Perbaikan: cek satuan peta lewat sistem koordinat data.
• Lupa membuat folder hasil. Skrip 31.2 dan 31.5 tidak membuatnya. Perbaikan: buat foldernya dulu, atau tambahkan os.makedirs seperti di Skrip 31.3.
• Menaruh except Exception di paling atas. Perbaikan: urutkan dari yang paling khusus ke yang paling umum.

Ringkasan dan latihan #

• ArcPy adalah perkakas Python untuk ArcGIS. ArcGIS Pro memakai Python 3, ArcMap 10.8 memakai Python 2.7.
• Pola skrip: lokasi data, sifat data, kursor per baris, lalu peta. Untuk rantai alat: tiap alat menulis ke ruang sementara.
• Lima skrip I4 di bab ini: membaca data dan proyek, rantai alat, batch, penanganan galat, dan proyek terpadu.
• Beda ArcMap yang utama: print dengan %, 10000.0, nama alat gaya _analysis, in_memory, arcpy.mapping, exist_ok yang tidak ada, dan mode "wb" untuk CSV.
• Skrip bagian ini hanya bisa dijalankan di Windows, dan bertanda [CEK].

Latihan. Ganti data pada Skrip 31.1 menjadi Sungai.shp. Cetak jumlah fitur dan daftar kolomnya saja. Bagian mana dari skrip yang boleh Anda hapus? Lalu, di Skrip 31.4, tambahkan catatan jumlah kiriman yang diterima di akhir log. Cukup satu baris print.

Tabel perbandingan: QGIS dan ArcPy #