Membuat grid/Fishnet dengan Python

Saat akan melakukan survei ataupun analisis terkadang memerlukan grid. Grid atau Fishnet dapat dibuat dengan menggunakan python dengan menggunakan library json, geojson, numpy, dan shapely

Mengimport library

sebelum menggunakan library maka perlu untuk membuka terlebih dahulu librarynya dengan cara mengimportnya, apabila belum terinstal maka perlu untuk menginstall terlebih dahulu dengan menggunkana perintah ‘pip install <nama library>’ berikut kode untuk mengimpor library yang dibutuhkan

import geojson
from shapely.geometry import shape
import numpy as np
import json

Membuka file dan menentukan parameter

kita membutuhkan batas untuk pembuatan grid, misalnya disini menggunakan batas berupa polygon seperti gambar berikut

untuk membuka polygon dengan ekstensi geojson dan menentukan panjang dan lebar grid dapat menggunakan kode berikut:

#membuka file geojson yang bernama 'input.geojson'
with open('input.geojson', 'r') as f:
  input = json.loads(f.read())

#menentukan parameter yaitu panjang x dan panjang y dengan satuan meter
xresm = 1000
yresm = 1000

Menentukan Bounding Box

Sebelum membuat grid perlu untuk menentukan batas kiri atas dan kanan bawah yang mencakup seluruh polygon input yang akan digunakan untuk membuat grid, dengan kode sebgai berikut:

input = geojson_input['features'][0]['geometry']
box = []
for i in (0,1): #0 adalah longitude berdasarkan urutan dalam tuple, 1 adalah latitude berdasarkan urutan dalam tuple
    #res = sorted(list(geojson.utils.coords(input)), key=lambda x:x[i])
    res = sorted(list(geojson.utils.coords(input)), key=lambda x:x[i])
    #sorting longitude dan latitude dari kecil ke besar
    box.append((res[0][i],res[-1][i])) #diambil nilai terkecil dan terbesar sebagai bounding box
ret = f"({box[0][0]},{box[1][1]},{box[0][1]},{box[1][0]})" #format bounding box menjadi x1,y1,x2,x3

Memecah bounding box menjadi koordinat terpisah

koordinat dipecah nilai koordiant latitude atas, longitude kiri, latitude bawah, dan longitude kanan

xy = ret[1:-1].split(',')
ulx, uly, lrx, lry = float(box[0][0]),float(box[1][1]),float(box[0][1]),float(box[1][0])

Mengkonversi nilai panjang dan lebar

panjang dan lebar yang ditentukan sebelumnya dalam satuan meter maka perlu dikonversi ke dalam satuan derajat karena input koordinat polygon merupakan koordinat geografis. Selain mengkonversi nilai panjang dan lebar juga perlu untuk menghitung setengah dari nilai panjang dan lebar untuk sebagai variabel perhitungan selanjutnya, dengan kode sebagai berikut:

# resolution 1 degree grid
xres = 0.00001*(xresm/1.11)
yres = -0.00001*(yresm/1.11)

# setengan dari resolusi untuk mencari koordinat tengah setiap kotak grid
dx = xres/2
dy = yres/2

Membuat grid dengan library numpy

pembuatan grid dengan menggunakan numpy meshgrid dan dengan input berupa numpy array, sehingga variabel yang telah ditentukan sebelumnya disusun dengan menggunakan numpy arrange. Dengan kode sebagai berikut:

# membuat grid dengan numpy meshgrid
xx, yy = np.meshgrid(
    np.arange(ulx+dx, lrx+dx, xres), #dibuat numpy array
    np.arange(uly+dy, lry+dy, yres), #dibuat numpy array
)

Memformat hasil numpy meshgrid menjadi geojson

Hasil dari numpy mesgrid berupa array seluruh koordinat dari x dan y, sehingga perlu disusun menjadi koordinat grid, dimana hasil dari numpy meshgrid ini masih berupa titik tengah dari kotak grid, maka perlu menggunakan nilai dari setengah resolusi tadi untuk menghitung koordinat grid. Dengan kode sebagai berikut:

grid = {"type":"FeatureCollection", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, "features":[]}
for x,y in zip(xx.ravel(), yy.ravel()): #loop setiap  x dan y dimasukkan ke dalam format geojson
    feature = {"type":"Feature","properties":{},"geometry":{"type":"MultiPolygon","coordinates":[]}}
    feature["geometry"]["coordinates"]=[[[[x-dx,y-dy],[x+dx,y-dy],[x+dx,y+dy],[x-dx,y+dy],[x-dx,y-dy]]]] #memformat koordinat untuk setiap kotak grid
    feature["properties"]={"name":"grid"}
    grid["features"].append(feature)
grid_json = json.dumps(grid)

Hasil dari proses diatas adalah sebagai berikut

Menyeleksi grid yang intersect dengan polygon input

Untuk menyeleksi polygon yang intersect dengan menggunakan library shapely. Semua variabel yang akan dicek perlu dibuat menjadi shapely shape. Cara kerja kode berikut ini yaitu dengan cara mencek setiap kotak grid apakah intersect dengan polygon input, apabila intersect maka dipertahankan, dan apabila tidak intersect maka dihapus.

s = shape(input)
grid_shape = geojson.loads(grid_json)
grid2 = {"type":"FeatureCollection", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, "features":[]}
#mengecek satu persatu polygon grid yang intersect dengan polygon input
for row in grid_shape['features']:
    r = shape(row['geometry'])
    if r.intersects(s):
        grid2["features"].append(row)

Menyimpan hasil menjadi geojson

Hasil yang masih berupa dictionary perlu diubah menjadi json dengan cara json.dumps(). Json dengan struktur berupa geojson disimpan menjadi file geojson dengan kode berikut

y = json.dumps(grid2) #membuat json dari dictionary
with open('outgrid.geojson', 'w') as f:
    f.write(str(y))
grid_json

Hasil akhir dari pembuatan grid adalah sebagai berikut:

Catatan: Apabila input ataupun outpun yang diinginkan berupa SHP atau ekstensi lainnya maka dapat menggunakan library ‘fiona’ atau library ‘ogr’

Leave a Reply