import rasterio
import numpy as np
import os
import glob
import [Link] as plt
from [Link] import ListedColormap
# Suppress division by zero and invalid value warnings
[Link](all='warn', divide='ignore', invalid='ignore')
# 1. Tentukan folder tempat script dijalankan (sama dengan lokasi .py)
script_dir = [Link]([Link](__file__))
# 2. Cari file TIFF di folder tersebut
tiff_files = [Link]([Link](script_dir, "*.tif")) +
[Link]([Link](script_dir, "*.tiff"))
# Cari file dengan 3 (RGB) atau 4 (RGBA) band
rgb_file = None
for tif in tiff_files:
try:
with [Link](tif) as src:
num_bands = [Link]
if num_bands in (3, 4):
rgb_file = tif
break
except [Link]:
continue # Lewati file yang tidak bisa dibaca
if not rgb_file:
print(f"Error: Tidak ditemukan file TIFF dengan 3 (RGB) atau 4 (RGBA) band di
folder {script_dir}.")
exit(1)
print(f"File TIFF yang digunakan: {rgb_file}")
# 3. Baca citra dan periksa deskripsi band
with [Link](rgb_file) as src:
rgb_data = [Link]() # Baca semua band
rgb_meta = [Link]()
num_bands = [Link]
band_descriptions = [Link]
print("Deskripsi band:", band_descriptions if band_descriptions else "Tidak ada
deskripsi band")
# Ekstrak band
if num_bands == 3:
red = rgb_data[0].astype(float)
green = rgb_data[1].astype(float)
blue = rgb_data[2].astype(float)
alpha = None
elif num_bands == 4:
red = rgb_data[0].astype(float)
green = rgb_data[1].astype(float)
blue = rgb_data[2].astype(float)
alpha = rgb_data[3].astype(float)
else:
print(f"Error: File {rgb_file} memiliki {num_bands} band, hanya mendukung 3
(RGB) atau 4 (RGBA) band.")
exit(1)
# Mask area dengan nodata=0 atau alpha=0 (jika RGBA)
mask = (red == 0) & (green == 0) & (blue == 0)
if alpha is not None:
mask = mask | (alpha == 0)
red = [Link](mask, [Link], red)
green = [Link](mask, [Link], green)
blue = [Link](mask, [Link], blue)
# Update metadata untuk output
rgb_meta_float = rgb_meta.copy()
rgb_meta_float.update({
"count": 1,
"dtype": rasterio.float32,
"nodata": [Link],
})
rgb_meta_uint8 = rgb_meta.copy()
rgb_meta_uint8.update({
"count": 1,
"dtype": rasterio.uint8,
"nodata": 0,
})
# 4. Hitung VARI tanpa kelandaian
denominator = green + red - blue
denom_threshold = 1e-8
vari = (green - red) / denominator
vari = [Link]([Link](vari) | [Link](vari) | ([Link](denominator) <
denom_threshold), [Link], vari)
# Stretch histogram VARI ke [-1, 1] menggunakan min/max
valid_vari = vari[~[Link](vari)]
stats = []
if valid_vari.size > 0:
vari_min, vari_max = [Link](valid_vari), [Link](valid_vari)
print(f"Rentang VARI sebelum stretching (min/max): {vari_min:.3f} hingga
{vari_max:.3f}")
[Link](f"VARI sebelum stretching - Min: {vari_min:.3f}, Max:
{vari_max:.3f}, Mean: {[Link](valid_vari):.3f}")
[Link](f"VARI sebelum stretching - Jumlah nilai unik:
{len([Link](valid_vari))}")
if vari_max > vari_min:
vari_stretched = -1 + 2 * (vari - vari_min) / (vari_max - vari_min)
vari_stretched = [Link]([Link](vari), [Link], vari_stretched)
vari_stretched = [Link](vari_stretched, -1, 1)
else:
vari_stretched = [Link]([Link](vari), [Link], 0.0)
valid_vari_stretched = vari_stretched[~[Link](vari_stretched)]
if valid_vari_stretched.size > 0:
[Link](f"VARI setelah stretching - Min:
{[Link](valid_vari_stretched):.3f}, Max: {[Link](valid_vari_stretched):.3f},
Mean: {[Link](valid_vari_stretched):.3f}")
[Link](f"VARI setelah stretching - Jumlah nilai unik:
{len([Link](valid_vari_stretched))}")
else:
[Link]("VARI setelah stretching - Tidak ada data valid")
else:
vari_stretched = np.full_like(vari, [Link])
[Link]("VARI sebelum stretching - Tidak ada data valid")
[Link]("VARI setelah stretching - Tidak ada data valid")
# Simpan [Link]
vari_path = [Link](script_dir, "[Link]")
if [Link](vari_path):
try:
[Link](vari_path)
except PermissionError:
print(f"Error: Tidak dapat menghapus {vari_path}. Pastikan file tidak
sedang digunakan.")
exit(1)
with [Link](vari_path, "w", **rgb_meta_float) as dst:
[Link](vari_stretched, 1)
# 5. Hitung TGI
tgi = green - 0.39 * red - 0.61 * blue
tgi = [Link]([Link](tgi), [Link], tgi)
# Stretch histogram TGI ke [-1, 1] menggunakan min/max
valid_tgi = tgi[~[Link](tgi)]
if valid_tgi.size > 0:
tgi_min, tgi_max = [Link](valid_tgi), [Link](valid_tgi)
print(f"Rentang TGI sebelum stretching (min/max): {tgi_min:.3f} hingga
{tgi_max:.3f}")
[Link](f"TGI sebelum stretching - Min: {tgi_min:.3f}, Max: {tgi_max:.3f},
Mean: {[Link](valid_tgi):.3f}")
[Link](f"TGI sebelum stretching - Jumlah nilai unik:
{len([Link](valid_tgi))}")
if tgi_max > tgi_min:
tgi_stretched = -1 + 2 * (tgi - tgi_min) / (tgi_max - tgi_min)
tgi_stretched = [Link]([Link](tgi), [Link], tgi_stretched)
tgi_stretched = [Link](tgi_stretched, -1, 1)
else:
tgi_stretched = [Link]([Link](tgi), [Link], 0.0)
valid_tgi_stretched = tgi_stretched[~[Link](tgi_stretched)]
if valid_tgi_stretched.size > 0:
[Link](f"TGI setelah stretching - Min:
{[Link](valid_tgi_stretched):.3f}, Max: {[Link](valid_tgi_stretched):.3f},
Mean: {[Link](valid_tgi_stretched):.3f}")
[Link](f"TGI setelah stretching - Jumlah nilai unik:
{len([Link](valid_tgi_stretched))}")
else:
[Link]("TGI setelah stretching - Tidak ada data valid")
else:
tgi_stretched = np.full_like(tgi, [Link])
[Link]("TGI sebelum stretching - Tidak ada data valid")
[Link]("TGI setelah stretching - Tidak ada data valid")
# Simpan [Link]
tgi_path = [Link](script_dir, "[Link]")
if [Link](tgi_path):
try:
[Link](tgi_path)
except PermissionError:
print(f"Error: Tidak dapat menghapus {tgi_path}. Pastikan file tidak sedang
digunakan.")
exit(1)
with [Link](tgi_path, "w", **rgb_meta_float) as dst:
[Link](tgi_stretched, 1)
# 6. Klasifikasi berdasarkan RGB dan VARI/TGI
vari_classes = np.zeros_like(vari, dtype=np.uint8)
tgi_classes = np.zeros_like(tgi, dtype=np.uint8)
# Klasifikasi air (class 1)
vari_classes[mask] = 1
tgi_classes[mask] = 1
# Klasifikasi bukaan tanah (class 2)
soil_mask = (~mask) & (red / (green + 1e-8) > 1.2) & (vari_stretched < -0.5) &
(tgi_stretched < -0.5)
vari_classes[soil_mask] = 2
tgi_classes[soil_mask] = 2
# Klasifikasi hutan (class 3-8)
forest_mask = (~mask) & (~soil_mask)
class_bounds = [Link](-1, 1, 7) # 6 kelas
for i, (lower, upper) in enumerate(zip(class_bounds[:-1], class_bounds[1:])):
class_id = i + 3
vari_classes[forest_mask & (vari_stretched >= lower) & (vari_stretched <
upper)] = class_id
tgi_classes[forest_mask & (tgi_stretched >= lower) & (tgi_stretched < upper)] =
class_id
vari_classes[forest_mask & (vari_stretched >= class_bounds[-1])] = 8
tgi_classes[forest_mask & (tgi_stretched >= class_bounds[-1])] = 8
# Set nodata areas yang tidak termasuk air ke 0
vari_classes[[Link](vari_stretched) & (~mask)] = 0
tgi_classes[[Link](tgi_stretched) & (~mask)] = 0
# Simpan VARI_class.tif
vari_class_path = [Link](script_dir, "VARI_class.tif")
if [Link](vari_class_path):
try:
[Link](vari_class_path)
except PermissionError:
print(f"Error: Tidak dapat menghapus {vari_class_path}. Pastikan file tidak
sedang digunakan.")
exit(1)
with [Link](vari_class_path, "w", **rgb_meta_uint8) as dst:
[Link](vari_classes, 1)
# Simpan TGI_class.tif
tgi_class_path = [Link](script_dir, "TGI_class.tif")
if [Link](tgi_class_path):
try:
[Link](tgi_class_path)
except PermissionError:
print(f"Error: Tidak dapat menghapus {tgi_class_path}. Pastikan file tidak
sedang digunakan.")
exit(1)
with [Link](tgi_class_path, "w", **rgb_meta_uint8) as dst:
[Link](tgi_classes, 1)
# 7. Visualisasi peta hasil
colors = ['#00B7EB', '#FF4500', '#228B22', '#32CD32', '#ADFF2F', '#FFFF00',
'#FFA500', '#FF69B4'] # Air: biru, Tanah: oranye, Hutan: hijau/kuning
custom_cmap = ListedColormap(colors)
fig, (ax1, ax2, ax3, ax4) = [Link](1, 4, figsize=(24, 5))
# Peta VARI
im1 = [Link](vari_stretched, cmap='RdYlGn', vmin=-1, vmax=1)
ax1.set_title('VARI (Stretched)')
[Link](im1, ax=ax1, label='Nilai VARI')
# Peta TGI
im2 = [Link](tgi_stretched, cmap='RdYlGn', vmin=-1, vmax=1)
ax2.set_title('TGI (Stretched)')
[Link](im2, ax=ax2, label='Nilai TGI')
# Peta kelas VARI
im3 = [Link](vari_classes, cmap=custom_cmap, vmin=1, vmax=8)
ax3.set_title('Klasifikasi VARI')
[Link](im3, ax=ax3, label='Kelas (1=Air, 2=Tanah, 3-8=Hutan)')
# Peta kelas TGI
im4 = [Link](tgi_classes, cmap=custom_cmap, vmin=1, vmax=8)
ax4.set_title('Klasifikasi TGI')
[Link](im4, ax=ax4, label='Kelas (1=Air, 2=Tanah, 3-8=Hutan)')
plt.tight_layout()
plot_path = [Link](script_dir, "VARI_TGI_plot.png")
[Link](plot_path, dpi=300, bbox_inches='tight')
[Link]()
# 8. Visualisasi histogram
fig, ((ax1, ax2), (ax3, ax4)) = [Link](2, 2, figsize=(12, 8))
# Histogram VARI sebelum stretching
if valid_vari.size > 0:
[Link](valid_vari, bins=100, color='blue', alpha=0.7)
ax1.set_title('Histogram VARI (Sebelum Stretching)')
ax1.set_xlabel('Nilai VARI')
ax1.set_ylabel('Frekuensi')
# Histogram VARI setelah stretching
if valid_vari_stretched.size > 0:
[Link](valid_vari_stretched, bins=100, color='blue', alpha=0.7)
ax2.set_title('Histogram VARI (Setelah Stretching)')
ax2.set_xlabel('Nilai VARI Stretched')
ax2.set_ylabel('Frekuensi')
# Histogram TGI sebelum stretching
if valid_tgi.size > 0:
[Link](valid_tgi, bins=100, color='orange', alpha=0.7)
ax3.set_title('Histogram TGI (Sebelum Stretching)')
ax3.set_xlabel('Nilai TGI')
ax3.set_ylabel('Frekuensi')
# Histogram TGI setelah stretching
if valid_tgi_stretched.size > 0:
[Link](valid_tgi_stretched, bins=100, color='orange', alpha=0.7)
ax4.set_title('Histogram TGI (Setelah Stretching)')
ax4.set_xlabel('Nilai TGI Stretched')
ax4.set_ylabel('Frekuensi')
plt.tight_layout()
hist_plot_path = [Link](script_dir, "VARI_TGI_histograms.png")
[Link](hist_plot_path, dpi=300, bbox_inches='tight')
[Link]()
# 9. Simpan statistik
stats_path = [Link](script_dir, "[Link]")
with open(stats_path, 'w') as f:
for stat in stats:
[Link](stat + '\n')
for i in range(1, 9):
vari_count = [Link](vari_classes == i)
tgi_count = [Link](tgi_classes == i)
[Link](f"Kelas {i} (VARI): {vari_count} piksel\n")
[Link](f"Kelas {i} (TGI): {tgi_count} piksel\n")
print(f"Hasil VARI disimpan di: {vari_path}")
print(f"Hasil TGI disimpan di: {tgi_path}")
print(f"Hasil klasifikasi VARI disimpan di: {vari_class_path}")
print(f"Hasil klasifikasi TGI disimpan di: {tgi_class_path}")
print(f"Plot peta disimpan di: {plot_path}")
print(f"Plot histogram disimpan di: {hist_plot_path}")
print(f"Statistik disimpan di: {stats_path}")
