Arcpy批量裁剪栅格数据并以图幅号命名文件分幅输出

之前应公司需求需要一个能批量裁剪输出的功能,上网研究了一下各路大神的方法,然后自己根据需要写了个程序,这里记录一下。


前言

基于ArcGIS的arcpy二次开发

一、数据准备

待裁剪影像、标准图框shp文件、ArcGIS10.2、Pycharm

二、代码

代码如下(示例):

import arcpy
import os, time

arcpy.CheckOutExtension("spatial")
typeall = {"tif": "tif", "img": "img", "bmp": "bmp", "jpg": "jpg", "png": "png"}
# arcpy.env.workspace = 'D:/yx'
# clip_Shp = 'D:/1/TFH.shp'
# fieldName = 'TFH'
# output_Path = 'D:/TF/ClipRasterResult'
# type1 = 'tif'
# workSpeed = 'HighSpeed'
# mosaic_first = 'NO'
arcpy.env.workspace = arcpy.GetParameterAsText(0)
type1 = arcpy.GetParameterAsText(1)
mosaic_first = arcpy.GetParameterAsText(2)
clip_Shp = arcpy.GetParameterAsText(3)
fieldName = arcpy.GetParameterAsText(4)
output_Path = arcpy.GetParameterAsText(5) + "\\ClipRasterResult"
workSpeed = arcpy.GetParameterAsText(6)
expression = arcpy.AddFieldDelimiters(clip_Shp, fieldName) + " <> ''"
type1 = typeall[type1]
rasters = arcpy.ListRasters("*", type1)


def create_dir(path):
    if os.path.exists(path):
        clear_folder(path)
    else:
        os.mkdir(path)


def clear_folder(root_path):
    files_list = os.listdir(root_path)
    for f in files_list:
        file_path = os.path.join(root_path, f)
        if os.path.isdir(file_path):
            clear_folder(file_path)
        else:
            os.remove(file_path)
    arcpy.AddMessage('clear ' + root_path)


def bool_point(X, Y, XMin, YMin, XMax, YMax):
    if XMin < X < XMax and YMin < Y < YMax:
        return True
    else:
        return False


def bool_extents(extent1, extent2):
    extent1_XMin = extent1.XMin
    extent1_XMax = extent1.XMax
    extent1_YMin = extent1.YMin
    extent1_YMax = extent1.YMax
    extent2_XMin = extent2.XMin
    extent2_XMax = extent2.XMax
    extent2_YMin = extent2.YMin
    extent2_YMax = extent2.YMax

    if (bool_point(extent2_XMin, extent2_YMax, extent1_XMin, extent1_YMin, extent1_XMax, extent1_YMax) == True) or (
            bool_point(extent2_XMax, extent2_YMax, extent1_XMin, extent1_YMin, extent1_XMax, extent1_YMax) == True) or (
            bool_point(extent2_XMax, extent2_YMin, extent1_XMin, extent1_YMin, extent1_XMax, extent1_YMax) == True) or (
            bool_point(extent2_XMin, extent2_YMin, extent1_XMin, extent1_YMin, extent1_XMax, extent1_YMax) == True):
        return True
    else:
        return False


def clip_method2(in_raster, in_extent, clip_shp, in_field_name, out_path, speed):
    for row in arcpy.SearchCursor(clip_shp, where_clause=expression):
    	t1 = time.time()
        mask = row.getValue("Shape")
        if bool_extents(in_extent, mask.extent):
            FID_name = str(row.getValue(in_field_name)) + '.tif'
            output_path = os.path.join(out_path + '/', FID_name)
            if os.path.exists(output_path):
                FID = str(row.getValue("FID"))
                (file_path, file_name) = os.path.split(in_raster)
                FID_name = str(file_name).replace('.tif', '') + '-' + FID + '-' + FID_name
                output_path = os.path.join(out_path + '\\', FID_name)
            if speed == "HighQuality":
                mask_raster = arcpy.sa.ExtractByMask(in_raster, mask)
                arcpy.CopyRaster_management(in_raster=mask_raster,
                                            out_rasterdataset=output_path,
                                            config_keyword="DEFAULTS",
                                            background_value=255,
                                            nodata_value=None,
                                            onebit_to_eightbit="",
                                            colormap_to_RGB="",
                                            pixel_type="8_BIT_UNSIGNED",
                                            scale_pixel_value=None,
                                            RGB_to_Colormap=None)
            if speed == "HighSpeed":
                arcpy.Clip_management(in_raster, '#', output_path, mask, None, 'ClippingGeometry', "MAINTAIN_EXTENT")
            t2 = time.time() - t1
            arcpy.AddMessage(FID_name + ' is generated successfully,  total time:' + str(round(t2, 2)))


if rasters:
    create_dir(output_Path)
    if mosaic_first == 'NO':
        for raster in rasters:
            r = arcpy.sa.Raster(raster)
            arcpy.AddMessage(raster + "--is being processed")
            clip_method2(raster, r.extent, clip_Shp, fieldName, output_Path, workSpeed)
            arcpy.AddMessage(raster + "--has done")
    else:
        rasterList = ";".join(rasters)
        r = arcpy.sa.Raster(rasters[0])
        desc = arcpy.Describe(rasters[0])
        outCoorSystem = desc.spatialReference
        dataType = desc.DataType
        cellWidth = r.meanCellWidth
        bandCount = r.bandCount
        arcpy.MosaicToNewRaster_management(rasterList, output_Path, "MosaicRaster.tif", outCoorSystem,"8_BIT_UNSIGNED", None, bandCount, "MAXIMUM", "FIRST")
        arcpy.AddMessage('raster is mosaiced successfully!')
        m = arcpy.sa.Raster(output_Path + "/MosaicRaster.tif")
        clip_method2(output_Path + "/MosaicRaster.tif", m.extent, clip_Shp, fieldName, output_Path, workSpeed)

else:
    arcpy.AddMessage('folder does not exist raster!')

其中注释部分为绝对路径,供测试时使用。
在ArcGIS中导入脚本并设置参数,如图所示:参数设置
脚本界面
输出结果

输出结果

  • 2
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值