基於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:\HDFFILES\ALLHDFS' #H:\HDFFILES\ALLHDFS
outDir=r'H:\HDFFILES\ALLTIFS'

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:\HDFFILES\ALLTIFS'
outDir=r'D:\HDFFILES\mosaic'
#查找每個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:\HDFFILES\ALLTIFS'
outDir=r'H:\HDFFILES\reproj'
originReferenceFile=r'H:\HDFFILES\ALLTIFS\A2006001.h26v05.tif'
referenceFile=r'H:\HDFFILES\reference.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\reproj'#輸入文件所在文件夾
outDir=r'D:\HDFFILES\Mask'#結果文件夾
arcpy.gp.overwriteOutput=1
arcpy.env.workspace = inDir  #所有柵格影像所在文件夾
rasters = arcpy.ListRasters("*", "tif")
mask= r'D:\00WORK\05QinlingProject\00DATA\QinlingBoundary\Dissolve.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格式時最大值發生了改變，不知道原因是什么，初步判斷為重采樣方法為最近鄰采樣，但是不是很確定