基于Python批量将HDF文件转为GeoTiff格式并进行拼接、投影转换和矢量裁剪

问题：通常使用的是MODIS的HDF格式数据，但是在近期使用GLASS产品(梁顺林团队)产品时，发现该HDF文件被GDAL库中函数读取后并无子数据集等信息，即包括的波段等等，中间有点麻烦，这篇文章主要记录基于Arcpy对于hdf产品批量的一条龙操作：

格式转换为TIF、拼接、投影转换、裁剪

注意：记录均设置为英文路径，以及不要以数字开头

ArcGIS也可以打开HDF影像

基于Python-GDAL将HDF文件转为GeoTiff格式

这里可能需要先手工转好一下以获取参数

https://www.jianshu.com/p/cd1e1080bd2b

基于Arcpy将HDF文件转为GeoTiff格式

import arcpy

arcpy.CheckOutExtension("Spatial")
arcpy.gp.overwriteOutput=1

inDir=r'H:HDFFILESALLHDFS' #H:HDFFILESALLHDFS
outDir=r'H:HDFFILESALLTIFS'

arcpy.env.workspace=inDir
rasters=arcpy.ListRasters("*","hdf")

for raster in rasters:
    print(raster)
    outName=outDir+'\'+raster[15:30]+'.tif'
    arcpy.ExtractSubDataset_management(raster,outName)
    print(outName)

基于Arcpy将Tiff文件进行批量拼接

import os
import arcpy


inDir=r'D:HDFFILESALLTIFS'
outDir=r'D:HDFFILESmosaic'
#查找每个tile所对应的文件序列
for year in range(2006,2007):#对2006年数据进行拼接，主要也是找出天数命名规律进行文件名设置
    print(year)  
    for eday in range(1,362,8):
        filename1=inDir+'\'+'A'+str(year)+str(eday).rjust(3,'0')+'.h27v05.tif'#影像tile1
        filename2=inDir+'\'+'A'+str(year)+str(eday).rjust(3,'0')+'.h26v05.tif'#影像tile2
        expression=filename1+';'+filename2
        outName=str(year)+str(eday).rjust(3,'0')+'_mosaic.tif'
        print(expression)
        arcpy.MosaicToNewRaster_management(expression,outDir,outName,"#","#", "#", "1", "#","#") 
        print(outName)

基于Arcpy将GeoTiff格式文件进行批量投影转换

先手工导出一个reference tif文件，以建立一个投影转换关系，再应用到其他文件上

import arcpy

inDir=r'H:HDFFILESALLTIFS'
outDir=r'H:HDFFILES
eproj'
originReferenceFile=r'H:HDFFILESALLTIFSA2006001.h26v05.tif'
referenceFile=r'H:HDFFILES
eference.tif'

arcpy.CheckOutExtension("spatial")
arcpy.gp.overwriteOutput=1

arcpy.env.workspace=inDir
rasters=arcpy.ListRasters("*","tif")

for raster in rasters:
    print(raster)
    outName=outDir+'\'+raster[0:15]+'_Reproject.tif'
    print(outName)
    arcpy.ProjectRaster_management(raster, outName,referenceFile, "#",
                                   '#','#','#',originReferenceFile)

基于Arcpy批量裁剪GeoTiff格式数据

需要先准备好矢量边界shp数据

import arcpy
arcpy.CheckOutExtension("spatial")
inDir=r'D:HDFFILES
eproj'#输入文件所在文件夹
outDir=r'D:HDFFILESMask'#结果文件夹
arcpy.gp.overwriteOutput=1
arcpy.env.workspace = inDir  #所有栅格影像所在文件夹
rasters = arcpy.ListRasters("*", "tif")
mask= r'D:0WORK5QinlingProject0DATAQinlingBoundaryDissolve.shp'  #用于提取的矢量掩膜
for raster in rasters:
    print(raster)
    outName= outDir+'\'+raster[0:7]+'.tif'#命名可适当更改
    print(outName)
    arcpy.gp.ExtractByMask_sa(raster, mask, outName)
print("OK")

经过以上步骤的处理后算是得到了研究区内的长时间序列影像，但是还是有一个问题，就是在HDF文件转为TIF格式时最大值发生了改变，不知道原因是什么，初步判断为重采样方法为最近邻采样，但是不是很确定