Dokumen ini berisi kode sel demi sel (cell-by-cell) yang siap Anda salin ke Google Colab untuk melakukan ekstraksi parameter Sentinel-2, pemfilteran non-lahan, pelabelan otomatis berbasis aturan ahli, dan pelatihan model Random Forest.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import reproject
import json, os, shutil, glob, warnings
warnings.filterwarnings('ignore')
LOCAL_FOLDER = '/content/konawe_data'
GEE_FOLDER = '/content/drive/MyDrive/GEE_Konawe'
os.makedirs(LOCAL_FOLDER, exist_ok=True)
print('Library siap!')from google.colab import drive
drive.mount('/content/drive')
os.makedirs(GEE_FOLDER, exist_ok=True)
print('Drive terhubung!')import geopandas as gpd
# Muat file batas wilayah Konawe
boundary_path = f'{GEE_FOLDER}/output_FINAL/konawe_boundary.geojson'
if not os.path.exists(boundary_path):
# Jika di subfolder output_FINAL belum ada, cari di folder utama GEE_Konawe
boundary_path = f'{GEE_FOLDER}/konawe_boundary.geojson'
gdf = gpd.read_file(boundary_path)
geom = gdf.geometry.union_all()
print("Batas wilayah Konawe berhasil dimuat!")Sel ini akan men-generate titik acak di daratan Konawe secara berulang (looping) hingga mendapatkan tepat 10.000 titik yang benar-benar merupakan lahan bervegetasi aktif (lolos filter NDVI > 0.25 DAN NDWI < 0) agar terhindar dari pemukiman, jalan, dan air.
from shapely.geometry import Point
# Path dataset Sentinel-2 di Drive Anda
sentinel_path = f'{GEE_FOLDER}/output_FINAL/Sentinel2_Konawe_2025.tif'
if not os.path.exists(sentinel_path):
sentinel_path = f'{GEE_FOLDER}/Sentinel2_Konawe_2025.tif'
# Dapatkan bounding box Konawe
minx, miny, maxx, maxy = geom.bounds
np.random.seed(42)
valid_points = []
attempts = 0
max_attempts = 100000
print("Men-generate 10.000 titik sampel acak bebas pemukiman & air...")
with rasterio.open(sentinel_path) as src:
while len(valid_points) < 10000 and attempts < max_attempts:
# Generate koordinat acak di dalam bbox
lon = np.random.uniform(minx, maxx)
lat = np.random.uniform(miny, maxy)
p = Point(lon, lat)
# 1. Pastikan titik berada di dalam daratan Konawe
if geom.contains(p):
# 2. Ekstrak nilai Sentinel-2 untuk cek NDVI & NDWI
coords = [(lon, lat)]
val = next(src.sample(coords))
b2, b3, b4, b8 = float(val[0]), float(val[1]), float(val[2]), float(val[3])
epsilon = 1e-8
ndvi = (b8 - b4) / (b8 + b4 + epsilon)
ndwi = (b3 - b8) / (b3 + b8 + epsilon)
# 3. Filter: Harus lahan bervegetasi (NDVI > 0.25) & bukan badan air (NDWI < 0)
if ndvi > 0.25 and ndwi < 0:
valid_points.append({
'longitude': lon,
'latitude': lat,
'B2_Blue': b2,
'B3_Green': b3,
'B4_Red': b4,
'B8_NIR': b8,
'NDVI': ndvi,
'NDWI': ndwi
})
attempts += 1
df_samples = pd.DataFrame(valid_points)
print(f"Selesai! Berhasil mengumpulkan tepat {len(df_samples)} titik sampel valid setelah {attempts} percobaan.")# Skalakan band ke rentang 0-1 untuk indeks tingkat lanjut
b2_s = df_samples['B2_Blue'] / 10000.0
b4_s = df_samples['B4_Red'] / 10000.0
b8_s = df_samples['B8_NIR'] / 10000.0
# EVI
evi_denom = b8_s + 6 * b4_s - 7.5 * b2_s + 1
df_samples['EVI'] = np.clip(np.where(evi_denom == 0, np.nan, 2.5 * (b8_s - b4_s) / evi_denom), -1, 1)
# SAVI
df_samples['SAVI'] = ((df_samples['B8_NIR'] - df_samples['B4_Red']) / (df_samples['B8_NIR'] + df_samples['B4_Red'] + 0.5)) * 1.5
# Tangani nilai kosong jika ada
df_samples = df_samples.fillna(0)
print("Indeks vegetasi EVI dan SAVI berhasil dihitung!")# Definisikan path dataset raster
env_rasters = {
'Elevasi_mdpl': 'DEM_Konawe_Aligned.tif',
'Slope_persen': 'Slope_Konawe_Aligned.tif',
'CurahHujan_mmth': 'CHIRPS_Rainfall_Konawe_2015_2025.tif',
'SQ1_Hara': 'sq1_aligned.tif',
'SQ2_Perakaran': 'sq2_aligned.tif',
'SQ3_Oksigen': 'sq3_aligned.tif',
'SQ4_Toksisitas': 'sq4_aligned.tif'
}
coords = list(zip(df_samples['longitude'], df_samples['latitude']))
for col_name, file_name in env_rasters.items():
# Cari di folder output_FINAL atau GEE_Konawe
path = f'{GEE_FOLDER}/output_FINAL/{file_name}'
if not os.path.exists(path):
path = f'{GEE_FOLDER}/{file_name}'
print(f"Mengekstrak fitur: {col_name}...")
with rasterio.open(path) as src:
df_samples[col_name] = [val[0] for val in src.sample(coords)]
# Isi data kosong/NoData dengan 0
df_samples = df_samples.fillna(0)
print("\nSemua fitur spasial berhasil diekstrak!")Sel ini menerapkan scoring threshold berbasis hukum minimum Liebig. Menghasilkan label klasifikasi 3-kelas.
# Fungsi scoring berdasarkan kriteria ilmiah
def score_dem(v):
s = np.ones_like(v)
s = np.where((v>=700)&(v<=1000), 2, s)
s = np.where(((v>=400)&(v<700))|((v>=0)&(v<10)), 3, s)
s = np.where((v>=10)&(v<=400), 4, s)
return s
def score_slope(v):
s = np.ones_like(v)
s = np.where((v>15)&(v<=30), 2, s)
s = np.where((v>8)&(v<=15), 3, s)
s = np.where(v<=8, 4, s)
return s
def score_rainfall(v):
s = np.ones_like(v)
s = np.where(((v>=1500)&(v<2000))|((v>4000)&(v<=4500)), 2, s)
s = np.where(((v>=2000)&(v<2500))|((v>3500)&(v<=4000)), 3, s)
s = np.where((v>=2500)&(v<=3500), 4, s)
return s
def score_sq(v):
# SQ 1 (Hara) s/d SQ4 (Toksisitas)
s = np.full_like(v, np.nan)
s = np.where(v>=4, 1, s)
s = np.where(v==3, 2, s)
s = np.where(v==2, 3, s)
s = np.where(v==1, 4, s)
return s
# Terapkan scoring
s_dem = score_dem(df_samples['Elevasi_mdpl'].values)
s_slope = score_slope(df_samples['Slope_persen'].values)
s_rain = score_rainfall(df_samples['CurahHujan_mmth'].values)
s_sq1 = score_sq(df_samples['SQ1_Hara'].values)
s_sq2 = score_sq(df_samples['SQ2_Perakaran'].values)
s_sq3 = score_sq(df_samples['SQ3_Oksigen'].values)
s_sq4 = score_sq(df_samples['SQ4_Toksisitas'].values)
# Liebig's Law of the Minimum (Ambil nilai terendah/terburuk sebagai pembatas)
label_4_kelas = np.fmin(np.fmin(np.fmin(s_dem, s_slope), s_rain),
np.fmin(np.fmin(s_sq1, s_sq2), np.fmin(s_sq3, s_sq4)))
# Konversi ke sistem 3-kelas
# 4 (S1) & 3 (S2) -> 2 (Sesuai)
# 2 (S3) -> 1 (Sesuai Marginal)
# 1 (N) -> 0 (Tidak Sesuai)
df_samples['Label_Skor'] = label_4_kelas
df_samples['Label_Kelas'] = np.where(label_4_kelas >= 3, 2, np.where(label_4_kelas == 2, 1, 0))
# Map ke label teks pembacaan
label_map = {2: 'Sesuai', 1: 'Sesuai Marginal', 0: 'Tidak Sesuai'}
df_samples['Label_Teks'] = df_samples['Label_Kelas'].map(label_map)
print("\nDistribusi Label Hasil Pseudo-Sampling (10.000 titik):")
print(df_samples['Label_Teks'].value_counts())# Simpan dataset lengkap
csv_output_path = f'{GEE_FOLDER}/output_FINAL/dataset_nilam_10k_labeled.csv'
df_samples.to_csv(csv_output_path, index=False)
print(f"Dataset 10.000 sampel berlabel berhasil disimpan ke: {csv_output_path}")Sel ini melatih model Random Forest menggunakan 15 fitur dengan pembagian data 80:20 (train/test) secara seimbang (stratified split).
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report
# Pilih 15 fitur final
features = [
'B2_Blue', 'B3_Green', 'B4_Red', 'B8_NIR',
'NDVI', 'EVI', 'NDWI', 'SAVI',
'Elevasi_mdpl', 'Slope_persen', 'CurahHujan_mmth',
'SQ1_Hara', 'SQ2_Perakaran', 'SQ3_Oksigen', 'SQ4_Toksisitas'
]
X = df_samples[features]
y = df_samples['Label_Kelas']
# Split train & test (80:20)
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42, stratify=y
)
print(f"Melatih model Random Forest pada {len(X_train)} sampel latih...")
# Konfigurasi model Random Forest yang aman & efisien
rf_model = RandomForestClassifier(
n_estimators=100,
max_depth=12,
random_state=42,
n_jobs=-1
)
rf_model.fit(X_train, y_train)
# Evaluasi hasil
y_pred = rf_model.predict(X_test)
print("\n[HASIL EVALUASI MODEL - DATA TESTING (2.000 Sampel)]")
print(classification_report(y_test, y_pred, target_names=['Tidak Sesuai', 'Sesuai Marginal', 'Sesuai']))
# Feature Importance
importances = pd.Series(rf_model.feature_importances_, index=features).sort_values(ascending=False)
print("\n[TINGKAT KEPENTINGAN FITUR (FEATURE IMPORTANCE)]")
for name, val in importances.items():
print(f" - {name}: {val*100:.2f}%")