# -*- coding: utf-8 -*-
# [SKRIP 2.1-ArcMap: Rantai alat: sempadan sungai, titik panas, dan NDVI per KPH (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Zonal Statistics as Table dan Float butuh ekstensi Spatial Analyst.
import arcpy
from arcpy.sa import Raster, Float, ZonalStatisticsAsTable

DATA = r"D:\Latihan\paket-i4"
HASIL = r"D:\Latihan\hasil_i4"
arcpy.env.workspace = DATA + r"\shp"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

# Keluaran sementara di ArcMap: ruang kerja in_memory (di Pro namanya memory)
arcpy.Buffer_analysis("Sungai.shp", r"in_memory\sempadan", "30 Meters", dissolve_option="ALL")
arcpy.Intersect_analysis(["Batas_KPH.shp", r"in_memory\sempadan"], r"in_memory\irisan")
arcpy.SpatialJoin_analysis("Batas_KPH.shp", "Hotspot_Nasional.shp", r"in_memory\kph_titik",
                           "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="INTERSECT")   # menambah Join_Count

sempadan = {}
with arcpy.da.SearchCursor(r"in_memory\irisan", ["NAMA_KPH", "SHAPE@AREA"]) as k:
    for nama, luas in k:
        sempadan[nama] = sempadan.get(nama, 0) + luas / 10000.0
with arcpy.da.SearchCursor(r"in_memory\kph_titik", ["NAMA_KPH", "Join_Count", "SHAPE@AREA"]) as k:
    for nama, jml, luas in k:
        print("%s: luas %.2f ha, sempadan %.2f ha, titik panas %d"
              % (nama, luas / 10000.0, sempadan.get(nama, 0), jml))

# NDVI = (NIR - Merah) / (NIR + Merah); band 4 = NIR, band 1 = Merah
nir = Float(Raster(DATA + r"\Citra_KPH.tif\Band_4"))
merah = Float(Raster(DATA + r"\Citra_KPH.tif\Band_1"))
ndvi = (nir - merah) / (nir + merah)
ndvi.save(HASIL + r"\NDVI_KPH.tif")
ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", ndvi, r"in_memory\ndvi_zonal", "DATA", "MEAN")
with arcpy.da.SearchCursor(r"in_memory\ndvi_zonal", ["NAMA_KPH", "MEAN"]) as k:
    for nama, rata in k:
        print("%s: NDVI rata-rata %.3f" % (nama, rata))
arcpy.Delete_management("in_memory")      # bersihkan seluruh ruang in_memory
arcpy.CheckInExtension("Spatial")