Bab 38 — GIS 1: Data Vektor: Peta yang Bisa Dihitung

Seri I5 — Aplikasi Pertama: Web, Peta, dan Dasbor

CatatanDi bab ini, Anda akan:

• Memahami data spasial: titik, garis, dan poligon
• Memakai Shapely (bentuk) dan GeoPandas (tabel berbentuk)
• Membedakan format GeoJSON, GeoPackage, dan Shapefile
• Memahami CRS (sistem koordinat) dan proyeksi, serta menghindari jebakan derajat
• Menghitung luas petak, jarak kolam ke sungai, buffer, overlay, dan spatial join

Semua skrip bab ini ada di folder kode/bab38/. Di versi web, Anda juga bisa mengunduh semuanya sekaligus (zip).

Sampai sekarang, “lokasi” di buku ini hanya nama: kolam A, petak H1. Padahal pertanyaan lapangan sering menyangkut letak: Berapa hektare sawah kita? Kolam mana yang paling dekat ke sungai? Petak mana yang masuk zona penyangga? GIS (Geographic Information System, sistem informasi geografis) menjawabnya dengan menyimpan bentuk dan letak benda di bumi sebagai data, lalu menghitungnya. Bab ini membahas data vektor, yaitu benda yang digambar dengan koordinat. (Data raster, yaitu kisi piksel seperti foto satelit, ada di Bab 39.)

Seluruh data bab ini tiruan buatan skrip dan letaknya fiktif, jadi angkanya hanya untuk belajar. Isinya enam petak (hutan, kebun, sawah), empat kolam, tiga kandang, satu sungai, dan dua dusun, tersimpan di kode/bab38/data/ sebagai berkas GeoJSON.

Tiga bentuk dasar data vektor.

Titik, garis, poligon: Shapely

Shapely adalah pustaka untuk membuat dan menghitung bentuk geometri. Ia tidak tahu apa-apa tentang bumi: ia hanya melihat angka pada bidang datar.

01_geometri_shapely.py

# Bab 38 - Titik, garis, dan poligon dengan Shapely
from shapely.geometry import LineString, Point, Polygon

# Koordinat di bawah ini dalam meter (bidang datar, satuan sama di kedua sumbu)
kolam = Point(120, 285)                                  # titik: satu lokasi
sungai = LineString([(0, 330), (200, 310), (400, 290)])   # garis: urutan titik
petak = Polygon([(50, 100), (250, 100), (250, 250), (50, 250)])  # poligon: area tertutup

print(kolam.geom_type, sungai.geom_type, petak.geom_type)
print("panjang sungai (m):", round(sungai.length, 1))
print("luas petak (m2)  :", petak.area)
print("keliling petak (m):", petak.length)
print("titik tengah petak:", petak.centroid)

print("kolam di dalam petak?  ", petak.contains(kolam))
print("sungai memotong petak? ", sungai.intersects(petak))
print("jarak kolam ke sungai (m):", round(kolam.distance(sungai), 1))
print("batas luar (minx, miny, maxx, maxy):", petak.bounds)

Point LineString Polygon
panjang sungai (m): 402.0
luas petak (m2)  : 30000.0
keliling petak (m): 700.0
titik tengah petak: POINT (150 175)
kolam di dalam petak?   False
sungai memotong petak?  False
jarak kolam ke sungai (m): 32.8
batas luar (minx, miny, maxx, maxy): (50.0, 100.0, 250.0, 250.0)

• Point untuk satu lokasi (kolam, kandang, pohon), LineString untuk urutan titik (sungai, jalan), Polygon untuk area tertutup (petak).
• Luas petak 30000 m² (200 m kali 150 m) dan kelilingnya 700 m. Sungai (dari tiga titik) panjangnya 402 m.
• contains, intersects, distance menjawab pertanyaan letak: kolam tidak di dalam petak, sungai tidak memotong petak, dan jarak kolam ke sungai 32,8 m.

CatatanINGAT

Shapely menghitung dalam satuan angka koordinat. Bila koordinat dalam meter, luasnya meter persegi. Bila dalam derajat, “luas”-nya derajat persegi, yang tidak berguna. Itu persoalan CRS, dibahas sebentar lagi.

