Peta Gratis Dengan Mengunduh Fitur Dari ArcGIS Rest Service (Bagian 2) – Script Python

Pada tulisan bagian 1 telah dijelaskan cara mengekstrak peta atau data spasial dari situs geoportal berbasis ArcGIS Rest Service dengan cara memanfaatkan fasilitas query nya secara langsung, nah pada tulisan bagian 2 ini Lintas Bumi akan menjelaskan alternatif lainnya yaitu dengan memanfaatkan script python.

Python merupakan salah satu bahasa pemrograman dasar yang menjadi fondasi software SIG seperti ArcGIS, QGIS, dan lain-lain. Python adalah bahasa pemrograman tingkat tinggi yang ditafsirkan untuk pemrograman tujuan umum. Dibuat oleh Guido van Rossum dan pertama kali dirilis pada tahun 1991, Python memiliki filosofi desain yang menekankan keterbacaan kode, terutama menggunakan spasi yang signifikan. Ini menyediakan konstruksi yang memungkinkan pemrograman yang jelas pada skala kecil dan besar. Di dalam ArcGIS bahkan ESRI membuat library khusus yang disebut ArcGIS API for Python.

Untuk metode kedua ini script python yang digunakan sebaiknya dijalankan melalui program Python yang menjadi bawaan ArcGIS, jika tidak anda harus menginstall Python beserta library yang dibutuhkan secara manual. Script ini Lintas Bumi kutip dari situs spatialtimes.com (thanks bro!). Script ini untuk ‘mengakali’ fitur yang jumlahnya lebih dari 1000, di mana dengan script ini proses unduhnya akan menjadi berkali-kali atau berulang setiap 1000 fitur. Lintas Bumi telah mencoba script ini di situs http://geoportal.menlhk.go.id/arcgis/rest/services milik KemenLHK.

Tambahan ‘senjata’ yang harus disiapkan untuk metode kedua ini selain yang sudah dijelaskan pada tulisan bagian 1 adalah;

  • Pastikan anda sudah menginstall arcgis dan python (kalau bisa minimal versi 2.7). Jika anda sudah menginstall ArcGIS biasanya ada di folder c:/Python27.
  • Jika tidak silahkan anda install sendiri Python dengan versi >= 2.7 beserta library yang diperlukan.
  • Software text editor (notepad, notepad++, dll.) – Lintas Bumi menggunakan notepad++
  • Sedikit pemahaman tentang Python.

 

#Langkah pertama: mengunjungi situs

Kunjungi situs yang akan dituju, pada contoh ini situs milik KLHK di atas, lalu klik Folder KLHK dan pada penjelasan kali ini Lintas Bumi bermaksud mengunduh peta lokasi kebakaran hutan tahun 2017 sehingga dipilih (klik) KLHK/BURN_AREA_2017 (MapServer) danpada Layer dipilih lagi Burn Area Indonesia 2017 (0). Pastikan layer yang diinginkan mempunyai format json (lihat bagian kiri atas gambar berikut) atau info di halaman tersebut.

 

#Langkah kedua: cek fungsi query

Jangan lupa cek fungsi query dari layer yang bersangkutan, lihat di bagian bawah halaman bahwa Supported Operations : Query aktif.

 

#Langkah ketiga: memastikan jumlah fitur layer

Sama dengan tulisan bagian 1, sebelum melakukan query anda harus memastikan terlebih dulu jumlah fitur yang ada di dalam Layer KLHK/BURN_AREA_2017, caranya klik Query di bawah tadi, masukan parameter Where: 1=1, Return Geomtery: False, dan Return Count Only: True.

Klik Query (GET), biarkan proses beberapa saat (tunggu), sampai muncul di bawahnya jumlah fitur (Count), sebagai contoh jumlahnya 2737.

 

#Langkah keempat: gunakan script python

Proses selanjutnya adalah mengekstrak data menggunakan Python, dengan script ini sebetulnya proses query dirubah dan dijalankan dalam kode Python, yang bersifat iterasi atau pengulangan, di mana proses query akan membuat file bertipe .json dan .shp per 1000 feature. Sehingga layer di atas dengan jumlah 2737 fitur, hasil unduhannya akan dibuat menjadi 3 buah file .shp yaitu 2 .shp pertama mengandung masing-masing 1000 fitur dan shp ketiga 737 fitur. Semua file .shp hasil unduh ini sebaiknya digabungkan secara manual. Silahkan copy kan dan edit kode Python berikut pada notepad++ (dengan Run As Administrator), modifikasi baris ke 11, 12, 13, dan 17.

#Name: Export ArcGIS Server Map Service Layer to Shapefile with Iterate
#Author: Bryan McIntosh
#Description: Python script that connects to an ArcGIS Server Map Service and downloads a single vector layer to shapefiles. If there are more features than AGS max allowed, it will iterate to extract all features.

