arcgis通过栅格影像数据批量裁剪shp


import arcpy
from arcpy import env

from util import get_imlist
import os
import math
import shutil

#输出定义polygon(extend)
polygons = r"polygons.shp"
# Set workspace
Raster_path = r'Input'
env.workspace = r"shapefile"
Out_path = r'Clip'
MBR_path = r'MBR'
BB_path = r'BoundingBox'
Raster_TYPE = 'tiff'
def check_dir(path):
    if os.path.exists(path):
        shutil.rmtree(path)
        os.mkdir(path)
    else:
        os.mkdir(path)


if __name__ == '__main__':

    check_dir(Out_path)
    check_dir(MBR_path)
    check_dir(BB_path)
    #系统默认为平面坐标系,如果栅格数据是经纬度需要定义投影坐标系
    #4326表示Coordinate Systems/Geographic Coodinate Systems/World/WGS 1984.prj
    sr = arcpy.SpatialReference(4326)
    in_features = "file.shp"
    rasters = get_imlist(Raster_path, Raster_TYPE)
    for raster in rasters:
        in_raster = os.path.join(Raster_path, raster)
        pnt_array = arcpy.Array()
        extent = arcpy.Raster(in_raster).extent
        pnt_array.add(extent.lowerLeft)
        pnt_array.add(extent.lowerRight)
        pnt_array.add(extent.upperRight)
        pnt_array.add(extent.upperLeft)

        #定义坐标参考系的polygon,输出裁剪的extend范围的polygon
        #poly = arcpy.Polygon(pnt_array,sr)
        arcpy.CreateFeatureclass_management(*os.path.split(polygons), geometry_type="Polygon",spatial_reference=sr)
        icursor = arcpy.da.InsertCursor(polygons, ["SHAPE@"])
        polygon = arcpy.Polygon(pnt_array,sr)
        icursor.insertRow([polygon])



        poly = arcpy.Polygon(pnt_array)
        filename = raster.strip('.tiff') + '.shp'
        out_feature_class = os.path.join(Out_path, filename)
        xy_tolerance = ""

        # Execute Clip
        inFeatures = arcpy.Clip_analysis(in_features, poly, out_feature_class, xy_tolerance)

        # make MBR
        outFeatureClass = os.path.join(MBR_path,filename)
        MB = arcpy.MinimumBoundingGeometry_management(inFeatures, outFeatureClass,
                                                 "ENVELOPE", "NONE")

        #get BoundingBox

        bb_file = open(os.path.join(BB_path,raster.strip('.tiff') + '.txt'), 'w')
        # Enter for loop for each feature
        #
        xMax = extent.XMax
        xMin = extent.XMin
        yMax = extent.YMax
        yMin = extent.YMin

        for row in arcpy.da.SearchCursor(MB, ["OID@", "SHAPE@"]):
            # Print the current multipoint's ID
            #
            print("Feature {}:".format(row[0]))
            partnum = 0

            # Step through each part of the feature
            bb_extent = row[1].extent
            xBB_min = math.ceil(bb_extent.XMin - xMin)
            xBB_max = math.floor(bb_extent.XMax - xMin)
            yBB_min = math.ceil(bb_extent.YMin - yMin)
            yBB_max = math.floor(bb_extent.YMax - yMin)
            result = r'{},{},{},{},0'.format(xBB_min,yBB_min,xBB_max,yBB_max)
            bb_file.write(result+'\n')
        bb_file.close()

转载自:https://blog.csdn.net/cugxyy6/article/details/79284533

You may also like...

退出移动版