进行批量掩模裁剪的准备工作:
1.将边界文件(线元素)转为掩膜文件(面元素)
2.将遥感图像和掩膜文件存储在同一文件夹下(路径一定要一致!!)
3.将需要进行操作的流域名汇总在EXEL中第一列
4.检查Arcgis的Arcgis administrator 中各个工具箱的授权情况,如果过期要重新授权
第一步:线元素转面元素
#-*- coding: UTF-8 -*- #识别中文
''' 这篇代码是在尝试将线文件通过调用FeatureToPolygon工具转化成面文件以便于后期的掩模裁剪时使用 '''import arcpy #导入模块arcpy.gp.overwriteOutput=1arcpy.env.workspace="F:\\budyko_ds\\" #输入环境outPath="F:\\budyko_ds\\ET_1\\" #输出路径files=arcpy.ListFeatureClasses() #列出工作空间中的要素类,受名称、要素类型和可选要素数据集的限制#必须先设置工作空间环境,之后才能使用多个列表函数for file in files:outFile=outPath+str(file[0:len(str(file))-8])+"polygon" #len用以计算字符长度,此处是在为输出文件命名,可以修改#调用management工具,括号中要素填写见:https://pro.arcgis.com/zh-cn/pro-app/tool-reference/data-management/feature-to-polygon.htmarcpy.FeatureToPolygon_management(file, outFile,"", "NO_ATTRIBUTES", "")print("All done,please!")
第二步就是批量流域的批量掩模裁剪:
#-*- coding: UTF-8 -*-
''' 这篇代码是用于从EXEL中读取众多流域的代码,进行裁减操作(在30年共360幅全球图上)并且输出存储 如果需要修改,特别注意lines(16、21、27、39、49、50、52、53、55、56、57) 嵌套的循环语句的逻辑是:先确定流域(EXEL中提取),再找他的面元素文件(.shp),最后针对每年每月的图像进行裁减并储存 '''
import arcpy
import glob
import os
import xlrdarcpy.CheckOutExtension("spatial")arcpy.gp.overwriteOutput=1arcpy.env.workspace = "F:/budyko_ds/P/"rasters = arcpy.ListRasters("*", "tif")#指定掩膜文件路径
path_file='F:/budyko_ds/P/'#存储路径中所有文件,以ls来表示,常用file
ls = os.listdir(path_file)
print(len(ls))data = xlrd.open_workbook('C:/Users/DELL/Desktop/name.xlsx') # 打开xls文件
table = data.sheets()[0] # 打开第一张表
nrows = table.nrows # 获取表的行数
#datalist用来存放数据
datalist_NAME=[]
#将table中第一列的数据读取并添加到data_list中
datalist_NAME.extend(table.col_values(0))for num in range(1,nrows): #遍历ls中存储的所有文件for i in ls:#模糊查找和datalist_NAME中存储的流域名相似的文件,作为掩膜文件if i.find(eval(datalist_NAME[num])+"_Boundary_polygon.shp")!=-1: #检查文件拓展名,读取.shp文件(注意,不是xml文件!!!)if os.path.splitext(i)[1] == ".shp":# if i.endswith(".shp"):mask = iprint maskfor raster in rasters:#判断存储最终图像的文件夹是否存在,是则存储,否则创建并存储;注意if和else后一定要加冒号,格式要对齐if os.path.exists(r"F:/Budyko_ds/P_extraction/"+eval(datalist_NAME[num])):out = (r"F:/Budyko_ds/P_extraction/"+eval(datalist_NAME[num])+"/clip_"+raster)else:os.makedirs(r"F:/Budyko_ds/P_extraction/"+eval(datalist_NAME[num]))out = "F:/Budyko_ds/P_extraction/"+eval(datalist_NAME[num])+"/clip_"+raster# 调用裁减工具arcpy.gp.ExtractByMask_sa(raster, mask, out)print("clip_"+raster+" has done")print(datalist_NAME[num]+"has done")print("All done")
经自己努力而攻克一个问题的感觉,畅快!!