Workflow Machine Learning Pada Citra Satelit dan Data Geospasial Menggunakan Pyspatialml
What if a country has an open and easy access to geospatial data…
Pesan sponsor: seluruh konten dan kode yang saya tulis akan selalu gratis dan dapat digunakan oleh siapa saja. Apabila Anda merasa terbantu, dukung saya melalui laman Saweria berikut. 10% dari dukungan Anda akan saya sumbangkan melalui Kitabisa. Terima kasih :)
Kalau dipikir-pikir, selama saya kuliah geologi dulu, pengenalan GIS dari kurikulum bisa dibilang… “cukup”. Yang penting kami bisa membedakan WGS84 dan UTM, bisa memplot stasiun pengamatan lapangan dari gawai GPS, bisa mendapatkan peta kontur untuk lapangan pemetaan, serta memanfaatkan citra satelit untuk picking kelurusan. I only know the bare minimum…
Padahal, kemampuan analisis geospasial adalah salah satu kemampuan paling dicari di berbagai industri saat ini.
Karena itulah, ketika mas Hazred Umar Fathan dari Dunia Tambang bertanya mengenai analisis prediktif pada data geospasial, saya agak bingung dalam menjawabnya. Pada saat itu saya langsung bilang bahwa saya belum pernah melakukannya, dan saya harus mempelajari dulu…
Tantangan Evergreen
Sebagai konteks, Mas Fathan membagikan tautan Evergreen: Generate new knowledge by predicting all Australian mineral deposits kepada saya. Seperti judulnya, Evergreen adalah kontes yang ambisius, dan menuntut peserta untuk mencoba memprediksi keberadaan endapan mineral di Australia. Kontes ini diadakan oleh Unearthed, sebuah organisasi inkubator startup teknologi di bidang industri sumber daya alam yang berbasis di Australia.
Peserta diberikan data ASTER yang mencitrakan gelombang elektromagnetik dari permukaan bumi. Dari panjang gelombang elektromagnetik tersebut, dapat disimpulkan unsur yang memancarkannya, seperti misalnya oksida besi, indeks silika, dan sebagainya. Selain itu, terdapat pula data-data geofisika seperti survei magnetik, survei radiometri, serta survei gravity. Semua data tersebut berformat GeoTIFF. Tidak ada yang fancy dari data yang diberikan. Mahasiswa dan praktisi geosains di Indonesia sudah biasa berinteraksi dengan data-data tersebut.
Sebagai ground truth, terdapat 1.863 koordinat endapan mineral berharga seperti emas, nikel, dan lainnya. Menariknya, semua data di atas adalah data terbuka yang dirilis oleh Pemerintah Australia melalui Geoscience Australia (Badan Geologi-nya Australia). Semua data ini dapat diakses melalui portal satu data Pemerintah Australia. Keren!
Melihat tantangan yang diberikan, menurut saya siapapun pemenangnya, akan sulit menemukan solusi yang sempurna karena kompleksnya kondisi bawah permukaan. Menyadari hal ini, Unearthed tidak hanya menggunakan leaderboard untuk menentukan pemenang kontes, tetapi lebih mementingkan kegunaan model/aplikasi yang dibuat peserta, serta hal-hal teknis seperti bagaimana cara peserta memastikan bahwa model yang dibuat tidak overfit, dan sebagainya.
Terlepas dari itu, kontes ini menunjukkan apa yang mungkin terjadi apabila data geosains di Indonesia dibuka seluas-luasnya untuk publik!
Pyspatialml dan cara prediksi dengan data raster
Tujuan tulisan ini bukan untuk memenangkan Evergreen, melainkan untuk menunjukkan penggunaan Pyspatialml dan Python dalam melakukan prediksi menggunakan data geospasial, khususnya data raster seperti citra satelit. Perlu diketahui, saat ini Pyspatialml (ditulis oleh Steven Pawley) masih dalam tahap pengembangan dan belum tersedia di PyPI: pengguna harus menginstall melalui repo GitHub. Adanya bug dalam alur kerja harus dimaklumi.
pip install git+https://github.com/stevenpawley/Pyspatialml Pyspatialml berpusat pada objek Raster(), yang mengumpulkan semua raster-raster yang akan digunakan dalam prediksi dalam satu objek (“stack”). Objek Raster() ini kompatibel dengan scikit-learn: setelah Raster() berisi citra satelit yang ingin kita gunakan siap, tahapan selanjutnya terasa seperti workflow machine learning pada umumnya.
Dalam tantangan Evergreen, kita diberikan data raster yang memiliki coverage sebesar area observasi, dan ground truth berupa data endapan yang merupakan titik koordinat. Lalu bagaimana kita dapat membangun model machine learning? Salah satu cara menyelesaikan masalah ini (yang menjadi filosofi Pyspatialml) adalah dengan:
- Mengekstrak nilai raster pada titik-titik ground truth tersebut
- Menggunakan titik ground truth sebagai target dan hasil ekstraksi nilai raster pada koordinat ground truth sebagai prediktor, kemudian melakukan training model seperti biasa
- Menggunakan model yang sudah di-train untuk melakukan prediksi pada keseluruhan raster
Gambar di bawah ini mengilustrasikan tahapan ekstraksi nilai raster pada koordinat titik ground truth.
Sekarang saatnya memulai analisis! Berikut package yang digunakan dalam tulisan ini.
from explore_australia.stamp import Stamp, get_coverages
import os
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as pltfrom pyspatialml import Raster
from pyspatialml.preprocessing import xy_coordinates, distance_to_corners
import rasterio.plot
from pyproj import CRSfrom sklearn.model_selection import cross_validate
from sklearn.ensemble import RandomForestClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
Mengunduh data
Unearthed telah menulis package untuk memudahkan pengunduhan data raster dari National Computational Infrastructure Australia. Dalam package tersebut, juga sudah disiapkan dependency yang dibutuhkan. Untuk menghindarkan konflik, disarankan untuk membuat virtual environment Anaconda dalam menggunakan package tersebut. Panduan lengkap cara instalasi dapat dilihat di repo GitHub explore_australia.
Untuk menggunakan package tersebut, kita harus:
- Membuat clone repo
explore_australiadari GitHub - Membuat virtual environment Anaconda (atau Miniconda)
- Menjalankan
setup.pydi clone repo melalui virtual environment tersebut
Untuk memilih wilayah yang akan kita unduh raster-nya, kita harus membuat objek Stamp() yang berisi koordinat titik tengah dari kotak, dan panjang sisi dari kotak tersebut. Saya menggunakan koordinat latitude -31.125° dan longitude 121.75° sebagai titik pusat kotak, dan 100 km sebagai panjang sisi. Stamp() ini menghasilkan segiempat berdimensi 100 x 100 km di sekitar Feysville, Western Australia. Ya, WA memang dikenal sebagai salah satu daerah penghasil sumberdaya alam — baik migas maupun minerba — terbesar di Australia.
Untuk mengunduh data raster, langkah selanjutnya adalah memasukkan stamp tersebut ke dalam fungsi get_coverages(). Fungsi ini akan mengunduh file GeoTIFF ke sistem Anda, dalam direktori dengan nama yang telah dimasukkan dalam argumen name (dalam hal ini nama direktorinya adalah ‘feysville’).
Memuat data sebagai objek Raster( )
Setelah berhasil menjalankan get_coverages(), semua file raster berformat GeoTIFF akan tersimpan di sistem Anda. Untuk memuat data sebagai objek Raster() di Pyspatialml, kita hanya perlu memasukkan semua nama-nama file raster yang ingin digunakan dalam pemodelan. Nama-nama ini bisa dikumpulkan dalam bentuk list.
Dalam kasus ini, entah kenapa saya tidak bisa mendapatkan hasil yang benar apabila menggunakan proyeksi asli dari data raster. Oleh karena itu, saya melakukan reproyeksi ke WGS84 (kode EPSG:4326) dengan method to_crs(). Anda bisa melakukan reproyeksi ke Coordinate Reference System (CRS) lainnya dengan menggunakan cara ini. Perlu diketahui bahwa tahapan ini membuat saya pusing, dan ini karena pemahaman GIS saya yang minim.
Sekarang saatnya memplot semua raster untuk mengecek apakah semua raster sudah dimuat atau belum. Pyspatialml sudah dilengkapi dengan method untuk memplot stack, tetapi karena keterbatasannya, saya membuat fungsi sendiri untuk memplot sebuah stack.
Dengan ini, tahap pemuatan data raster ke dalam stack sudah beres. Tentunya kita bisa melakukan penambahan/pengurangan layer, subsetting, dan bahkan melakukan operasi matematika antar layer dalam stack. Ke depannya, kita hanya perlu berurusan dengan stack ini.
Feature Engineering geolokasi
Seperti yang sudah dijelaskan sebelumnya, pada akhirnya kita akan meng-ekstrak nilai raster pada titik-titik koordinat ground truth. Artinya, kita akan kehilangan feature posisi dari raster-raster tersebut. Padahal menurut Hukum Pertama Geografi, “everything is related to everything else, but near things are more related than distant things”.
Untungnya, Pyspatialml telah menyiapkan fungsi xy_coordinates() untuk feature engineering koordinat X-Y, dan distance_to_corners() untuk feature engineering jarak dari sudut kotak. Hasilnya akan tersimpan dalam bentuk raster, dan dapat digunakan dalam training model kita nanti. Tambahkan raster-raster baru ini ke stack utama kita dengan menggunakan method append().
Sekarang raster kita sudah siap! Saatnya menyiapkan data ground truth untuk melakukan prediksi.
Memuat data endapan mineral
Untuk memuat data endapan mineral, kita akan menggunakan geopandas. geopandas adalah dataframe mirip pandas, tapi dengan tambahan informasi sistem koordinat CRS. Dari sisi penggunaan, selain embel-embel ‘geo’ pada namanya, hampir tidak ada perbedaan antara dataframe pandas dan geopandas.
Untuk memuat data GeoJSON, kita dapat menggunakan method read_file. Saya melakukan subsetting menggunakan koordinat sudut bounding box dari raster, dan menyisakan hanya titik-titik yang berada dalam area observasi. Dari 1.863 titik di seluruh Australia, terdapat 150 titik endapan di dalam area observasi kita.
Untuk memastikan data endapan telah termuat dengan benar dan memiliki sistem koordinat yang sama dengan raster, saya melakukan pengeplotan titik di atas salah satu raster yaitu raster anomali gravity Bouguer. Dapat kita lihat bahwa titik-titik endapan kebanyakan berada di atas lokasi dengan nilai anomali tinggi.
Selain itu, dapat dilihat bahwa satu titik endapan bisa terdiri dari beberapa mineral yang disimpan dalam bentuk string dengan pemisah semikolon (;) di data asli. Banyak taktik yang dapat dilakukan tergantung kebutuhan model, tetapi saya memilih untuk melakukan one-hot encoding untuk setiap mineral. Misalnya, apabila suatu titik terdiri dari Ni;Co, maka saya akan membuat dua kolom baru untuk flag is_Ni dan is_Co.
Mengekstrak nilai raster pada koordinat titik endapan
Setelah berhasil memuat data, kita harus menyiapkan data dalam format yang bisa diproses oleh scikit-learn. Untuk itu, kita harus menyiapkan dataframe geopandas berisi flag endapan mineral dan nilai raster-raster pada koordinat endapan tersebut.
Pyspatialml punya fungsi extract_vector untuk melakukan hal ini.
<class 'geopandas.geodataframe.GeoDataFrame'>
Int64Index: 150 entries, 0 to 149
Data columns (total 33 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 is_Au 150 non-null float64
1 is_Ni 150 non-null float64
2 is_Ag 150 non-null float64
3 is_Co 150 non-null float64
4 is_Cu 150 non-null float64
5 is_PGE 150 non-null float64
6 is_Pb 150 non-null float64
7 is_Sn 150 non-null float64
8 is_Ta 150 non-null float64
9 is_Zn 150 non-null float64
10 thermal_infrared_silica_index 150 non-null float32
11 ferric_oxide_content 150 non-null float32
12 kaolin_group_index 150 non-null float32
13 ferrous_iron_content 150 non-null float32
14 aloh_groun_content 150 non-null float32
15 thermal_infrared_gypsum_index 150 non-null float32
16 tir_quartz_index 150 non-null float32
17 bouger_gravity_anomaly 150 non-null float32
18 isostatic_residual_gravity_anomaly 150 non-null float32
19 total_magnetic_intensity 150 non-null float32
20 variable_reduction_to_pole 150 non-null float32
21 filtered_uranium_ppm 150 non-null float32
22 filtered_potassium_pct 150 non-null float32
23 filtered_thorium_ppm 150 non-null float32
24 filtered_terrestrial_dose 150 non-null float32
25 x_coordinates 150 non-null float32
26 y_coordinates 150 non-null float32
27 top_left 150 non-null float32
28 top_right 150 non-null float32
29 bottom_left 150 non-null float32
30 bottom_right 150 non-null float32
31 centre_indices 150 non-null float32
32 geometry 150 non-null geometryDataframe sudah siap, saatnya melakukan workflow ML pada umumnya menggunakan scikit-learn!
Machine Learning dan CV pada data spasial
Tahap berikutnya relatif mudah apabila Anda sudah familiar dengan scikit-learn. Kita cukup membagi prediktor dan target, kemudian melakukan training pada data-data tersebut. Dalam contoh ini, saya ingin menjadikan emas sebagai target (variabel y), dan nilai-nilai raster sebagai prediktor (variabel X). Algoritma random forest digunakan dalam contoh ini. Ya, model kita sudah siap…
Namun, ada satu hal yang penting dan agak berbeda apabila kita menganalisis data spasial, yaitu di tahapan cross validation (CV). Apabila kita mengambil data train dan validation secara random seperti CV pada umumnya, bisa jadi kita mendapatkan hasil yang terlalu optimis. Padahal, kebetulan poin-poin train dan validation letaknya berdekatan (sehingga berkorelasi).
Untuk mengatasinya, kita dapat melakukan k-means clustering, lalu melakukan validation fold antar cluster. Hal ini dapat meminimalisir korelasi spasial antara poin-poin train dan validation dan membuat hasil CV lebih valid. Gambar oleh Schratz dkk (2018) mengilustrasikan masalah ini, dan bagaimana clustering dapat meminimalisir hal tersebut.
Untuk implementasinya, pertama-tama kita melakukan clustering kepada semua titik endapan berdasarkan koordinat x-y nya. Artinya, kita mengelompokkan titik-titik yang berdekatan. Kemudian, kita gunakan cluster itu sebagai grouping dalam cross_validate.
0.7999999999999999Setelah cukup yakin dengan model yang kita buat, saatnya kita melakukan prediksi pada keseluruhan raster. Objek Stack() didesain untuk workflow scikit-learn, dan kita bisa menggunakan method predict_proba untuk membuat raster berisi probabilitas kehadiran target (dalam hal ini emas) dengan nilai antara 0–1.
result_probs = stack_fin.predict_proba(estimator=clf, nodata=0)
Dari raster di atas, sepertinya lokasi memegang peranan penting dalam analisis kita. Terlihat bahwa nilai probabilitas kehadiran emas di bagian utara wilayah observasi relatif lebih tinggi. Plot feature importance juga mengkonfirmasi hipotesis tersebut. Namun, saya berekspektasi bahwa pengukuran geofisika yang lebih klasik seperti anomali gravity lebih penting dibanding hal-hal tersebut.
Tentunya model yang saya buat dalam tulisan ini perlu dievaluasi lebih mendalam. Saya tidak melakukan seleksi feature, padahal ada feature yang berkorelasi antara satu dan lainnya. Selain itu, karena saya tidak familiar dengan eksplorasi emas, saya tidak tahu citra satelit apa yang seharusnya kita perhatikan lebih mendalam. Hal ini menunjukkan bahwa domain knowledge sangat penting, terlepas dari metode yang kita gunakan.
Kesimpulan dan saran
Rupanya seru juga bermain-main dengan data geospasial! Selain itu, saya salut dengan keterbukaan data geosains Pemerintah Australia, karena dengan ini, setiap orang dapat mencoba dan berkontribusi terhadap kemajuan industri! Berikut beberapa catatan yang saya rangkum dari proses ini:
- Pemahaman mengenai sistem koordinat serta proyeksi sangat penting sebelum kita memproses data geospasial
- Pemahaman mengenai korelasi spasial sangat diperlukan dalam ML pada data geospasial
- Penambahan fitur dari shapefile/vektor data geologi (sesar, litologi, dll.) pastinya akan menjadikan model lebih robust
- Contoh dalam tulisan ini hanya untuk mendemonstrasikan workflow ML pada data raster. Evaluasi hasil ML model secara mendalam dan dengan domain knowledge diperlukan untuk menghasilkan model yang robust
Pyspatialmlsangat praktis digunakan, tetapi banyak batasan dalam penggunaannya. Apabila diperlukan kebebasan lebih, maka sepertinya workflow di luarPyspatialmlharus digunakan. Misalnya, proses raster secara manual denganRasterio- Workflow yang sama dapat diaplikasikan dalam konteks lainnya dalam geosains, misalnya pada data seismik 3D, maupun di luar geosains seperti dalam analisis ekonometrika dan/atau business intelligence
Terima kasih telah membaca tulisan ini! Seperti biasa, notebook yang saya gunakan dapat diakses dalam repo GitHub berikut. Tentunya kritik, komentar, dan pertanyaan Anda sangat ditunggu. Sampai jumpa!
