BAB 12: Kotak Perkakas Tambahan: Pustaka Ilmiah

Studi kasus: "Rata-rata dari sembilan puluh ribu piksel" #

Kepala Seksi meminta rata-rata nilai sebuah citra. Citra itu 300 kali 300 piksel, jadi ada 90.000 angka. Dengan perulangan for, Anda harus menulis dan menunggu. Ada pustaka yang menghitung semuanya sekaligus dalam satu baris. Bab ini memperkenalkan empat pustaka yang paling sering muncul di skrip GIS.

Konsep: pustaka dalam tiga kalimat #

Sebuah pustaka adalah kotak perkakas tambahan yang ditulis orang lain dan boleh Anda pakai. Python sendiri sudah membawa perkakas dasar, dan pustaka menambahkan perkakas khusus: hitungan angka massal, pembaca raster, pembuat grafik. Anda memanggilnya dengan import, seperti di Bab 8.

Istilah baru di bagian ini:

• Pustaka (library): kotak perkakas tambahan.
• Array: tabel angka berbentuk kotak, dengan baris dan kolom, dikelola oleh numpy.
• Bentuk (shape): ukuran sebuah array, misalnya 3 baris dan 4 kolom.

numpy: kalkulator untuk deretan angka #

import numpy as np

luas = np.array([4.0, 4.0, 8.0, 2.5])
print(luas.sum(), luas.mean())
print(luas[luas > 3])
print(luas * 2)
a = np.arange(12).reshape(3, 4)
print(a.shape, a[1, 2], a.sum(axis=0))

Penjelasannya:

• import numpy as np memanggil numpy dan memberinya nama pendek np. Itu kebiasaan umum.
• np.array([...]) membuat array dari sebuah daftar.
• luas.sum() dan luas.mean() menjumlahkan dan merata-ratakan seluruh isi sekaligus, tanpa for.
• luas[luas > 3] menyaring: hanya nilai yang lebih dari 3.
• luas * 2 mengalikan setiap isi dengan 2. Daftar biasa tidak bisa melakukannya.
• np.arange(12).reshape(3, 4) membuat angka 0 sampai 11, lalu menatanya menjadi 3 baris dan 4 kolom.
• a.shape memberi bentuknya. a[1, 2] mengambil baris nomor 1, kolom nomor 2. a.sum(axis=0) menjumlahkan tiap kolom.

18.5 4.625
[4. 4. 8.]
[ 8.  8. 16.  5.]
(3, 4) 6 [12 15 18 21]

osgeo (GDAL): pembaca raster #

GDAL adalah pembaca dan penulis hampir semua jenis berkas peta. Di Python, ia diimpor dari paket osgeo. Raster dibaca menjadi array numpy, jadi dua pustaka ini hampir selalu berpasangan.

Skrip yang dibedah: membaca citra di i2_04.

Tujuan. Membuka citra empat band, mengubahnya menjadi array, lalu menyusunnya menjadi tabel satu baris per piksel.

Bahan. Paket data paket-i2, berkas Citra_Lahan.tif. Jalankan di Python bawaan QGIS.

Potongan skrip (baris 15 sampai 19 di berkas aslinya):

ds = gdal.Open(os.path.join(PAKET_I2, "Citra_Lahan.tif"))
gt, proj = ds.GetGeoTransform(), ds.GetProjection()
citra = ds.ReadAsArray().astype("float64")                  # bentuk (4, tinggi, lebar)
_, tinggi, lebar = citra.shape
piksel = citra.reshape(4, -1).T                             # satu baris = satu piksel, empat kolom = empat band

Penjelasan baris demi baris. Nomor baris dihitung dari potongan di atas.

Hasil. Untuk melihat bentuknya, tambahkan beberapa print di ujung potongan. Berikut hasil di komputer penulis dengan PAKET_I2 menunjuk folder paket-i2:

4 (4, 300, 300) (90000, 4)
(312000.0, 1.0, 0.0, 9996300.0, 0.0, -1.0)
60.1

• Baris pertama: citra punya 4 band, bentuk array (4, 300, 300), dan tabel piksel (90000, 4).
• Baris kedua: gt berisi koordinat sudut kiri-atas 312000.0 dan 9996300.0, serta ukuran piksel 1 meter (1.0 dan -1.0).
• Baris ketiga: rata-rata band pertama, dibulatkan dua desimal.

Layar Anda juga bisa menampilkan peringatan FutureWarning dari GDAL tentang gdal.UseExceptions(). Itu hanya peringatan versi mendatang. Skrip tetap berjalan.

Bagian yang boleh diubah.

• Lokasi citra (baris 1). Ganti PAKET_I2 dan nama berkas dengan citra Anda.
• Angka 4 di reshape(4, -1) (baris 5). Samakan dengan jumlah band citra Anda, yaitu ds.RasterCount.
• float64 (baris 3). Bila citra Anda besar, float32 menghemat memori. [CEK]

Jebakan.

