ArcPy 批处理之 [ hdf转tif ]; [ Con函数 ]; 镶嵌至新栅格 [ Mosaic to New Raster ]; 重投影[ Reproject ]; 按掩膜提取[ Extract by Mask ]; [ 按条件乘积 ]; 以表格显示分区统计[ Zonal Statistics As Table ];汇总属性表

发布时间 2023-04-17 22:33:12作者: GEOZHO

一、 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 -