GeoPandas: tabel yang punya bentuk

GeoPandas menambahkan satu kolom khusus, geometry, ke DataFrame pandas (Bab 23). Semua yang Anda tahu tentang pandas (menyaring, groupby, menggabung) tetap berlaku. Mesin yang membaca dan menulis berkas peta di belakang layar adalah pyogrio (pengganti fiona, yang belum tersedia untuk Python 3.14 saat bab ini ditulis).

02_geodataframe_format.py

# Bab 38 - GeoDataFrame dan tiga format berkas: GeoJSON, GeoPackage, Shapefile
import warnings
from pathlib import Path

import geopandas as gpd
from shapely.geometry import Point

tmp = Path(__file__).parent / "_tmp"
tmp.mkdir(exist_ok=True)

# GeoDataFrame = tabel pandas + satu kolom bernama geometry + CRS
kandang = gpd.GeoDataFrame(
    {"kode": ["D1", "D2"], "jenis": ["sapi", "ayam"], "populasi_ekor_hari_ini": [12, 200]},
    geometry=[Point(106.8168, -6.9955), Point(106.8118, -6.9945)],   # (bujur, lintang)
    crs="EPSG:4326",
)
print(kandang)
print("CRS:", kandang.crs)

kandang.to_file(tmp / "kandang.geojson", driver="GeoJSON")
kandang.to_file(tmp / "kandang.gpkg", layer="kandang", driver="GPKG")
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    kandang.to_file(tmp / "kandang.shp")                          # Shapefile
    print("peringatan Shapefile:", str(w[0].message))

print("berkas Shapefile:", sorted(p.name for p in tmp.glob("kandang.*") if p.suffix not in (".geojson", ".gpkg")))
print("kolom dari Shapefile  :", list(gpd.read_file(tmp / "kandang.shp").columns))
print("kolom dari GeoPackage :", list(gpd.read_file(tmp / "kandang.gpkg").columns))

  kode jenis  populasi_ekor_hari_ini                  geometry
0   D1  sapi                      12  POINT (106.8168 -6.9955)
1   D2  ayam                     200  POINT (106.8118 -6.9945)
CRS: EPSG:4326
peringatan Shapefile: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
berkas Shapefile: ['kandang.cpg', 'kandang.dbf', 'kandang.prj', 'kandang.shp', 'kandang.shx']
kolom dari Shapefile  : ['kode', 'jenis', 'populasi_e', 'geometry']
kolom dari GeoPackage : ['kode', 'jenis', 'populasi_ekor_hari_ini', 'geometry']

Tiga format yang paling sering Anda temui:

Format	Berkas	Kelebihan	Kekurangan
GeoJSON	satu .geojson (teks)	Mudah dibaca manusia dan web; selalu EPSG:4326	Boros ruang; lambat untuk data besar
GeoPackage	satu .gpkg (SQLite)	Satu berkas, banyak lapisan, nama kolom bebas, standar terbuka	Tidak terbaca manusia
Shapefile	lima berkas atau lebih (.shp, .shx, .dbf, .prj, .cpg)	Dibaca hampir semua perangkat lunak lama	Nama kolom maksimum 10 huruf

Keluarannya menunjukkan batasan Shapefile: kolom populasi_ekor_hari_ini dipotong menjadi populasi_e (dan Anda diberi peringatan), sedangkan GeoPackage menyimpannya utuh. Untuk data baru, pilih GeoPackage.

CRS: derajat atau meter?

Bumi bulat, peta datar. CRS (Coordinate Reference System) menjelaskan arti angka koordinat. Dua yang paling penting:

• EPSG:4326 (WGS 84): koordinat bujur dan lintang dalam derajat. Dipakai GPS dan web.
• Proyeksi UTM, misalnya EPSG:32748 (UTM zona 48 Selatan): koordinat dalam meter di bidang datar. Cocok untuk luas dan jarak di satu wilayah. Indonesia terbagi dalam beberapa zona UTM.

03_crs_proyeksi.py

