# [SKRIP 3.1-Pro: Batch: NDVI semua ubin dan potongan CHM per KPH dengan ArcPy]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh ekstensi Spatial Analyst.
import os
import arcpy
from arcpy.sa import Raster, Float, ExtractByMask

DATA = r"D:\Latihan\paket-i4"
HASIL = r"D:\Latihan\hasil_i4"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
os.makedirs(HASIL + r"\NDVI_ubin", exist_ok=True)

# A. Semua ubin dalam satu folder (band 1 = NIR, band 2 = Merah)
arcpy.env.workspace = DATA + r"\Ortofoto_Mentah"
daftar = arcpy.ListRasters("*", "TIF")
for i, nama in enumerate(daftar, start=1):
    b1 = Float(Raster(os.path.join(DATA, "Ortofoto_Mentah", nama, "Band_1")))
    b2 = Float(Raster(os.path.join(DATA, "Ortofoto_Mentah", nama, "Band_2")))
    ((b1 - b2) / (b1 + b2)).save(os.path.join(HASIL, "NDVI_ubin", "NDVI_" + nama))
    print("%d/%d selesai: %s" % (i, len(daftar), nama))

# B. Satu putaran per KPH: pilih satu poligon, potong CHM dengannya
chm = Raster(DATA + r"\DSM_Drone.tif") - Raster(DATA + r"\DTM_Drone.tif")
batas = DATA + r"\shp\Batas_KPH.shp"
with arcpy.da.SearchCursor(batas, ["NAMA_KPH"]) as k:
    nama_kph = [b[0] for b in k]
for nama in nama_kph:
    arcpy.analysis.Select(batas, r"memory\satu_kph", "NAMA_KPH = '%s'" % nama)
    keluar = os.path.join(HASIL, "CHM_" + nama.replace(" ", "_") + ".tif")
    ExtractByMask(chm, r"memory\satu_kph").save(keluar)
    print("CHM tersimpan untuk", nama)
arcpy.CheckInExtension("Spatial")