# -- coding: utf-8 --

import arcpy
import os
import glob
import re
import datetime

arcpy.CheckOutExtension('Spatial')

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

    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"))
        print(rasters,shps)

        for shp in shps:
            for ras in rasters:
                nameT = os.path.basename(ras) + ".tif"
                outname = os.path.join(self.outws, nameT)
                out_extract = arcpy.sa.ExtractByMask(ras, shp)
                out_extract.save(outname)
                print(os.path.basename(ras) + " ---- down")

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

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

    def reclassify_rasters(self, index):
        folder = arcpy.env.workspace[-4:]
        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")

        if rasterlist is not None:
            for raster in rasterlist:
                inRaster = arcpy.Raster(raster)
                print(inRaster)
                time = inRaster.name[:8]
                print(time)
                out = os.path.join(output_path, time + "_" + index + ".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")
        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\wurenlongchang_ndvi"
    outws = r"F:\GEE\wurenlongchang"
    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, "ndvi")
    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, "reclassifyndvi")
    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, "polygonndvi")
    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()
文件读不了怎么办