ArcPy 批处理之 [ hdf转tif ]； [ Con函数 ]；镶嵌至新栅格 [ Mosaic to New Raster ]；重投影[ Reproject ]；按掩膜提取[ Ex

标签：Statistics 掩膜 Data list data 栅格 file import arcpy

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

标签：Statistics,掩膜,Data,list,data,栅格,file,import,arcpy
From： https://www.cnblogs.com/geozho/p/17327789.html

ArcPy 批处理之 [ hdf转tif ]； [ Con函数 ]；镶嵌至新栅格 [ Mosaic to New Raster ]；重投影[ Reproject ]；按掩膜提取[ Ex

相关文章

赞助商

阅读排行

ArcPy 批处理之 [ hdf转tif ]； [ Con函数 ]； 镶嵌至新栅格 [ Mosaic to New Raster ]； 重投影[ Reproject ]； 按掩膜提取[ Ex

相关文章

赞助商

阅读排行

ArcPy 批处理之 [ hdf转tif ]； [ Con函数 ]；镶嵌至新栅格 [ Mosaic to New Raster ]；重投影[ Reproject ]；按掩膜提取[ Ex