# -*- coding: utf-8 -*-
# [SKRIP 7.1-ArcMap: Proyek mini terpadu dengan ArcPy: ringkasan per KPH (Python 2.7)]
# Penulis: Badar Mubarok Yogaswara
# [CEK] Cocokkan dengan versi Anda. Butuh Spatial Analyst.
import csv
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")

try:
    # I1: sempadan 30 m dan titik panas
    arcpy.Buffer_analysis("Sungai.shp", r"in_memory\sem", "30 Meters", dissolve_option="ALL")
    arcpy.Intersect_analysis(["Batas_KPH.shp", r"in_memory\sem"], r"in_memory\iris")
    arcpy.SpatialJoin_analysis("Batas_KPH.shp", "Hotspot_Nasional.shp", r"in_memory\ringkas",
                               "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="INTERSECT")
    # I2: NDVI dan CHM, lalu rata-rata per KPH
    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)
    chm = Raster(DATA + r"\DSM_Drone.tif") - Raster(DATA + r"\DTM_Drone.tif")
    ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", ndvi, r"in_memory\z_ndvi", "DATA", "MEAN")
    ZonalStatisticsAsTable("Batas_KPH.shp", "NAMA_KPH", chm, r"in_memory\z_chm", "DATA", "MEAN")

    def tabel(jalur, kolom):
        with arcpy.da.SearchCursor(jalur, ["NAMA_KPH", kolom]) as k:
            return dict((n, v) for n, v in k)

    sempadan = {}
    with arcpy.da.SearchCursor(r"in_memory\iris", ["NAMA_KPH", "SHAPE@AREA"]) as k:
        for n, luas in k:
            sempadan[n] = sempadan.get(n, 0) + luas / 10000.0
    n_rata, c_rata = tabel(r"in_memory\z_ndvi", "MEAN"), tabel(r"in_memory\z_chm", "MEAN")
    with open(HASIL + r"\ringkasan_kph_arcmap.csv", "wb") as f:      # Python 2: mode "wb" untuk csv
        w = csv.writer(f)
        w.writerow(["NAMA_KPH", "LUAS_HA", "SEMPADAN_HA", "HOTSPOT", "NDVI_RATA", "CHM_RATA_M"])
        with arcpy.da.SearchCursor(r"in_memory\ringkas", ["NAMA_KPH", "SHAPE@AREA", "Join_Count"]) as k:
            for n, luas, jml in k:
                w.writerow([n, round(luas / 10000.0, 2), round(sempadan.get(n, 0), 2), jml,
                            round(n_rata[n], 3), round(c_rata[n], 2)])
    print("Ringkasan tersimpan.")
except arcpy.ExecuteError:
    print("Alat gagal: %s" % arcpy.GetMessages(2))
finally:
    arcpy.CheckInExtension("Spatial")