【python+GIS】读取多波段hdf影像并将其中的波段转为tif

一、应用背景

以前数据量小的时候，用arcgis进行HDF文件转TIF的操作，但数据量一大，就很费时间，因此来找一找批处理方法。

二、参考文献

《读取HDF或者NetCDF格式的栅格数据》
《GDAL学习三----关于读取HDF数据并保存为tif数据》
这个是参考文献的读取hdf的原代码：

from osgeo import gdal

root_ds = gdal.Open('example.hdf')
# 返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息
# tuple中的第一个元素描述的是数据子集的全路径
ds_list = root_ds.GetSubDatasets()

band_1 = gdal.Open(ds_list[11][0])  # 取出第12个数据子集（MODIS反射率产品的第一个波段）
arr_bnd1 = band_1.ReadAsArray()  # 将数据集中的数据转为ndarray

# 创建输出数据集，转为GeoTIFF进行写入
out_file = 'sr_band1.tif'
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.CreateCopy(out_file, band_1)
out_ds.GetRasterBand(1).WriteArray(arr_bnd1)
out_ds.FlushCache()

# 关闭数据集
out_ds = None
root_ds = None

三、读取HDF中的一个波段，并保存

咱要修改为读一个文件夹里的所有的HDF文件，且读取一个波段的代码：

from osgeo import gdal
import os 

input_folder = r"D:\DATA\ArcGIS\sabal\1.NDVI\1.initail_data" #多个hdf文件所在文件夹路径
input_folder_list=os.listdir(input_folder) #读取文件夹里所有未见
hdf_files=list() #创建一个只装hdf格式的列表
for filename in input_folder_list:  #遍历
    if os.path.splitext(filename)[1] == '.hdf':  #不管文件名里面有多少个hdf，都只认最后一个hdf
        hdf_files.append(filename)  #将文件夹里的hdf文件加入只有hdf的列表
output_folder=r'D:\DATA\ArcGIS\sabal\1.NDVI\2.hdf_to_tif'  #设置输出文件夹

for i in range(0,len(hdf_files)): 
    root_ds = gdal.Open(input_folder+'\\'+hdf_files[i]) #打开
    ds_list = root_ds.GetSubDatasets() #读取hdf里的所有数据集

    band_1 = gdal.Open(ds_list[0][0])  # 取出第1个数据子集（MOD13产品的第一个波段NDVI）
    arr_bnd1 = band_1.ReadAsArray()  # 将数据集中的数据转为ndarray

    # 创建输出数据集，转为GeoTIFF进行写入
    out_file = output_folder+'\\'+'NDVI_'+hdf_files[i][9:23]+'.tif' #创建输出数据集的名称，[9:23]是该hdf文件的第10-23个字符
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.CreateCopy(out_file, band_1)
    out_ds.GetRasterBand(1).WriteArray(arr_bnd1)
    out_ds.FlushCache()

    # 关闭数据集
    out_ds = None
    root_ds = None

结果如下：

将转出图像和原HDF图像进行比较，投影没变，数值相等：

四、读取HDF的所有波段，并保存

相比较读取一个波段，增加了一个循环实现遍历所有波段：

from osgeo import gdal
import os 

input_folder = r"D:\DATA\ArcGIS\sabal\1.NDVI\1.initail_data" #多个hdf文件所在文件夹路径
input_folder_list=os.listdir(input_folder) #读取文件夹里所有未见
hdf_files=list() #创建一个只装hdf格式的列表
for filename in input_folder_list:  #遍历
    if os.path.splitext(filename)[1] == '.hdf':  #不管文件名里面有多少个hdf，都只认最后一个hdf
        hdf_files.append(filename)  #将文件夹里的hdf文件加入只有hdf的列表
output_folder=r'D:\DATA\ArcGIS\sabal\1.NDVI\2.hdf_to_tif'  #设置输出文件夹

for i in range(0,len(hdf_files)): 
    root_ds = gdal.Open(input_folder+'\\'+hdf_files[i]) #打开
    ds_list = root_ds.GetSubDatasets() #读取hdf里的所有数据集
    for j in range(0,len(ds_list)):
        band_1 = gdal.Open(ds_list[j][0])  # 取出第j个波段
        arr_bnd1 = band_1.ReadAsArray()  # 将数据集中的数据转为ndarray

    # 创建输出数据集，转为GeoTIFF进行写入
        out_file = output_folder+'\\'+hdf_files[i][9:23]+' - '+ds_list[j][1][12:-42]+'.tif' #创建输出数据集的名称，[9:23]是该hdf文件的第10-23个字符
        driver = gdal.GetDriverByName('GTiff')
        out_ds = driver.CreateCopy(out_file, band_1)
        out_ds.GetRasterBand(1).WriteArray(arr_bnd1)
        out_ds.FlushCache()

    # 关闭数据集
    out_ds = None
    root_ds = None

结果如下所示：