# Bab 38 - Sistem koordinat (CRS): derajat atau meter?
import geopandas as gpd
from pyproj import CRS, Transformer

petak = gpd.read_file("data/petak.geojson")
print("CRS berkas:", petak.crs)
print("satuan sumbu:", [s.unit_name for s in petak.crs.axis_info])
print("pojok petak H1 (bujur, lintang):", [round(v, 5) for v in petak.geometry.iloc[0].exterior.coords[0]])
print("'luas' H1 dalam derajat persegi:", petak.geometry.iloc[0].area)   # tidak ada artinya!

utm = petak.estimate_utm_crs()                        # menebak zona UTM dari letak data
print("zona UTM yang cocok:", utm.to_string(), "|", utm.name)
datar = petak.to_crs(utm)
print("satuan sumbu UTM:", [s.unit_name for s in datar.crs.axis_info])
print("luas H1 (m2):", round(datar.geometry.iloc[0].area))

# Pyproj: pindah satu titik dari 4326 ke UTM. always_xy=True = urutan (bujur, lintang).
pindah = Transformer.from_crs(CRS.from_epsg(4326), CRS.from_epsg(32748), always_xy=True)
timur, utara = pindah.transform(106.8125, -6.9950)
print("titik di UTM 48S: timur =", round(timur, 1), "| utara =", round(utara, 1))

CRS berkas: EPSG:4326
satuan sumbu: ['degree', 'degree']
pojok petak H1 (bujur, lintang): [106.81361, -6.99516]
'luas' H1 dalam derajat persegi: 7.363732839988365e-06
zona UTM yang cocok: EPSG:32748 | WGS 84 / UTM zone 48S
satuan sumbu UTM: ['metre', 'metre']
luas H1 (m2): 90001
titik di UTM 48S: timur = 700227.2 | utara = 9226418.1

• File kita berasal dari EPSG:4326, sumbunya derajat. “Luas” H1 di situ 7.36e-06 derajat persegi: angka tanpa arti praktis.
• estimate_utm_crs() menebak zona UTM yang cocok dari letak data: EPSG:32748.
• Setelah to_crs(...), sumbu berubah menjadi meter dan luas H1 menjadi 90001 m² (kira-kira 9 hektare; selisih 1 m² berasal dari pembulatan koordinat GeoJSON ke 7 desimal, sekitar 1 cm).
• pyproj bisa memindahkan satu titik sendiri. always_xy=True memastikan urutannya (bujur, lintang), karena EPSG:4326 resminya berurutan (lintang, bujur) dan sering membingungkan.

PeringatanAWAS

Jangan menghitung luas atau jarak langsung dari data berderajat. Ubah dulu ke CRS yang berunit meter dengan to_crs. Ini kesalahan paling sering di GIS. Setiap berkas membawa CRS-nya sendiri (gdf.crs): periksa sebelum menggabung.

Luas petak dan jarak kolam ke sungai

04_luas_jarak.py

# Bab 38 - Luas petak dan jarak kolam ke sungai
import geopandas as gpd

petak = gpd.read_file("data/petak.geojson").to_crs(32748)    # meter
kolam = gpd.read_file("data/kolam.geojson").to_crs(32748)
sungai = gpd.read_file("data/sungai.geojson").to_crs(32748)

petak["luas_ha"] = (petak.area / 10_000).round(2)
print(petak[["kode", "sektor", "komoditas", "luas_ha"]].to_string(index=False))
print()
print("luas per sektor (ha):")
print(petak.groupby("sektor")["luas_ha"].sum().to_string())

garis_sungai = sungai.geometry.iloc[0]
kolam["jarak_sungai_m"] = kolam.distance(garis_sungai).round(1)
print()
print(kolam[["kode", "jenis", "jumlah_ikan", "jarak_sungai_m"]].to_string(index=False))
dekat = kolam[kolam["jarak_sungai_m"] < 25]
print("kolam kurang dari 25 m dari sungai:", list(dekat["kode"]))

