# [SKRIP 1.1-Pro: Memotong citra dan menghitung NDVI dengan ArcPy]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh ekstensi Spatial Analyst atau Image Analyst.
import arcpy
from arcpy.sa import Raster, Float, ExtractByMask

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\data_kph"      # ganti dengan folder Anda

# 1. Masukan: band Red, band NIR, dan batas area (CRS harus sama, lihat 1.2)
red = Raster("B04.tif")
nir = Raster("B08.tif")
batas = "Batas_KPH.shp"

# 2. Potong kedua band sesuai batas
red_clip = ExtractByMask(red, batas)
nir_clip = ExtractByMask(nir, batas)

# 3. NDVI = (NIR - RED) / (NIR + RED)
ndvi = (Float(nir_clip) - Float(red_clip)) / (Float(nir_clip) + Float(red_clip))
ndvi.save("NDVI_KPH.tif")

arcpy.CheckInExtension("Spatial")
print("Selesai: NDVI_KPH.tif tersimpan.")