# -- coding: utf-8 --
import arcpy
import os
import glob
import re
import datetime
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True

class BatchClip:
    def __init__(self, inws, outws, mask):
        self.inws = inws
        self.outws = outws
        self.mask = mask
        print(self.inws)

    def clip(self):
        if not os.path.exists(self.outws):
            os.mkdir(self.outws)
        rasters = glob.glob(os.path.join(self.inws, "*.tif"))
    
        shps = glob.glob(os.path.join(self.mask, "*.shp"))
        
        inwsname = self.inws[-5:]

         
        for shp in shps:
            shp_name = os.path.basename(shp)
            print(shp_name)
            
            for ras in rasters:
                ras_name = os.path.basename(ras)
                shpname = shp_name[-8:-4]
                rastername = ras_name[:8]
                nameT = inwsname + shpname+"_"+ rastername + ".tif"
                outname = os.path.join(self.outws, nameT)
                out_extract = arcpy.sa.ExtractByMask(ras, shp)
                out_extract.save(outname)
                

        print(" ----Batchclip down ----")

class RasterReclassify:
    def __init__(self, workspace):
        # Set workspace
        self.workspace = workspace 
        arcpy.env.workspace = self.workspace

    def reclassify_rasters(self, index):
        folder = arcpy.env.workspace[-4:]
        print(folder)
        folder_name = "reclassify" + folder
        folder_path = os.path.dirname(arcpy.env.workspace)
        output_path = os.path.join(folder_path, folder_name)

        if not os.path.exists(output_path):
            os.mkdir(output_path)
        print(output_path)

        rasterlist = arcpy.ListRasters("*", "tif")
        print("raster:",rasterlist)

        if rasterlist is not None:
            for raster in rasterlist:
                inRaster = arcpy.Raster(raster)
                print(inRaster)
                time = inRaster.name[5:18]
                print(time)
                out = os.path.join(output_path, time + "_"  + ".tif")
                print(out)

            # if index == "ndvi":
            #     # NDVI
            #     outCon = arcpy.sa.Con(inRaster <= 0.19, 1,
            #                           arcpy.sa.Con((inRaster > 0.19) & (inRaster <= 0.39), 2,
            #                                        arcpy.sa.Con((inRaster > 0.39) & (inRaster <= 0.55), 3,
            #                                                     arcpy.sa.Con((inRaster > 0.55) & (inRaster <= 0.75), 4,
            #                                                                  arcpy.sa.Con((inRaster > 0.75), 5)))))

            # elif index == "gndvi":
                # GNDVI
                outCon = arcpy.sa.Con(inRaster <= 0.2, 1,
                arcpy.sa.Con((inRaster > 0.2) & (inRaster <= 0.4), 2,
                arcpy.sa.Con((inRaster > 0.4) & (inRaster <= 0.6), 3,
                arcpy.sa.Con((inRaster > 0.6) & (inRaster <= 0.8), 4,
                arcpy.sa.Con((inRaster > 0.8), 5)))))

            # elif index == "vci":
            #     # VCI
            #     outCon = arcpy.sa.Con(inRaster <= 0.15, 1,
            #                           arcpy.sa.Con((inRaster > 0.15) & (inRaster <= 0.3), 2,
            #                                        arcpy.sa.Con((inRaster > 0.3) & (inRaster <= 0.45), 3,
            #                                                     arcpy.sa.Con((inRaster > 0.45) & (inRaster <= 0.6), 4,
            #                                                                  arcpy.sa.Con((inRaster > 0.6), 5)))))

            # elif index == "chlor":
            #     # Chlor
            #     outCon = arcpy.sa.Con(inRaster <= 12, 1,
            #                           arcpy.sa.Con((inRaster > 12) & (inRaster <= 21), 2,
            #                                        arcpy.sa.Con((inRaster > 21) & (inRaster <= 32), 3,
            #                                                     arcpy.sa.Con((inRaster > 32) & (inRaster <= 43), 4,
            #                                                                  arcpy.sa.Con((inRaster > 43), 5)))))

                outCon.save(out)
print("---raclassify down---")

