# -- coding: utf-8 --
import arcpy
import os
import glob
import re
import datetime
import pandas as pd
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.outws)

    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:]
        print(inwsname)

         
        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.join(arcpy.env.workspace)  #dirname
        raclassify_output_path = os.path.join(folder_path, folder_name)
        # print(raclassify_output_path)

        if not os.path.exists(raclassify_output_path):
            os.mkdir(raclassify_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(raclassify_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
        # print(arcpy.env.workspace)
        
        # 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[:13]
            # 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, input_workspace, output_workspace):
        self.input_workspace = input_workspace
        self.output_workspace = output_workspace
        

    def dissolve_shapefiles(self):
        arcpy.env.workspace = self.input_workspace
        shps = arcpy.ListFiles("*.shp")
        
        if shps is not None:
            for shp in shps:
                catchment_name = shp[:13]
                print(catchment_name)
                disslove_out_name = os.path.join(self.output_workspace,catchment_name)
                # print(disslove_out_name)   
                arcpy.Dissolve_management(shp, disslove_out_name, "gridcode")
                # print("Finished dissolving")


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

    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(self.input_polygon, "*.shp"))
        print(shps)
    

        for shp in shps:
            catchment_name = shp[-12:-4]
            landid = shp[-17:-13]
            print(landid)
            landId = int(landid)
            print(landId)

            if landid is not None:
                filtered_df = self.df[self.df['ID'] == landId]

                df_crop = filtered_df[u'种植作物'].iloc[0]
                crop = '"' + df_crop + '"'

                # crop = df_crop.encode('unicode_escape').decode('utf-8')


                df_landName = filtered_df[u'项目名称'].iloc[0]
                landName = '"' + df_landName + '"'

                # landName = df_landName.encode('unicode_escape').decode('utf-8')
                
                print(crop)
                print(landName)
            

            date = int(catchment_name)

            arcpy.DeleteField_management(shp, ["id"])

            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', landid, expression_type="python")

            if 'landName' not in field_names:
                arcpy.AddField_management(shp, 'landName', 'TEXT')
            arcpy.CalculateField_management(shp, 'landName', expression=landName, expression_type="python")

            if 'crop' not in field_names:
                arcpy.AddField_management(shp, 'crop', 'TEXT')
            arcpy.CalculateField_management(shp, 'crop', expression=crop, expression_type="python")

            if 'Area' not in field_names:
                arcpy.AddField_management(shp, 'Area', 'Double')
            arcpy.CalculateField_management(shp, 'Area', "!SHAPE.area!", expression_type="python")

            if 'Area_mu' not in field_names:
                arcpy.AddField_management(shp, 'Area_mu', 'Double')
            arcpy.CalculateField_management(shp, 'Area_mu', "round( !Area! *0.0015,4)", expression_type="python")

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

class Merge:
    def __init__(self,workfloder,Merge_output):
        self.workfloder = workfloder
        self.Merge_output = Merge_output
    
    def merge(self):
        arcpy.env.workspace = self.workfloder  #change this
        merge_name = self.workfloder[-4:]
        if not os.path.exists(self.Merge_output):
            os.mkdir(self.Merge_output)
        
        output_merge_path = os.path.join(self.Merge_output, merge_name + ".shp")
        shplist =  arcpy.ListFeatureClasses('*.shp')
        arcpy.Merge_management(shplist, output_merge_path)

    


def main():
    # Set input and output paths
    inws = r"F:\GEE\chenlongwen\tif\chenlongwen_ndvi"
    outws = r"F:\GEE\chenlongwen\growth"
    mask = r"F:\GEE\chenlongwen\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()

    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)
    raclassify_output_path = os.path.join(outws, "reclassify" + folder)  # 定义 raclassify_output_path 变量
    converter = RasterToPolygonConverter(raclassify_output_path, 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)
    if not os.path.exists(output_dissolve_path):
        os.mkdir(output_dissolve_path)
    dissolve = Dissolve(output_polygon_path, output_dissolve_path)
    dissolve.dissolve_shapefiles()

    # Calculate fields
    df = pd.read_excel(ur'F:\GEE\z_other\\遥感地块信息.xlsx')
    field_calculator = FieldCalculator(output_dissolve_path,df)
    field_calculator.calculate_fields()

    #Merge  
    folder_name = "merge" 
    output_merge_path = os.path.join(outws, folder_name)
    meger = Merge(output_dissolve_path,output_merge_path)
    meger.merge()


arcpy.ClearWorkspaceCache_management()

if __name__ == "__main__":
    main()

如何使用index