# [SKRIP 1.1-ArcMap: Memotong citra dan menghitung NDVI dengan ArcPy]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh ekstensi Spatial Analyst.
from __future__ import print_function   # agar print() berlaku di Python 2.7
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

red = Raster("B04.tif")
nir = Raster("B08.tif")
batas = "Batas_KPH.shp"

red_clip = ExtractByMask(red, batas)
nir_clip = ExtractByMask(nir, batas)

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.")