• Salah alamat. Bila alamat citra salah, gdal.Open mengembalikan None, dan baris berikutnya gagal dengan AttributeError karena None tidak punya GetGeoTransform. Periksa alamat lebih dulu.
• Jumlah band tidak cocok. reshape(4, -1) gagal bila band bukan empat. Pakai ds.RasterCount.
• Memori. Citra besar menjadi array besar di memori. Citra 300 kali 300 empat band masih kecil, tetapi citra satelit penuh bisa memenuhi memori.
• Indeks band dari 1 di GDAL, dari 0 di numpy. ds.GetRasterBand(1) adalah band pertama, tetapi citra[0] juga band pertama.

matplotlib: pembuat grafik #

import os
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

plt.bar(["Alpha", "Beta", "Gamma"], [4.0, 4.0, 8.0])
plt.ylabel("Luas (ha)")
plt.savefig("luas.png", dpi=100)
print(os.path.exists("luas.png"))

Penjelasannya:

• matplotlib.use("Agg") memilih mode tanpa jendela: grafik hanya disimpan ke berkas. Aman dijalankan di QGIS dan di terminal.
• plt.bar(nama, nilai) menggambar batang. plt.ylabel memberi judul sumbu tegak.
• plt.savefig("luas.png", dpi=100) menyimpan gambar. dpi mengatur ketajaman.

True

Baris terakhir membuktikan berkas luas.png sudah ada. Buka berkas itu dengan penampil gambar untuk melihat grafiknya.

scipy dan scikit-learn: statistik dan pembelajaran mesin #

scipy berisi banyak perkakas ilmiah. Salah satunya curve_fit, yang mencari garis terbaik untuk sekumpulan titik. Buku M2 memakainya pada skrip interpolasi.

import numpy as np
from scipy.optimize import curve_fit

x = np.array([1.0, 2.0, 3.0, 4.0])
y = np.array([2.1, 3.9, 6.2, 7.8])


def garis(x, a, b):
    return a * x + b


(a, b), _ = curve_fit(garis, x, y)
print(round(float(a), 2), round(float(b), 2))

Penjelasannya: fungsi garis adalah rumus garis lurus dengan dua angka yang belum diketahui, a dan b. curve_fit mencari a dan b yang paling cocok dengan empat titik data. Hasilnya kemiringan dan titik potong.

1.94 0.15

scikit-learn (nama modulnya sklearn) berisi alat pembelajaran mesin, misalnya klasifikasi hutan acak. Pustaka itu tidak ada di Python bawaan QGIS 4.0.2. Karena itu Buku M3 menyertakan hutan_mini.py, versi ringkas yang hanya memakai numpy, sehingga skrip tetap jalan tanpa scikit-learn.

Bagaimana mengenali pustaka di skrip orang lain #

Lihat baris import di bagian atas skrip. Bacalah:

• import numpy as np: ada hitungan angka massal. Cari np..
• from osgeo import gdal: ada pembacaan raster. Cari gdal..
• import matplotlib...: ada grafik.
• from scipy... import ...: ada statistik atau optimisasi.
• import psycopg2: ada koneksi ke basis data PostGIS (Buku M2).

Anda tidak perlu menguasai semuanya. Pada skrip yang Anda pelajari di Bagian B, perkakas ini akan muncul satu per satu, dengan penjelasan seperlunya.

Latihan variasi. Dari array [4.0, 4.0, 8.0, 2.5], hitung jumlah luas petak yang lebih dari 3 hektare memakai numpy, tanpa for.

Contoh jawaban:

import numpy as np

luas = np.array([4.0, 4.0, 8.0, 2.5])
print(luas[luas > 3].sum())

16.0

Cek paham #

1. Mengapa luas * 2 bekerja pada array numpy, tetapi tidak pada daftar biasa?
2. Pustaka mana yang membaca berkas raster?
3. Apa fungsi matplotlib.use("Agg")?

Jawaban.

1. numpy memperlakukan array sebagai satu kesatuan, sehingga operasi berlaku untuk setiap isinya. Daftar biasa hanya menggandakan daftarnya.
2. GDAL, yang diimpor dari paket osgeo.
3. Memilih mode tanpa jendela, sehingga grafik hanya disimpan ke berkas.

Kesalahan umum #

• Memakai for untuk hitungan massal. Perbaikan: pakai numpy bila datanya besar.
• Salah menjalankan. Skrip yang memakai osgeo gagal dengan ModuleNotFoundError di Python biasa. Perbaikan: jalankan dengan Python bawaan QGIS.
• Mencampur indeks band. Perbaikan: ingat GDAL mulai dari 1, numpy mulai dari 0.

Ringkasan dan latihan #

• Pustaka adalah kotak perkakas tambahan, dipanggil dengan import.
• numpy untuk hitungan massal, GDAL (osgeo) untuk raster, matplotlib untuk grafik, scipy dan sklearn untuk statistik dan pembelajaran mesin.
• Raster dibaca GDAL menjadi array numpy, dihitung, lalu ditulis kembali.

Latihan. Dengan numpy, buat array berisi lima nilai suhu, lalu cetak rata-rata, nilai terbesar, dan jumlah nilai di atas rata-rata.

Tabel perbandingan: empat pustaka #