ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

HDF转TIFF实战:GDAL与ENVI完整转换指南

HDF转TIFF实战:GDAL与ENVI完整转换指南 简介针对遥感与地理数据处理中常见的HDF转TIFF需求这份资源提供可直接使用的Python转换脚本面向需要将HDF文件导入ENVI、GIS或其他仅支持TIFF格式软件的开发者、科研人员与学生。脚本基于Python编写可读取HDF数据集并写出带空间参考的GeoTIFF适合处理单波段或多波段栅格数据便于后续在ENVI等平台中开展分析。压缩包内共1个文件为py格式的源码脚本整体仅1KB代码结构精简适合快速查看与二次修改。目前已有2285人学习下载。读者拿到脚本后可结合自身HDF文件中的数据集名称、边界范围及坐标系信息调整参数实现格式转换脚本同时展示了h5py与rasterio两个关键库的配合用法对理解HDF内部结构、栅格空间参考写入及Python批量转换流程也有直接参考价值。1. 拿到 hdf 文件先别急着双击打开转 tif 前先弄懂这件事做遥感数据处理的工程师几乎都经历过这样的场景从 NASA 或 LP DAAC 下载回来的 MODIS 产品是一堆.hdf文件文件名里带了一长串产品和日期标识例如MOD09GA.A2024185.h11v04.006.2024190072813.hdf。双击打开普通看图软件不认ArcGIS 直接拖进去也经常只看到一条黑带或者根本没有渲染。用户拿到的标题里同时出现了envi和Python说明实际需求是从这份 HDF 里提取出能被 GIS 和深度学习框架直接读取的 GeoTIFF而且希望这条路既能靠 ENVI 的可视化解决也能用 Python 批量自动化。这个需求的难点不在“转换”动作本身而在于 HDF 是容器格式一个文件里塞着多个科学数据集子数据集投影、缩放因子、填充值都得先读出来再写进 GeoTIFF 的元数据里。这篇文章就从文件结构讲起把 hdf 转 tif 的完整路径拆开ENVI 和 Python 两条线都会覆盖到。2. hdf 到 tif 的本质先拆 HDF-EOS 子数据集再谈投影与比例因子2.1 为什么不能像打开 tif 一样直接打开 hdf这里说的 hdf 是 Hierarchical Data Format和华为鸿蒙系统里那个 HDF 硬件驱动框架完全是两回事别混在一起。遥感领域最常见的是 HDF4 和基于它的 HDF-EOS 扩展格式MODIS、ASTER、部分 VIIRS 产品都走这条路。一个 HDF 文件内部采用树状结构组织根节点下面挂着多个科学数据集每个数据集单独存储了一层栅格数据还带自己的属性信息比如scale_factor、add_offset、_FillValue、long_name等。标准 GeoTIFF 是单栅格文件即使有多波段也是统一存储在一个文件里。所以 hdf 转 tif 并不是拿个文件转换工具套一下就完事真正的动作是先列出 HDF 里的子数据集清单选定需要的波段把它的投影信息和地理变换参数一并读取出来再写入一个新的 GeoTIFF 文件。如果直接对整个 HDF 文件做转换GDAL 会返回一大串子数据集路径而不是一个可直接使用的栅格。2.2 用 gdalinfo 查看 hdf 的子数据集清单GDAL 自带gdalinfo命令这是排查 HDF 结构最先要用的工具。安装 GDAL 的方式很多Windows 上可以用 OSGeo4W 或 condaLinux 直接用apt install gdal-bin或者conda install -c conda-forge gdal。在终端里对目标文件执行gdalinfo MOD09GA.A2024185.h11v04.006.2024190072813.hdf输出内容很长核心在Subdatasets:这一段。典型输出如下Subdatasets: SUBDATASET_1_NAMEHDF4_EOS:EOS_GRID:MOD09GA.A2024185.h11v04.006.2024190072813.hdf:MODIS_Grid_1km_2D:sur_refl_b01_1 SUBDATASET_1_DESC[200x164] sur_refl_b01_1 MODIS_Grid_1km_2D (16-bit unsigned integer) SUBDATASET_2_NAMEHDF4_EOS:EOS_GRID:MOD09GA.A2024185.h11v04.006.2024190072813.hdf:MODIS_Grid_1km_2D:sur_refl_b02_1 SUBDATASET_2_DESC[200x164] sur_refl_b02_1 MODIS_Grid_1km_2D (16-bit unsigned integer)SUBDATASET_1_NAME里的完整字符串就是后续 Python 或gdal_translate真正要打开的数据集标识。[200x164]是栅格的行列数16-bit unsigned integer是存储类型。这段信息直接决定了后续怎么设计转换脚本。对于想要一次看清所有波段的情况可以只过滤子数据集列表gdalinfo MOD09GA.A2024185.h11v04.006.2024190072813.hdf | grep -E SUBDATASET_.*NAME这里grep的正则匹配了所有子数据集名方便拷贝进脚本或做人工核对。注意不同 MODIS 产品的子数据集命名规律不一样MOD13Q1 是MODIS_Grid_16DAY_250m_500m_VIMOD09GA 是MODIS_Grid_1km_2D写代码时最好不要写死路径前缀。2.3 读懂 MODIS 产品的行列/投影/缩放因子为 hdf转tif 做准备拿到子数据集名称后还需要读它的投影和属性。对单个子数据集执行gdalinfo HDF4_EOS:EOS_GRID:MOD09GA.A2024185.h11v04.006.2024190072813.hdf:MODIS_Grid_1km_2D:sur_refl_b01_1输出会包含关键的地理参考信息。MODIS 陆地产品的标准投影是正弦曲线投影SinusoidalGDAL 通常将其识别为SIN或ESRI:54008同时会给出Origin (...),Pixel Size (...)这两行参数。它们决定了转换后 GeoTIFF 的仿射变换矩阵。同样的命令能读到波段属性例如Band 1 Block200x164 TypeInt16, ColorInterpUndefined Min0.000 Max8000.000 NoData Value-28672 Scale Factor0.000100 Offset0.000000Scale Factor和Offset是 hdf 转 tif 时最容易丢的信息。MODIS 反射率产品内部一个像元值可能是 8000但真实反射率是乘上 0.0001 后的 0.8。很多初学者直接把原始整数写进 tif后续做植被指数计算时结果完全错乱。提前跑一次gdalinfo把这里的信息摸清楚后面写 Python 脚本时才知道要不要乘缩放因子。3. 用 Python 批量 hdf转tifGDAL API 与栅格读写细节3.1 安装 GDAL 依赖与验证子数据集可读Python 环境下做栅格转换最可靠的方式是 GDAL 的 Python 绑定常用导入名是from osgeo import gdal。安装时用 conda 不容易出现 DLL 缺失问题conda install -c conda-forge gdal如果项目用的是 venv也可以用pip install gdal但 Windows 上 pip 版本经常需要和已安装的 GDAL 二进制版本严格对应报错率比 conda 高不少因此我一般会建议优先 conda 环境。安装完成后先跑一个最小验证确认能打开 HDF 并列出子数据集from osgeo import gdal file_path MOD09GA.A2024185.h11v04.006.2024190072813.hdf ds gdal.Open(file_path) for i in range(ds.RasterCount): sub ds.GetSubDatasets()[i] print(sub)ds.GetSubDatasets()返回一个列表每个元素是(subdataset_path, description)的元组。如果打印结果为 0 条说明当前 GDAL 版本缺少 HDF4 驱动需要额外安装libgdal-hdf4或者重新编译启用 HDF4 支持的 GDAL。这一步不通过后面所有脚本都无从谈起。3.2 单文件 hdf转tif 的标准 Python 脚本下面的脚本完成一个基本但完整的工作读取指定的 HDF 子数据集把投影和地理变换信息一并复制到输出的 GeoTIFF 中。from osgeo import gdal def hdf_sub_to_tif(hdf_path, subdataset_name, out_tif): gdal.UseExceptions() # 直接用子数据集完整路径打开 src_ds gdal.Open(subdataset_name) if src_ds is None: raise RuntimeError(f无法打开子数据集: {subdataset_name}) # 获取仿射变换参数和投影 geotransform src_ds.GetGeoTransform() projection src_ds.GetProjection() if geotransform is None or projection is None: print(警告: 当前子数据集缺少地理参考信息) # 创建输出文件复制波段数、行列数、数据类型 driver gdal.GetDriverByName(GTiff) dst_ds driver.Create( out_tif, src_ds.RasterXSize, src_ds.RasterYSize, src_ds.RasterCount, gdal.GDT_Float32, options[COMPRESSLZW, TILEDYES] ) dst_ds.SetGeoTransform(geotransform) dst_ds.SetProjection(projection) # 逐波段复制数据随后清理 for band_idx in range(1, src_ds.RasterCount 1): src_band src_ds.GetRasterBand(band_idx) data src_band.ReadAsArray() dst_band dst_ds.GetRasterBand(band_idx) dst_band.WriteArray(data) src_ds None dst_ds None hdf_path MOD09GA.A2024185.h11v04.006.2024190072813.hdf subdataset_path HDF4_EOS:EOS_GRID:MOD09GA.A2024185.h11v04.006.2024190072813.hdf:MODIS_Grid_1km_2D:sur_refl_b01_1 hdf_sub_to_tif(hdf_path, subdataset_path, sur_refl_b01_1.tif)脚本的逻辑分四步。先通过gdal.Open打开某个子数据集因为 HDF 作为一个容器无法直接转换为单文件 GeoTIFF必须定位到内部某一个科学数据集。接着用GetGeoTransform()和GetProjection()把原始位置信息取出来这两行决定了输出 tif 在 GIS 软件里能不能落到正确的地理位置上。driver.Create中的COMPRESSLZW是一种无损压缩选项对遥感影像效果好TILEDYES则让 tif 按块存储后续读取局部窗口时性能更好。gdal.GDT_Float32意味着输出为 32 位浮点能容纳小数标度避免整数化导致精度丢失。3.3 处理填充值、比例因子与 NoData 的三个必调参数实际项目中直接把原始波段值写进 tif 往往不够还必须处理以下三个参数。第一是填充值。MODIS 产品的无效像元通常用-28672之类的大负数表示直接参与计算会污染结果。在写入前应读取子数据集属性里的 NoData 值并同步设置到输出 tif 的波段上src_band src_ds.GetRasterBand(1) nodata src_band.GetNoDataValue() # 处理完数据后在输出波段上设置 if nodata is not None: dst_band.SetNoDataValue(nodata)第二是比例因子。前面提到 MODIS 反射率数据存的是扩大 10000 倍的整数想要真实反射率就逐像元乘scale_factor。在脚本里做这个乘法时要使用data.astype(np.float32) * scale避免整型乘法溢出。第三是在写入时保持数据类型一致性。如果输出选GDT_Float32而源数据是Int16WriteArray时 GDAL 会自动做类型转换反而是数据结构复杂时手动把src_ds.GetRasterBand(...).DataType读出来传给Create更不容易丢信息。参数取值可以按实际需求调整如果输出仅用于深度学习分割不需要保留 NoData可直接用源数据ReadAsArray()后的数值范围做 0-255 拉伸此时COMPRESS可以换为DEFLATE以提升压缩比如果输出要进 ArcGIS 做分析则必须保留 NoData。三个参数的取舍顺序是先确认存什么值再定数据类型最后调整压缩方式。4. ENVI 下的 hdf转tif可视化检查与批量导出需要注意的坑4.1 ENVI 打开 hdf 的两种常见方式很多从业者拿到 HDF 第一反应是打开 ENVI因为 ENVI 对 HDF-EOS 的支持比较成熟。常规做法是直接File Open As在文件类型列表里选择EOS或HDF4然后定位到目标文件。ENVI 会列出文件内的所有科学数据集勾选需要的波段就能打开到视窗中。第二种方式更推荐先拖拽.hdf文件到 ENVI 的图层窗口如果 ENVI 无法自动识别再手动指定Open As EOS HDF4打开。这里有个常见的坑是 ENVI 版本差异较新的 ENVI 默认用栅格管理器方式解析而旧版需要额外的 HDF 补丁否则打开后只显示一个空栅格或警告“未找到有效的 HDF 数据集”。ENVI 打开 HDF 后的一个明显优势是能直接叠加在已知底图上快速判断投影是否正确。MODIS 产品默认的 Sinusoidal 投影在 ENVI 里如果显示异常通常是因为 ENVI 缺少对应的投影定义文件此时需要在投影选择对话框中手工指定Sinusoidal (Sphere)或自定义ESRI:54008投影参数。4.2 在 ENVI 中叠加投影与导出 GeoTIFF当 HDF 在 ENVI 中正常显示后导出 GeoTIFF 的操作是File Export Export to Raster TIFF。在导出面板中ENVI 会询问是否保留原始投影这里务必选择Apply Map Info或保持默认的地理参考选项否则生成的是不带坐标的普通 tif后续在 ArcGIS 里打开就是一张没有位置的图。导出前建议先做一次坐标系检查右键图层名选择Edit Map Info确认Projection、X/Y coordinate和Pixel Size和原始 HDF 的gdalinfo输出一致。这一步在 ENVI 里看似多余却能在源头避免投影丢失。如果需要在导出时同时处理比例因子ENVI 的做法是在打开 HDF 时进入Data Manager选中波段后右键查看属性找到Scale Factor字段。部分 ENVI 版本不会自动应用该值需要手动对波段做Band Math乘法运算例如写公式float(b1) * 0.0001。导出的 tif 再做后续分析时数值就已经是真实反射率了。4.3 用 ENVI 做批量 hdf转tif 的边界在哪里ENVI 的图形界面适合单文件或少量文件的转换因为可视化检查方便。但当文件数量达到几十上百时逐个打开再导出的效率很低此时有两个替代方案。方案一是 ENVI 的 Batch 模式通过File Batch命令配合 IDL 脚本调用ENVI_Open_File和ENVI_Output_To_File完成批量导出。IDL 脚本的可控性比手工操作强但前提是熟悉 IDL 语法且 ENVI 的批处理模块需要单独的许可。方案二是回到 Python 路线ENVI 负责抽检和验证Python 负责批处理。实际项目里我一般会先用 ENVI 打开两三个文件目视检查投影和范围再把完整的文件列表交给 Python 脚本统一转 tif最后再用 ENVI 抽查输出结果。ENVI 在 hdf转tif 这个场景里扮演的角色是质检工具而不是批量生产工具把批量压缩的活交给 GDAL效率会明显提升。5. 从单文件到批处理hdf转tif 的自动化脚本与结果验证5.1 批量转换脚本与内存控制处理一个月的 MODIS 产品文件数量通常是几十个。最简单的批量思路是遍历目录下的所有.hdf文件每个文件读取指定子数据集并输出为 tif。但无节制的循环会带来内存问题尤其是ReadAsArray()读取大影像时单波段数据一旦超过内存上限进程就会退出。常见做法是按文件为单位循环并在每个文件处理完后用None释放对象引用。import glob import os from osgeo import gdal gdal.UseExceptions() def convert_all_hdfs(input_dir, output_dir, subdataset_idx0): os.makedirs(output_dir, exist_okTrue) hdf_files glob.glob(os.path.join(input_dir, *.hdf)) for hdf_path in hdf_files: ds gdal.Open(hdf_path) subdatasets ds.GetSubDatasets() if subdataset_idx len(subdatasets): print(f{hdf_path} 中没有第 {subdataset_idx} 个子数据集) ds None continue sub_path, desc subdatasets[subdataset_idx] # 用 basename 加编号命名避免重名 out_name os.path.splitext(os.path.basename(hdf_path))[0] f_band{subdataset_idx}.tif out_path os.path.join(output_dir, out_name) # 用 gdal.Warp 实现重投影式转换内存占用比逐波段读写更低 warp_options gdal.WarpOptions( formatGTiff, creationOptions[COMPRESSLZW, BIGTIFFIF_SAFER], resampleAlgnear ) gdal.Warp(out_path, sub_path, optionswarp_options) ds None print(f完成: {out_path})脚本里用了gdal.Warp而不是手写波段循环原因在于Warp在底层做了数据流切块能把内存峰值压到较低水平。处理超大 HDF 影像时会发现单波段ReadAsArray()直接占掉几个 GB 内存而Warp分块读写后内存占用能降低一半以上。BIGTIFFIF_SAFER表示当输出文件超过 4GB 时自动升级为 BigTIFF 格式避免大文件写成后无法打开。resampleAlgnear保留原始像元值不做插值适合分类和反射率这种要保持原始采样语义的数据。如果子数据集序号固定subdataset_idx传 0 或 1 即可需要一次转多个波段时改为双层循环就能扩展。5.2 验证输出 tif 的投影、范围和像元值批量转换完成后不能只扫一眼文件名就收工一个可靠的验证方式是回到gdalinfo检查输出文件的关键字段gdalinfo sur_refl_b01_1.tif | grep -E ^(Size|Origin|Pixel Size|Coordinate System)能同时看到Size is 200, 164、Origin、Pixel Size和Coordinate System输出说明投影信息写入成功。更进一步可以对比转换前后波段的统计值来验证像元值没有被意外改动from osgeo import gdal import numpy as np def validate_tif(tif_path, expected_nodataNone): ds gdal.Open(tif_path) band ds.GetRasterBand(1) data band.ReadAsArray() print(fshape: {data.shape}) print(fmin: {data.min()}, max: {data.max()}) if expected_nodata is not None: print(fnodata 数量: {(data expected_nodata).sum()}) ds None这里的min/max应和源 HDF 子数据集里的波段范围对得上。如果差了一个量级多半是比例因子没乘如果完全一致但缺少 NoData 设置则检查源文件_FillValue是否在转换过程中丢失。验证手段虽然基础但能在问题扩散到下游分析之前暴露出来。实际项目中把验证脚本挂在整个转换流程的末尾循环跑一遍通常五秒钟就能定位到出问题的文件。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进