Bab 39 — GIS 2: Data Raster: Foto Udara yang Bisa Dihitung
Seri I5 — Aplikasi Pertama: Web, Peta, dan Dasbor
Semua skrip bab ini ada di folder kode/bab39/. Di versi web, Anda juga bisa mengunduh semuanya sekaligus (zip).
Bab 38 menggambar lahan dengan garis. Tetapi foto satelit atau foto drone tidak berisi garis: isinya kisi kotak-kotak kecil, masing-masing punya satu angka. Itulah raster. Bayangkan selembar kain kotak-kotak di atas kebun Anda; tiap kotak dicat sesuai warna rata-rata tanah di bawahnya. Kotak itu disebut piksel.
Tiga istilah pokok:
- Piksel: satu kotak pada kisi, dengan satu nilai per band.
- Resolusi: ukuran sisi satu piksel di lapangan. Resolusi 10 m berarti satu piksel mewakili 10 m kali 10 m (100 m²). Makin kecil angkanya, makin halus citranya, dan makin besar berkasnya.
- Band: satu “lapisan” nilai. Foto biasa punya tiga band (merah, hijau, biru). Citra satelit punya band lain yang tak terlihat mata, misalnya inframerah dekat (NIR).
Tumbuhan sehat memantulkan banyak inframerah dekat dan menyerap banyak cahaya merah. Air hampir tidak memantulkan inframerah. Perbedaan ini yang dipakai untuk mengukur “kehijauan”.
Dilarang mengunduh citra besar untuk buku ini, maka kita membuat citra satelit tiruan dengan skrip: kisi 100 kolom kali 80 baris, piksel 10 m, menutupi area 1 km × 800 m yang sama dengan Bab 38, dengan dua band (merah dan NIR) dan sedikit derau acak. Tiap petak sengaja diberi tingkat kehijauan tertentu.
01_buat_citra_tiruan.py
# Bab 39 - Membuat citra satelit tiruan (2 band: merah dan inframerah dekat)
from pathlib import Path
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize
from rasterio.transform import from_origin
tmp = Path(__file__).parent / "_tmp"
tmp.mkdir(exist_ok=True)
# Kisi 100 kolom x 80 baris, piksel 10 m, pojok kiri-atas di (700000, 9226800) UTM 48S
transform = from_origin(700000, 9226800, 10, 10)
bentuk = (80, 100)
petak = gpd.read_file("data/petak.geojson").to_crs(32748)
kolam = gpd.read_file("data/kolam.geojson").to_crs(32748)
sungai = gpd.read_file("data/sungai.geojson").to_crs(32748)
# Peta "tutupan lahan": tiap piksel diberi nilai NDVI yang diinginkan
ndvi_target = np.full(bentuk, 0.15) # tanah terbuka
nilai_petak = {"H1": 0.80, "H2": 0.78, "K1": 0.60, "K2": 0.55, "S1": 0.45, "S2": 0.70}
for _, baris in petak.iterrows():
mask = rasterize([(baris.geometry, 1)], out_shape=bentuk, transform=transform, dtype="uint8")
ndvi_target[mask == 1] = nilai_petak[baris["kode"]]
air = list(kolam.geometry) + list(sungai.buffer(8).geometry) # kolam dan badan sungai
ndvi_target[rasterize(air, out_shape=bentuk, transform=transform, dtype="uint8") == 1] = -0.20
# Dari NDVI target, susun band merah dan NIR (ditambah sedikit derau acak)
rng = np.random.default_rng(42)
ndvi_target = ndvi_target + rng.normal(0, 0.03, bentuk)
merah = rng.uniform(0.05, 0.09, bentuk)
nir = merah * (1 + ndvi_target) / (1 - ndvi_target)
citra = np.stack([merah, nir]).astype("float32")
profil = dict(driver="GTiff", height=80, width=100, count=2, dtype="float32",
crs="EPSG:32748", transform=transform)
with rasterio.open(tmp / "citra_tiruan.tif", "w", **profil) as tujuan:
tujuan.write(citra)
tujuan.set_band_description(1, "merah")
tujuan.set_band_description(2, "nir")
print("tersimpan:", (tmp / "citra_tiruan.tif").name, "| bentuk:", citra.shape)
# Bukti: hasilnya sama persis dengan berkas yang disediakan di data/
with rasterio.open("data/citra_tiruan.tif") as ada:
print("identik dengan data/citra_tiruan.tif:", np.array_equal(ada.read(), citra))tersimpan: citra_tiruan.tif | bentuk: (2, 80, 100)
identik dengan data/citra_tiruan.tif: TrueSkrip ini membuat citra di _tmp/ dan membuktikan bahwa hasilnya identik dengan berkas data/citra_tiruan.tif yang disediakan. Skrip-skrip berikutnya membaca berkas yang disediakan itu. Dua alat baru terlihat di sini: rasterize mengubah poligon menjadi piksel, dan from_origin(x, y, lebar, tinggi) menentukan letak kisi di bumi (pojok kiri-atas dan ukuran piksel).
Membaca raster dengan Rasterio
Rasterio membaca dan menulis raster (memakai GDAL di belakangnya). Isi tiap band dibaca sebagai array NumPy (Bab 23).
02_baca_raster.py
# Bab 39 - Membaca raster: piksel, resolusi, band, dan koordinat
import rasterio
with rasterio.open("data/citra_tiruan.tif") as citra:
print("ukuran (lebar x tinggi):", citra.width, "x", citra.height, "piksel")
print("jumlah band:", citra.count, "| deskripsi:", citra.descriptions)
print("tipe data:", citra.dtypes[0])
print("CRS:", citra.crs)
print("resolusi (m):", citra.res)
print("batas (kiri, bawah, kanan, atas):", tuple(citra.bounds))
merah = citra.read(1) # band 1 -> array NumPy dua dimensi
nir = citra.read(2)
print("\nbentuk array:", merah.shape, "(baris, kolom)")
print("merah: min %.3f maks %.3f" % (merah.min(), merah.max()))
print("nir : min %.3f maks %.3f" % (nir.min(), nir.max()))
# Piksel (baris 10, kolom 20) -> koordinat pusatnya di peta, dan sebaliknya
x, y = citra.xy(10, 20)
print("\npusat piksel (10, 20):", x, y)
print("koordinat itu kembali ke (baris, kolom):", citra.index(x, y))
# Nilai di satu titik, mis. di dalam kolam A (timur 700122, utara 9226285)
print("nilai [merah, nir] di titik kolam A:", [round(float(v), 3) for v in next(citra.sample([(700122, 9226285)]))])ukuran (lebar x tinggi): 100 x 80 piksel
jumlah band: 2 | deskripsi: ('merah', 'nir')
tipe data: float32
CRS: EPSG:32748
resolusi (m): (10.0, 10.0)
batas (kiri, bawah, kanan, atas): (700000.0, 9226000.0, 701000.0, 9226800.0)
bentuk array: (80, 100) (baris, kolom)
merah: min 0.050 maks 0.090
nir : min 0.030 maks 1.314
pusat piksel (10, 20): 700205.0 9226695.0
koordinat itu kembali ke (baris, kolom): (10, 20)
nilai [merah, nir] di titik kolam A: [0.057, 0.038]- Citra 100 × 80 piksel, dua band bernama
merahdannir, tipefloat32, CRS EPSG:32748 (meter), resolusi 10 m. boundsmenunjukkan letaknya: dari easting 700000 sampai 701000 dan northing 9226000 sampai 9226800.- Array berbentuk
(80, 100): baris dulu, baru kolom. Itu urutan NumPy, kebalikan dari (x, y) yang biasa kita pikirkan pada peta. xy(baris, kolom)mengubah nomor piksel menjadi koordinat pusatnya;index(x, y)membalikkannya.samplemembaca nilai di satu titik: di dalam kolam A, merah 0,057 dan NIR 0,038. NIR lebih rendah daripada merah, ciri khas air.
NDVI: mengukur kehijauan
NDVI (Normalized Difference Vegetation Index) = (NIR − Merah) ÷ (NIR + Merah). Hasilnya antara −1 dan 1: air bernilai negatif, tanah terbuka dekat 0,1 sampai 0,2, dan tumbuhan lebat mendekati 0,8 ke atas.
03_ndvi.py
# Bab 39 - NDVI: indeks kehijauan tumbuhan
from pathlib import Path
import numpy as np
import rasterio
tmp = Path(__file__).parent / "_tmp"
tmp.mkdir(exist_ok=True)
with rasterio.open("data/citra_tiruan.tif") as citra:
merah, nir = citra.read(1), citra.read(2)
profil = citra.profile
# NDVI = (NIR - Merah) / (NIR + Merah), nilainya antara -1 dan 1
ndvi = (nir - merah) / (nir + merah)
print("NDVI: min %.2f | rata-rata %.2f | maks %.2f" % (ndvi.min(), ndvi.mean(), ndvi.max()))
# Kelas sederhana (batas angka hanya contoh belajar)
kelas = np.digitize(ndvi, [0.0, 0.3, 0.6]) # 0=air, 1=terbuka, 2=sedang, 3=rapat
nama = ["air", "tanah terbuka", "vegetasi sedang", "vegetasi rapat"]
for i, n in enumerate(nama):
jumlah = int((kelas == i).sum())
print(f"{n:16s}: {jumlah:5d} piksel = {jumlah * 100 / 10_000:5.1f} ha")
# Simpan NDVI sebagai GeoTIFF satu band
profil.update(count=1, dtype="float32", nodata=-9999.0)
with rasterio.open(tmp / "ndvi.tif", "w", **profil) as tujuan:
tujuan.write(ndvi.astype("float32"), 1)
print("tersimpan:", "ndvi.tif")NDVI: min -0.29 | rata-rata 0.33 | maks 0.89
air : 216 piksel = 2.2 ha
tanah terbuka : 4964 piksel = 49.6 ha
vegetasi sedang : 777 piksel = 7.8 ha
vegetasi rapat : 2043 piksel = 20.4 ha
tersimpan: ndvi.tifRata-rata NDVI seluruh citra 0,33 dan terkecil −0,29 (air). Kelas dibuat dengan np.digitize; batas 0, 0,3, dan 0,6 hanya contoh belajar (nilai sebenarnya bergantung pada sensor dan tanaman). Misalnya vegetasi rapat 2043 piksel × 100 m² = 20,4 ha. Operasi (nir - merah) / (nir + merah) bekerja pada seluruh kisi sekaligus, tanpa perulangan: kekuatan NumPy. Hasilnya disimpan sebagai GeoTIFF (nodata ditandai −9999).
Memotong citra dengan batas lahan
Biasanya kita hanya butuh citra di dalam lahan tertentu. rasterio.mask.mask memotongnya dengan poligon.
04_potong_batas_lahan.py
# Bab 39 - Memotong citra dengan batas lahan (clip)
from pathlib import Path
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.mask import mask
tmp = Path(__file__).parent / "_tmp"
tmp.mkdir(exist_ok=True)
petak = gpd.read_file("data/petak.geojson")
with rasterio.open("data/citra_tiruan.tif") as citra:
s1 = petak[petak["kode"] == "S1"].to_crs(citra.crs) # CRS vektor = CRS raster
potongan, transform = mask(citra, s1.geometry, crop=True, nodata=np.nan)
profil = citra.profile
merah, nir = potongan
ndvi = (nir - merah) / (nir + merah)
print("bentuk potongan (band, baris, kolom):", potongan.shape)
print("piksel di dalam S1:", int(np.isfinite(ndvi).sum()))
print("NDVI rata-rata S1: %.2f" % np.nanmean(ndvi))
profil.update(height=potongan.shape[1], width=potongan.shape[2], transform=transform, nodata=np.nan)
with rasterio.open(tmp / "potongan_S1.tif", "w", **profil) as tujuan:
tujuan.write(potongan)
print("tersimpan:", "potongan_S1.tif")bentuk potongan (band, baris, kolom): (2, 17, 21)
piksel di dalam S1: 300
NDVI rata-rata S1: 0.45
tersimpan: potongan_S1.tifPotongan berukuran 17 baris kali 21 kolom: kotak yang memuat S1 (aslinya 15 baris × 20 kolom piksel), sedikit lebih besar karena batas kotak dibulatkan keluar (penyebab pastinya tidak saya periksa). Yang benar-benar di dalam petak: 300 piksel = 3 hektare, sama dengan luas S1 di Bab 38. NDVI rata-ratanya 0,45. Perhatikan .to_crs(citra.crs): poligon harus berCRS sama dengan raster. Bila lupa:
08_salah_crs_mask.py
# Bab 39 - Jebakan: batas lahan berbeda CRS dengan raster
import geopandas as gpd
import rasterio
from rasterio.mask import mask
petak = gpd.read_file("data/petak.geojson") # EPSG:4326 (derajat), lupa to_crs!
with rasterio.open("data/citra_tiruan.tif") as citra: # raster EPSG:32748 (meter)
print("CRS raster:", citra.crs, "| CRS vektor:", petak.crs)
mask(citra, petak.geometry, crop=True)CRS raster: EPSG:32748 | CRS vektor: EPSG:4326
Traceback (most recent call last):
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\mask.py", line 80, in raster_geometry_mask
window = geometry_window(dataset, shapes, pad_x=pad_x, pad_y=pad_y)
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\features.py", line 522, in geometry_window
bounding_window = bounding_window.intersection(raster_window)
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\windows.py", line 787, in intersection
return intersection([self, other])
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\windows.py", line 125, in wrapper
return function(*args[0])
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\windows.py", line 239, in intersection
return functools.reduce(_intersection, windows)
~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\windows.py", line 257, in _intersection
raise WindowError(f"Intersection is empty {w1} {w2}")
rasterio.errors.WindowError: Intersection is empty Window(col_off=-69990, row_off=922680, width=1, height=1) Window(col_off=0, row_off=0, width=100, height=80)
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "08_salah_crs_mask.py", line 9, in <module>
mask(citra, petak.geometry, crop=True)
~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\mask.py", line 180, in mask
shape_mask, transform, window = raster_geometry_mask(
~~~~~~~~~~~~~~~~~~~~^
dataset, shapes, all_touched=all_touched, invert=invert, crop=crop,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
pad=pad, pad_width=pad_width)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "<folder-proyek>\.venv\Lib\site-packages\rasterio\mask.py", line 86, in raster_geometry_mask
raise ValueError('Input shapes do not overlap raster.')
ValueError: Input shapes do not overlap raster.Rasterio menolak dengan Input shapes do not overlap raster: koordinat derajat (sekitar 106) tidak tumpang tindih dengan koordinat meter (sekitar 700000). Syukurlah di sini ada error. Bahkan tanpa error, hasil seperti ini selalu patut dicurigai.
Reproject dan resample
Reproject mengubah CRS raster; resample mengubah ukuran piksel. Keduanya memaksa komputer menghitung ulang nilai piksel, dan ada pilihan cara: bilinear (rata-rata berbobot dari tetangga, baik untuk data kontinu seperti NDVI), average, atau nearest (menyalin tetangga terdekat, wajib untuk data kelas seperti peta tutupan lahan, supaya tidak muncul “kelas 1,5”).
05_reproject_resample.py
# Bab 39 - Reproject (ganti CRS) dan resample (ganti resolusi)
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import calculate_default_transform, reproject
with rasterio.open("data/citra_tiruan.tif") as citra:
nir = citra.read(2)
# Reproject: UTM (meter) -> WGS 84 (derajat)
transform, lebar, tinggi = calculate_default_transform(
citra.crs, "EPSG:4326", citra.width, citra.height, *citra.bounds)
hasil = np.zeros((tinggi, lebar), dtype="float32")
reproject(source=nir, destination=hasil, src_transform=citra.transform, src_crs=citra.crs,
dst_transform=transform, dst_crs="EPSG:4326", resampling=Resampling.bilinear)
print("sebelum:", citra.width, "x", citra.height, "piksel, satuan meter")
print("sesudah:", lebar, "x", tinggi, "piksel, ukuran piksel (derajat): %.6f" % transform.a)
# Resample: baca langsung dengan resolusi 3x lebih kasar (30 m), nilai dirata-ratakan
kasar = citra.read(2, out_shape=(citra.height // 3, citra.width // 3), resampling=Resampling.average)
print("\nresolusi 10 m :", nir.shape, "piksel")
print("resolusi 30 m :", kasar.shape, "piksel")
print("rata-rata NIR 10 m vs 30 m: %.4f vs %.4f" % (nir.mean(), kasar.mean()))sebelum: 100 x 80 piksel, satuan meter
sesudah: 100 x 80 piksel, ukuran piksel (derajat): 0.000090
resolusi 10 m : (80, 100) piksel
resolusi 30 m : (26, 33) piksel
rata-rata NIR 10 m vs 30 m: 0.2249 vs 0.2249Setelah reproject ke EPSG:4326, ukuran piksel menjadi sekitar 0,00009 derajat (kira-kira 10 m). Resolusi 10 m menjadi 30 m membuat kisi 26 × 33 piksel (sepertiga ukuran, dibulatkan ke bawah). Rata-rata NIR sebelum dan sesudah sama (0,2249): Resampling.average menjaga nilai rata-rata.
Statistik zonal
Statistik zonal merangkum nilai raster di dalam tiap zona (poligon): rata-rata, terkecil, terbesar. Kita membuat masker untuk tiap poligon dengan geometry_mask, lalu menghitung dengan NumPy.
06_statistik_zonal.py
# Bab 39 - Statistik zonal: ringkasan nilai raster di dalam tiap poligon
import geopandas as gpd
import rasterio
from rasterio.features import geometry_mask
with rasterio.open("data/citra_tiruan.tif") as citra:
merah, nir = citra.read(1), citra.read(2)
ndvi = (nir - merah) / (nir + merah)
petak = gpd.read_file("data/petak.geojson").to_crs(citra.crs)
kolam = gpd.read_file("data/kolam.geojson").to_crs(citra.crs)
def zonal(zona):
hasil = []
for _, z in zona.iterrows():
# invert=True: True di DALAM poligon, False di luarnya
dalam = geometry_mask([z.geometry], out_shape=ndvi.shape, transform=citra.transform, invert=True)
nilai = ndvi[dalam]
hasil.append((z["kode"], int(dalam.sum()), nilai.mean(), nilai.min(), nilai.max()))
return hasil
print("zona piksel NDVI_rata NDVI_min NDVI_maks")
for kode, n, rata, kecil, besar in zonal(petak):
print(f"petak {kode:>2s} {n:4d} {rata:9.2f} {kecil:8.2f} {besar:9.2f}")
for kode, n, rata, kecil, besar in zonal(kolam):
print(f"kolam {kode:>2s} {n:4d} {rata:9.2f} {kecil:8.2f} {besar:9.2f}")zona piksel NDVI_rata NDVI_min NDVI_maks
petak H1 900 0.80 0.71 0.89
petak H2 695 0.78 0.67 0.88
petak K1 400 0.60 0.51 0.68
petak K2 300 0.55 0.45 0.63
petak S1 300 0.45 0.36 0.55
petak S2 225 0.70 0.62 0.76
kolam A 12 -0.20 -0.24 -0.15
kolam B 12 -0.21 -0.25 -0.15
kolam C 12 -0.20 -0.25 -0.14
kolam D 16 -0.19 -0.26 -0.14Hutan H1 dan H2 hijau lebat (0,80 dan 0,78), sawah S1 paling rendah (0,45; pada data tiruan ini nilai itu sengaja dibuat rendah, seolah padinya masih muda), dan keempat kolam berNDVI negatif (sekitar −0,2): air. Jumlah piksel pun cocok dengan luas: H1 900 piksel × 100 m² = 9 ha.
Proyek Bersambung I5: NDVI per petak
Hasil zonal dimasukkan kembali ke tabel petak dan ditandai.
07_proyek_ndvi_petak.py
# Bab 39 - Proyek bersambung I5: NDVI per petak, disimpan ke GeoPackage
from pathlib import Path
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import geometry_mask
tmp = Path(__file__).parent / "_tmp"
tmp.mkdir(exist_ok=True)
with rasterio.open("data/citra_tiruan.tif") as citra:
merah, nir = citra.read(1), citra.read(2)
ndvi = (nir - merah) / (nir + merah)
petak = gpd.read_file("data/petak.geojson").to_crs(citra.crs)
rata = []
for geom in petak.geometry:
# invert=True: True DI DALAM poligon
dalam = geometry_mask([geom], out_shape=ndvi.shape, transform=citra.transform, invert=True)
rata.append(round(float(ndvi[dalam].mean()), 2))
petak["ndvi_rata"] = rata
petak["status"] = np.where(petak["ndvi_rata"] < 0.5, "perlu dicek", "baik") # batas 0,5 = contoh
petak["luas_ha"] = (petak.area / 10_000).round(2)
print(petak[["kode", "sektor", "komoditas", "luas_ha", "ndvi_rata", "status"]].to_string(index=False))
petak.to_file(tmp / "petak_ndvi.gpkg", layer="petak_ndvi", driver="GPKG")
print("\ntersimpan: petak_ndvi.gpkg | kolom:", list(gpd.read_file(tmp / "petak_ndvi.gpkg").columns))kode sektor komoditas luas_ha ndvi_rata status
H1 hutan Jati 9.00 0.80 baik
H2 hutan Meranti 6.93 0.78 baik
K1 kebun kopi 4.00 0.60 baik
K2 kebun kakao 3.00 0.55 baik
S1 sawah padi 3.00 0.45 perlu dicek
S2 sawah padi 2.25 0.70 baik
tersimpan: petak_ndvi.gpkg | kolom: ['kode', 'sektor', 'komoditas', 'ndvi_rata', 'status', 'luas_ha', 'geometry']Hanya S1 (0,45) berstatus “perlu dicek” karena NDVI di bawah batas 0,5, yang di sini hanya angka contoh. Tabel itu disimpan ke GeoPackage: data vektor dan hasil raster akhirnya bersatu. Ingat bahwa citra ini tiruan; pada citra nyata, NDVI perlu dibandingkan dengan musim tanam dan pengamatan lapangan.
Latihan
- Di
03_ndvi.py, ubah batas kelas menjadi[0.0, 0.2, 0.5]dan lihat luas tiap kelas berubah. - Di
04_potong_batas_lahan.py, potong petakK1(kebun kopi) dan cetak NDVI rata-ratanya. - Di
06_statistik_zonal.py, tambahkan kolom simpangan baku (values.std()). Petak mana yang paling seragam? - Tantangan: buat kolom
jumlah_piksel_di_bawah_0_5untuk tiap petak dan tampilkan persentasenya. Bandingkan dengan status “perlu dicek”.
Ringkasan Bab
- Raster = kisi piksel; resolusi = ukuran piksel di lapangan; band = lapisan nilai.
- Rasterio membaca band sebagai array NumPy;
xydanindexberpindah antara piksel dan koordinat. - NDVI = (NIR − merah) / (NIR + merah): air negatif, tanah terbuka rendah, tumbuhan lebat tinggi.
maskmemotong dengan poligon berCRS sama;reprojectdan resample menghitung ulang piksel (pilih metode sesuai jenis data).- Statistik zonal merangkum raster per poligon; hasilnya bisa digabung ke tabel vektor.
Selanjutnya, di Bab 40: peta statis dan interaktif, PostGIS, dan mengendalikan QGIS dari Python.