kode sektor komoditas  luas_ha
  H1  hutan      Jati     9.00
  H2  hutan   Meranti     6.93
  K1  kebun      kopi     4.00
  K2  kebun     kakao     3.00
  S1  sawah      padi     3.00
  S2  sawah      padi     2.25

luas per sektor (ha):
sektor
hutan    15.93
kebun     7.00
sawah     5.25

kode jenis  jumlah_ikan  jarak_sungai_m
   A  lele          800            15.4
   B  nila          500            29.8
   C  lele          650            20.0
   D  nila          300           116.4
kolam kurang dari 25 m dari sungai: ['A', 'C']

Hutan 15,93 ha (9,00 + 6,93), kebun 7,00 ha, sawah 5,25 ha. kolam.distance(garis) menghitung jarak terdekat dari tiap kolam ke sungai: A dan C kurang dari 25 m, kolam D jauh (116,4 m). Pada contoh nyata, kolam dekat sungai berguna untuk air, tetapi berisiko banjir; angka batas 25 m di sini hanya contoh.

Buffer dan overlay

Buffer membuat area di sekitar suatu bentuk dengan jarak tertentu. Overlay menumpuk dua lapisan dan menghitung potongannya: intersection (bagian yang beririsan) atau difference (bagian yang tidak).

05_buffer_overlay.py

# Bab 38 - Buffer (zona penyangga) dan overlay (tumpang susun)
import geopandas as gpd

petak = gpd.read_file("data/petak.geojson").to_crs(32748)
sungai = gpd.read_file("data/sungai.geojson").to_crs(32748)

# Buffer: area dalam jarak 60 m dari sungai (angka 60 hanya contoh belajar)
zona = gpd.GeoDataFrame({"zona": ["penyangga"]}, geometry=sungai.buffer(60), crs=sungai.crs)
print("luas zona penyangga (ha):", round(zona.area.iloc[0] / 10_000, 2))

# Overlay intersection: potongan petak yang masuk zona
masuk = gpd.overlay(petak, zona, how="intersection")
masuk["luas_masuk_m2"] = masuk.area.round(0)
print("\npetak yang masuk zona penyangga:")
print(masuk[["kode", "sektor", "luas_masuk_m2"]].to_string(index=False))

# Overlay difference: bagian petak yang di luar zona
aman = gpd.overlay(petak, zona, how="difference")
aman["luas_aman_ha"] = (aman.area / 10_000).round(2)
print("\nluas petak di luar zona (ha):")
print(aman[["kode", "luas_aman_ha"]].to_string(index=False))

luas zona penyangga (ha): 13.24

petak yang masuk zona penyangga:
kode sektor  luas_masuk_m2
  K2  kebun            4.0
  S1  sawah          141.0
  S2  sawah         2333.0

luas petak di luar zona (ha):
kode  luas_aman_ha
  H1          9.00
  H2          6.93
  K1          4.00
  K2          3.00
  S1          2.99
  S2          2.02

Zona 60 m di kiri-kanan sungai luasnya 13,24 ha. Tiga petak menjorok masuk: S2 paling banyak (2333 m²), S1 hanya 141 m², dan K2 cuma 4 m². Baris terakhir menunjukkan sisanya: S2 tinggal 2,02 ha (dari 2,25) dan S1 2,99 ha. Angka 60 m hanyalah contoh; aturan sempadan sungai yang sebenarnya ditentukan peraturan setempat.

Peta lokasi usaha (data tiruan): petak, kolam, kandang, sungai, dan zona 60 m.

Spatial join: menggabung menurut letak

Di pandas, merge menggabung dua tabel menurut nilai kolom yang sama. Di GeoPandas, spatial join (sjoin) menggabung menurut hubungan letak: “titik ini di dalam poligon yang mana?”.

06_spatial_join.py

# Bab 38 - Spatial join: menggabungkan tabel berdasarkan letak
import warnings

import geopandas as gpd

dusun = gpd.read_file("data/dusun.geojson").to_crs(32748)
kolam = gpd.read_file("data/kolam.geojson").to_crs(32748)
kandang = gpd.read_file("data/kandang.geojson").to_crs(32748)

