ArcPy 批处理之 [ hdf转tif ]; [ Con函数 ]; 镶嵌至新栅格 [ Mosaic to New Raster ]; 重投影[ Reproject ]; 按掩膜提取[ Extract by Mask ]; [ 按条件乘积 ]; 以表格显示分区统计[ Zonal Statistics As Table ];汇总属性表
一、 ArcPy 批量将文件夹内的 *.hdf 文件转为 *.tif 文件:
#encoding:utf-8
## hdf2tif import arcpy import os inPath = r'E:\Data\S00_DataHdf\\' outPath= r'E:\Data\S01_DataTif\\' for dirpath, dirnames, filenames in os.walk(inPath): for file in filenames: if file.endswith('.hdf'): arcpy.ExtractSubDataset_management(inPath+file, outPath+os.path.splitext(file)[0][11:26]+".tif", "1") # 11:26为截取作为新文件名的字符串部分;第三个参数0,1,2为第1、2、3个波段数,此处为第二个波段 print("OK")
二、 ArcPy 批量采用Con函数对栅格文件数值区间处理:
#encoding:utf-8
## Con 去除异常值 import os import arcpy from arcpy.sa import * arcpy.env.overwriteOutput=1 arcpy.env.workspace = r'E:\Data\S01_DataTif' files = arcpy.ListRasters() outPath = r'E:\Data\S02_Data_Con\\' for file in files: if os.path.splitext(file)[1]=='.tif': out_raster=Con(file, Float(file)*0.0001,"","VALUE <= 32700") out_raster.save(outPath+file[1:5]+file[9:]) # 如原数据名称太长(>13个字符?)将导致写入失败,需要调整命名长度 print("OK")
三、 ArcPy 批量将同一时间段的栅格文件合并,即镶嵌至新栅格 (新文件的数据类型依据自身数据而定):
#encoding:utf-8 import arcpy import os import os.path arcpy.env.overwriteOutput = True arcpy.env.workspace = r'E:\Data\S02_Data_Con' output_location = r'E:\Data\S03_Data_Mosaic\\' rasters = arcpy.ListRasters("*", 'tif') Date_List=list(set(rasters)) for i in range(0,len(Date_List)): #整理Date_List,提取时间序列 Date_List[i]=Date_List[i][:4] Date_List=list(set(Date_List)) #删除重复时间序列。 for i in range(0, len(Date_List)):#按照相同时间序列判定,拼接所有同一时间序列内的tif文件。 strList='' for j in range(0, len(rasters)):#选取某一时间序列,判定符合要求的文件名添加到字符串中。 if rasters[j][:4]==Date_List[i]: strList+=rasters[j]+";" tifName=Date_List[i]+ ".tif" #设置对应时间序列的输出文件名 ,而后进行拼接。 arcpy.MosaicToNewRaster_management(strList, output_location, tifName, "#", "32 bit float", "", "1", "MAXIMUM", "MATCH") print(Date_List[i]+" OK") print("OKK")
四、 ArcPy 批量重投影[ Reproject ]:
#encoding:utf-8 """ 批量重投影: Coordinate_System参数可参照相应矢量文件的.prj文件,将双引号改为单引号即可 """ import os import arcpy from arcpy.sa import * arcpy.env.overwriteOutput=1 arcpy.env.workspace = r'E:\Data\S03_Data_Mosaic' files = arcpy.ListRasters() outPath = r'E:\Data\S04_Data_RePrj\\' for file in files: Coordinate_System = "PROJCS['Krasovsky_1940_Albers',GEOGCS['GCS_Krasovsky_1940',DATUM['D_Krasovsky_1940',SPHEROID['Krasovsky_1940',6378245.0,298.3]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',105.0],PARAMETER['standard_parallel_1',25.0],PARAMETER['standard_parallel_2',47.0],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]" cellsize = 500 arcpy.ProjectRaster_management(file, outPath+file, Coordinate_System, "BILINEAR", cellsize, "", "", "") print(file+" OK") print("OKK")
五、 ArcPy 按掩膜提取[ Extract by Mask ] :
#encoding:utf-8 """ 批量掩膜处理 """ import os import arcpy from arcpy import env from arcpy.sa import * arcpy.env.overwriteOutput=1 inPath = r'E:\Data\S04_Data_RePrj' outPath= r'E:\Data\S05_Data_Mask\\' maskfile = r'E:\Data\shp\China.shp' arcpy.env.workspace = inPath files = arcpy.ListRasters() for file in files: result = ExtractByMask(file, maskfile) result.save(outPath+file) print(file+"OK") print("OKK")
六、 ArcPy 将不同文件夹中的数据按条件相乘(此处以年份为例) :
#encoding:utf-8
import arcpy from arcpy import env from arcpy.sa import * import os import os.path arcpy.env.overwriteOutput = True extension_list = [".tif",".tiff"] output_location = r'E:\Data\S06_Data_Times\\' dir_1 = r'E:\Data1\Annual_Index' data_1 = [] for file in os.listdir(dir_1): for extension in extension_list: if file.endswith(extension): data_1.append(file) # 2005.tif data_1[i][:4] dir_2 = r'E:\Data2\Annual_Index2' extension_list = [".tif",".tiff"] data_2 = [] for file in os.listdir(dir_2): for extension in extension_list: if file.endswith(extension): data_2.append(file) # Index2005.tif data_2[i][5:9] for i in range(0, len(data_1)):# Data1 for j in range(0, len(data_2)): # Data2 OutRas=0 if data_1[i][:4]==data_2[j][5:9]: OutRas = Times(Raster(dir_1+"\\"+data_1[i]),Raster(dir_2+"\\"+data_2[j])) OutRas.save(output_location+data_2[j][:2]+"Index"+data_2[j][5:9]+ ".tif") print(data_1[i]+" ok") print("All is well")
七、 ArcPy 批量以表格显示分区统计[ Zonal Statistics As Table ],最终保存为*.dbf 和*.xls 格式:
# -*- coding: utf-8 -*- """ ZonalStatisticsAsTable "DATA":{是否忽略空值,DATA表示不计算空值} """ import arcpy from arcpy import env from arcpy.sa import * import os in_path = r'E:\Data\S06_Data_Times' out_path = r'E:\Data\S07_Data_Stats\\' shp = r'E:\Data\shp\China.shp' zoneField = "CityID" env.workspace = in_path files = arcpy.ListRasters() for file in files: outTable= out_path+file[:-4]+".dbf" outExcel= out_path+file[:-4]+".xls" outZSaT = ZonalStatisticsAsTable(shp,zoneField,file,outTable,"DATA", "MEAN") arcpy.TableToExcel_conversion(outTable, outExcel) print(file+" OK") print("OKK")
八、 ArcPy 批量汇总属性表:
#encoding:utf-8 import arcpy from arcpy import env from arcpy.sa import * import os import os.path import xlrd import xlwt arcpy.env.overwriteOutput = True dir = r'E:\Data\S07_Data_Stats' output_location = r'E:\Data\S08_Summary\\' extension_list = [".xls"] data_list = [] for file in os.listdir(dir): for extension in extension_list: if file.endswith(extension): data_list.append(file) excel_1 = xlrd.open_workbook(dir+"\\"+data_list[1]) sheet_1 = excel_1.sheet_by_index(0) cols_ = sheet_1.ncols row_max = 0 row_max_n = 0 for i in range(0, len(data_list)): excel = xlrd.open_workbook(dir+"\\"+data_list[i]) # 打开excel文件 sheet = excel.sheet_by_index(0) if row_max < sheet.nrows: row_max = sheet.nrows row_max_n = i #最多列出现的文件序号 # 创建一个workbook对象(Excel文件) workbook = xlwt.Workbook(encoding='utf-8',style_compression=0) worksheet = workbook.add_sheet('data',cell_overwrite_ok=True) # 向sheet页中添加数据:worksheet.write(行,列,值) worksheet.write(0,0,'CityID') # 第1行第1列写入数据 for j in range(0,row_max+1): worksheet.write(j+1,0,j)# 写入编号 for i in range(0, len(data_list)): excel = xlrd.open_workbook(dir+"\\"+data_list[i]) # 打开excel文件 sheet = excel.sheet_by_index(0) worksheet.write(0,i+1,data_list[i][:2]+"_"+data_list[i][5:-4]) # 第1行第1列写入数据 for j in range(1,sheet.nrows): dat_id = int(sheet.cell_value(j,1)) # 第2列为CityID dat_value = round(sheet.cell_value(j,4),3) # 第4列为Mean值 if dat_value > 0: worksheet.write(dat_id+1,i+1,dat_value) # 将以上内容保存到指定的文件中 workbook.save(output_location+"ModisNPPStat.xls") print("ok")
- END -
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
· 上周热点回顾(2.17-2.23)