class RasterToPolygonConverter:
    def __init__(self, input_raster_path, output_shp_path):
        self.input_raster_path = input_raster_path
        self.output_shp_path = output_shp_path

    def convert_raster_to_polygon(self):
        # Set workspace
        arcpy.env.workspace = self.input_raster_path
        
        # Create output folder if it doesn’t exist
        if not os.path.exists(self.output_shp_path):
            os.mkdir(self.output_shp_path)

        # List raster files
        rasters = arcpy.ListRasters("*", "tif")
        print("rasterstopolygon:",rasters)
        for raster in rasters:
            # Extract catchment name from raster filename
            catchment_name = raster[:8]
            print(catchment_name)

            # Set output shapefile path
            shp_out_name = os.path.join(self.output_shp_path, catchment_name)
            print(shp_out_name)

            # Convert raster to polygon
            arcpy.RasterToPolygon_conversion(raster, shp_out_name, "true", "Value")
            print("finished")
print("---polygon down---")

class Dissolve:
    def init(self):
        pass

    def dissolve_shapefiles(self, input_workspace, output_workspace):
        arcpy.env.workspace = input_workspace
        shps = arcpy.ListFiles("*.shp")
        
        if shps is not None:
            for shp in shps:
                catchment_name = shp[:8]
                shp_out_name = os.path.join(output_workspace, catchment_name)

                arcpy.Dissolve_management(shp, shp_out_name, "gridcode")
                print("Finished dissolving", catchment_name)


class FieldCalculator:
    def __init__(self, input_polygon):
        self.input_polygon = input_polygon

    def calculate_fields(self):
        current_year = datetime.datetime.now().year
        if not os.path.exists(self.input_polygon):
            return

        path_list = os.listdir(self.input_polygon)

        for i in range(len(path_list)):
            landid = path_list[i][-4:]
            print(landid)

            path = os.path.join(self.input_polygon, path_list[i])
            shps = glob.glob(os.path.join(path, "*.shp"))

            for shp in shps:
                print(shp)
                catchment_name = shp[-12:-4]

                print(catchment_name)

                date = int(catchment_name)

                # arcpy.DeleteField_management(shp, [“LandId”])

                field_names = [f.name for f in arcpy.ListFields(shp)]

                if 'date' not in field_names:
                    arcpy.AddField_management(shp, 'date', 'Long')
                    arcpy.CalculateField_management(shp, 'date', expression=date, expression_type="python", code_block="")

                if 'landId' not in field_names:
                    arcpy.AddField_management(shp, 'landId', 'Short’')
                    arcpy.CalculateField_management(shp, 'landId', 0, expression_type="python")

                if 'landName' not in field_names:
                    arcpy.AddField_management(shp, 'landName', 'TEXT')

                if 'crop' not in field_names:
                    arcpy.AddField_management(shp, 'crop', 'TEXT')

                if 'Area' not in field_names:
                    arcpy.AddField_management(shp, 'Area', 'Double')

                if 'Area_mu' not in field_names:
                    arcpy.AddField_management(shp, 'Area_mu', 'Double')

                if 'year' not in field_names:
                    arcpy.AddField_management(shp, 'year', 'Double')
                arcpy.CalculateField_management(shp, 'year', current_year, expression_type="python")
print("---dbf down---")


def main():
    # Set input and output paths
    inws = r"F:\GEE\wurenlongchang\tif\wurenlongchang_GDVI"
    outws = r"F:\GEE\wurenlongchang\mask"
    mask = r"F:\GEE\wurenlongchang\1_shp"

    # Batch clip rasters
    batch_clip = BatchClip(inws, outws, mask)
    batch_clip.clip()

    # Reclassify rasters
    reclassify_workspace = os.path.join(outws)
    rr = RasterReclassify(reclassify_workspace)
    rr.reclassify_rasters("ndvi")
    # rr.reclassify_rasters("gndvi")
    # # rr.reclassify_rasters("vci")
    # rr.reclassify_rasters("chlor")

    # Convert rasters to polygons
    polygon_workspace = os.path.join(outws)
    folder = polygon_workspace[-4:]
    folder_name = "polygon" + folder
    output_polygon_path = os.path.join(outws, folder_name)
    converter = RasterToPolygonConverter(polygon_workspace, output_polygon_path)
    converter.convert_raster_to_polygon()

    # Dissolve polygons
    dissolve_workspace = os.path.join(outws)
    folder = dissolve_workspace[-4:]
    folder_name = "dissolve" + folder
    output_dissolve_path = os.path.join(outws, folder_name)
    dissolve = Dissolve()
    dissolve.dissolve_shapefiles(reclassify_workspace, output_dissolve_path)

    # Calculate fields
    field_calculator = FieldCalculator(output_dissolve_path)
    field_calculator.calculate_fields()

if __name__ == "__main__":
    main()
class RasterReclassify的输出文件作为class RasterToPolygonConverter:的输入文件如何处理