# Titik pusat tiap kolam jatuh di dusun mana? (predicate="within": di dalam)
kolam_titik = kolam.copy()
kolam_titik["geometry"] = kolam.centroid
gabung = gpd.sjoin(kolam_titik, dusun[["nama", "geometry"]], predicate="within")
print(gabung[["kode", "jumlah_ikan", "nama"]].to_string(index=False))
print("\nikan per dusun:")
print(gabung.groupby("nama")["jumlah_ikan"].sum().to_string())

# Kandang terdekat dari tiap kolam (sjoin_nearest) beserta jaraknya
dekat = gpd.sjoin_nearest(kolam_titik, kandang[["kode", "geometry"]], distance_col="jarak_m",
                          lsuffix="kolam", rsuffix="kandang")
print("\nkandang terdekat dari tiap kolam:")
print(dekat[["kode_kolam", "kode_kandang", "jarak_m"]].round(0).to_string(index=False))

# Jebakan: CRS berbeda -> hanya peringatan, hasilnya KOSONG
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    salah = gpd.sjoin(gpd.read_file("data/kolam.geojson"), dusun)    # kiri 4326, kanan UTM
print("\nCRS beda ->", str(w[0].message).splitlines()[0])
print("jumlah baris hasil:", len(salah))

kode  jumlah_ikan  nama
   A          800  Hulu
   B          500  Hulu
   C          650 Hilir
   D          300 Hilir

ikan per dusun:
nama
Hilir     950
Hulu     1300

kandang terdekat dari tiap kolam:
kode_kolam kode_kandang  jarak_m
         A           D2    466.0
         B           D2    442.0
         C           D1    403.0
         D           D1    534.0

CRS beda -> CRS mismatch between the CRS of left geometries and the CRS of right geometries.
jumlah baris hasil: 0

• Kolam A dan B jatuh di dusun Hulu, C dan D di Hilir, sehingga ikan per dusun: Hulu 1300, Hilir 950.
• sjoin_nearest mencari yang terdekat; kandang D2 terdekat dari kolam A dan B, sedangkan D1 dari C dan D.
• Jebakan di akhir: menggabung data berCRS beda tidak menimbulkan error, hanya peringatan, dan hasilnya kosong (0 baris). Selalu to_crs dulu, dan baca peringatan.

Memasang GIS: pip atau conda-forge?

Di komputer penulis, semua pustaka GIS terpasang dengan pip karena tersedia sebagai wheel (paket jadi) yang sudah membawa pustaka C di dalamnya. Tidak selalu begitu. Pustaka dasarnya, GDAL, GEOS, dan PROJ, ditulis dalam C/C++, dan di beberapa sistem pip tidak menemukan wheel yang cocok sehingga pemasangan gagal. conda-forge (saluran paket komunitas untuk conda atau mamba) memasang pustaka C itu bersama paket Python-nya sehingga biasanya lebih mulus:

	pip	conda-forge
Perintah	pip install geopandas rasterio	conda install -c conda-forge geopandas rasterio
Pustaka C (GDAL, PROJ)	Ikut di wheel, bila ada untuk sistem Anda	Dipasang sebagai paket sendiri, versinya serasi
Kelebihan	Ringan, satu cara dengan paket lain	Paling andal untuk GDAL, banyak alat GIS lain
Kekurangan	Bisa gagal bila tidak ada wheel	Perlu memasang conda/mamba; lingkungan terpisah dari venv

BahayaBAGIAN TEKNIS (boleh dilewati)

Perintah conda di atas tidak saya jalankan: conda belum terpasang di komputer penulis, jadi bagian conda-forge belum diuji di sini. Perbandingan ini berdasar cara kerja umum kedua alat. Bila pip install geopandas berhasil di komputer Anda, Anda tidak perlu conda.

Proyek Bersambung I5: peta lokasi usaha

Kita satukan semuanya: membaca lima lapisan, mengubahnya ke UTM, menyimpannya sebagai satu GeoPackage berisi lima lapisan, merangkum luas per sektor, dan membuat peta cepat dengan .plot(). GeoPackage ini bahan untuk Bab 40 (PostGIS dan QGIS).

