python实现hdf/h5文件转tif（可批量）

所用数据：AMSR_U2_L3_DailySnow_B02_20121230.he5 雪水当量数据产品

环境：python2.7  并安装h5py模块

H5数据介绍：

HDF(Hierarchical Data Format)，设计用于存储和组织大量数据的文件格式。

h5文件中有两个核心的概念：组“group”和数据集“dataset”。组可以理解为文件夹，数据集理解为文件，文件夹中可以递归地包含文件夹、文件等。

一、查看数据信息

用python查看没搞明白，我是用matlab看的

二、转tif

查阅了一些文章和自己摸索，发现以下两种方法都可以

方法一：

我用之前nc转tif的代码改写的

先看函数

def dayextract_hdf(hdf_file, output_dir):
    fp = h5py.File(hdf_file, 'r')


    Lat = fp['HDFEOS']['GRIDS']['Northern Hemisphere']["YDim"][:]
    Lon = fp['HDFEOS']['GRIDS']['Northern Hemisphere']["XDim"][:]


    # 影像的左上角和右下角坐标
    # LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]

    # 分辨率计算
    N_Lat = len(Lat)
    N_Lon = len(Lon)
    Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
    cols = N_Lon
    rows = N_Lat

    # 设置影像的显示范围
    # -Lat_Res一定要是-的
    # geo = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
    geo = (-9036843.07384, 25067.52586, 0, 9036843.07384, 0, -25067.52586)

    # 构造projection
    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(3408)  #查看文件投影代码为3408
    src_srs_wkt = src_srs.ExportToWkt()   # 给新建图层赋予投影信息
    dayarray = fp['HDFEOS']['GRIDS']['Northern Hemisphere']['Data Fields']['SWE_NorthernDaily'][:]
    # print(len(dayarray))

    out_file =os.path.join(output_dir, hdf_file[11:-4] + ".tif")
    print(out_file)

    arcpy.env.overwriteOutput = 1  # 输出文件夹里面已经有内容的，就覆盖掉
    arcpy.CheckOutExtension("Spatial")


    # 生成tif，写入
    driver = gdal.GetDriverByName("Gtiff")

    outdataset = driver.Create(out_file, cols, rows, 1, gdal.GDT_Float32)
    outdataset.SetGeoTransform(geo)
    outdataset.SetProjection(src_srs_wkt)

    band = outdataset.GetRasterBand(1)
    band.WriteArray(data_array)
    fp.close()

之后再批量操作

    # 读取所有HDF数据
    hdf_file = r"D:\test\nc"
    tiff_dir = r"D:\test\tif"
    hdfFileList = getFileName(hdf_file, {'.he5'})
    for hdfFile in hdfFileList:
    # nc数据转换为tif数据
        print(hdfFile)
        hdf_to_tif(hdfFile, tiff_dir)

这里如果xy或是经纬度读取有问题的话，可以直接把h5文件拖进arcmap点击图层属性查看图像的范围坐标，可以手动输入geo。同理，rows和cols也可以在arcmap中查看直接赋值

方法二：

#-*- coding: UTF-8 -*-                                                                  #

import arcpy
import glob
import os

import h5py
import xlrd

arcpy.env.overwriteOutput = 1                  #输出文件夹里面已经有内容的，就覆盖掉
arcpy.CheckOutExtension("Spatial")             #检查ArcGIS扩展模块是否可用
inPath = 'D:\\test\\nc\\'
outDir = 'D:\\test\\tiff\\'

def h5_to_tiff(inPath, outDir):
	file = os.listdir(inPath)  # 读取路径中所有文件，以ls来表示，常用file
	for i in file:  # 遍历ls中的文件
		arcpy.env.workspace = inPath + i
		# 脚本中最为常用的环境变量设置就是arcpy.env.workspace，该变量用于定义当前脚本的工作目录（或者称为工作空间）
		# print(arcpy.env.workspace)
		arcpy.env.scratchWorkspace = inPath + i  # （为什么定义临时空间变量作用暂不明）
		hdfList = arcpy.ListRasters('*', 'he5')  # 按名称和栅格类型返回工作空间中的栅格列表(遍历指定类型的文件)，需先设置工作空间环境

		# 判断存储最终图像的文件夹是否存在，是则存储，否则创建并存储；注意if和else后一定要加冒号，格式要对齐
		if os.path.exists(outDir + i[:-4]):
			for hdf in hdfList:
				# 此Hdf文件有两个波段，全部输出，相应的ExtractSubDataset_management 中也要提取对应波段
				for number in range(0, 4, 2):  # 有几个波段就写几个,我这里需要输出第1个和第3个波段数据
					outPath = outDir + str(i[:-4]) + '\\'
					# 根据对 HDF 数据集创建新栅格数据集,语法 ExtractSubDataset(in_raster, out_raster, {subdataset_index})
					out = arcpy.ExtractSubDataset_management(hdf, outPath + hdf[0:-4] + str(number) + ".tif", number)
	else:
		os.makedirs(outDir + i[:-4])
		for hdf in hdfList:
			for number in range(0, 4, 2):
				outPath = outDir + str(i[:-4]) + '\\'
				out = arcpy.ExtractSubDataset_management(hdf, outPath + hdf[0:-4] + str(number) + ".tif", number)
	print(outPath)

h5_to_tiff(inPath, outDir)

这段代码参考：

https://blog.csdn.net/qq_38882446/article/details/103461759

来源：是一个橙子呀