import urllib2,json,os,arcpy,itertools
ws = os.getcwd() + os.sep

#Set koneksi ke ArcGIS Server, map service, layer ID, dan jumlah maksimal request ke server (1000 adalah default jika tidak diketahui). 
#Contoh alamat koneksi ke geoportal kemenlhk adalah http://geoportal.menlhk.go.id/arcgis/rest/services/KLHK/KHG/MapServer/0/query, maka dibagi menjadi 3, tulisan query di belakang tidak disertakan, sehingga menjadi;

serviceURL = "http://geoportal.menlhk.go.id/arcgis/rest/services"
serviceMap = "/KLHK/BURN_AREA_2017/MapServer"
serviceLayerID = 0
serviceMaxRequest = 1000

#set nama file json dan shp output
dataOutputName = "Burn2017"

def defServiceGetIDs():
    IDsRequest = serviceURL + serviceMap + "/" + str(serviceLayerID) + "/query?where=1%3D1&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=true&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&f=pjson"
    IDsResponse = urllib2.urlopen(IDsRequest)
    IDsJSON = json.loads(IDsResponse.read())
    IDsSorted = sorted(IDsJSON['objectIds'])
    return IDsSorted

def defGroupList(n, iterable):
    args = [iter(iterable)] * n
    return ([e for e in t if e != None] for t in itertools.izip_longest(*args))
    
def defQueryExtractRequests(idMin, idMax):
    myQuery = "&where=objectid+>%3D+" + idMin + "+and+objectid+<%3D+" + idMax
    myParams = "query?geometryType=esriGeometryEnvelope&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&returnZ=false&returnM=false&returnDistinctValues=false&returnTrueCurves=false&f=pjson"
    myRequest = serviceURL + serviceMap + "/" + str(serviceLayerID) + "/" + myParams + myQuery
    response = urllib2.urlopen(myRequest)
    myJSON = response.read()
  
    # Write response to json text file
    foo = open(dataOutputName + idMin + ".json", "w+")
    foo.write(myJSON);
    foo.close()
  
    # Create Feature Class
    arcpy.JSONToFeatures_conversion(dataOutputName + idMin + ".json", ws + dataOutputName + idMin + ".shp")
    
#**MAIN**#
#Get all objectIDs (OIDs) for the layer (there is no server limit for this request)
AllObjectIDs = defServiceGetIDs()

#Divide the OIDs into chunks since there is a limit to map queries (assumed limit stored in serviceMaxRequest variable)
ObjectID_Groups = list(defGroupList(serviceMaxRequest, AllObjectIDs))

#Membuat shapefile for masing-masing geojson
for ObjectID_Group in ObjectID_Groups:
    idMin = str(ObjectID_Group[0])
    idMax = str(ObjectID_Group[-1])
    defQueryExtractRequests(idMin, idMax)
  
#Gabungkan shp hasil dengan perintah Append di aplikasi SIG

Jika kode sudah disesuaikan, silahkan simpan (save) sebagai unduhpeta.py atau nama lain dengan extensi .py (Python), dan sangat disarankan simpan di folder Python ArcGIS sebagai contoh di c:/Program Files/Python27 atau c:/Python27.

 

#Langkah kelima: jalankan kode python dan unduh

Setelah dipastikan kodenya ok dan lengkap, silahkan anda jalankan command prompt windows (cari dengan kata ‘cmd’ di Windows) dan masuk ke folder Python ArcGIS, jalankan perintah (ketikan) python unduhpeta.py. Tunggu beberapa saat, lama proses bergantung jumlah fitur yang diunduh.

Jika proses sudah selesai, tidak ada pesan error dan command prompt akan kembali ke kondisi semula. Lama proses unduh untuk 2737 feature dengan koneksi 4G yang dicoba Lintas Bumi adalah sekitar 4 menit.

 

#Langkah keenam: buka file hasil unduh

Berikut adalah hasil unduhan, terlihat ada 3 file .json, dan 3 file .shp seperti yang diperkirakan. Berarti unduhan sukses dan silahkan anda buka file .shp di software SIG dan anda bisa menggabungkannya menjadi satu shp baru sesuai keinginan anda.

Lintas Bumi mencoba membuka file hasil unduh di software QGIS atau anda bisa membukanya di software SIG lainnya, dan berikut adalah contoh .shp yang ketiga beserta tabel atributnya yang berjumlah 737.

Hasil pengetesan pada beberapa layer berukuran besar terkadang masih gagal karena time out (lebih dari 60 detik), namun tak ada salahnya selalu berusaha. Selamat mencoba, semoga sukses.

Bagikan ini:

Tinggalkan Balasan

Alamat email Anda tidak akan dipublikasikan. Ruas yang wajib ditandai *