07_proyek_peta_lokasi.py

# Bab 38 - Proyek bersambung I5: peta lokasi usaha, disimpan ke GeoPackage
from pathlib import Path

import geopandas as gpd
import matplotlib.pyplot as plt

tmp = Path(__file__).parent / "_tmp"
tmp.mkdir(exist_ok=True)
lapisan = ["petak", "kolam", "kandang", "sungai", "dusun"]
data = {nama: gpd.read_file(f"data/{nama}.geojson").to_crs(32748) for nama in lapisan}

# Satu berkas GeoPackage berisi banyak lapisan (layer)
gpkg = tmp / "lokasi_usaha.gpkg"
gpkg.unlink(missing_ok=True)
for nama, gdf in data.items():
    gdf.to_file(gpkg, layer=nama, driver="GPKG")
print("lapisan di dalam GeoPackage:", gpd.list_layers(gpkg)["name"].tolist())

# Ringkasan per sektor: luas petak dan jumlah titik
petak = data["petak"]
print("\nluas per sektor (ha):")
luas_ha = petak.assign(luas_ha=petak.area / 10_000)
print(luas_ha.groupby("sektor")["luas_ha"].sum().round(2).to_string())
print("jumlah kolam:", len(data["kolam"]), "| jumlah kandang:", len(data["kandang"]))

# Peta cepat dengan .plot()
fig, ax = plt.subplots(figsize=(7, 5.6))
petak.plot(ax=ax, column="sektor", legend=True, alpha=0.6, edgecolor="black")
data["sungai"].plot(ax=ax, color="tab:blue", linewidth=2)
data["kolam"].plot(ax=ax, color="navy")
data["kandang"].plot(ax=ax, color="red", markersize=40)
ax.set_title("Peta lokasi usaha (data tiruan)")
fig.savefig(tmp / "peta_cepat.png", dpi=100)
print("peta tersimpan:", (tmp / "peta_cepat.png").name)

lapisan di dalam GeoPackage: ['petak', 'kolam', 'kandang', 'sungai', 'dusun']

luas per sektor (ha):
sektor
hutan    15.93
kebun     7.00
sawah     5.25
jumlah kolam: 4 | jumlah kandang: 3
peta tersimpan: peta_cepat.png

gpd.list_layers membuktikan satu berkas memuat lima lapisan. Peta cepatnya tersimpan di _tmp/peta_cepat.png (peta yang lebih rapi kita buat di Bab 40).

Latihan

1. Di 01_geometri_shapely.py, buat poligon segitiga dengan titik (0,0), (100,0), (0,60) dan cetak luasnya. (Cocokkan dengan rumus setengah alas kali tinggi.)
2. Di 04_luas_jarak.py, tambahkan kolom jarak_sungai_m untuk kandang dan cetak kandang yang paling jauh dari sungai.
3. Di 05_buffer_overlay.py, ganti 60 m menjadi 100 m. Petak mana saja yang masuk, dan berapa luas yang masuk di S2?
4. Tantangan: hitung, untuk setiap dusun, luas total petak yang berada di dalamnya memakai overlay(petak, dusun, how="intersection") dan groupby. Petak H2 melintasi batas dusun: periksa bahwa jumlah luasnya tetap 6,93 ha.

Ringkasan Bab

• Data vektor = titik, garis, poligon. Shapely menghitung bentuknya; GeoPandas menaruhnya dalam tabel.
• Pilih GeoPackage untuk data baru; Shapefile memotong nama kolom dan butuh banyak berkas.
• CRS menentukan arti koordinat. Untuk luas dan jarak, ubah ke CRS bermeter (mis. UTM) dengan to_crs.
• area, distance, buffer, overlay, sjoin, dan sjoin_nearest menjawab pertanyaan lapangan tentang letak.
• CRS berbeda pada sjoin atau overlay hanya memberi peringatan dan hasil keliru: waspada.
• pip cukup bila ada wheel; conda-forge lebih andal untuk GDAL (belum diuji di sini).

Selanjutnya, di Bab 39: data raster: foto satelit tiruan, piksel, dan NDVI.