如何用h5py存储带GIS可识别地理空间信息的PRISM栅格为HDF5？

核心原因 #

GIS软件（如ArcGIS、QGIS，底层依赖GDAL）识别HDF5栅格的空间信息，需要遵循GDAL兼容的元数据规范或CF（Climate and Forecast）元数据标准，而非仅简单存储坐标数组。MODIS HDF4的元数据结构和HDF5的解析逻辑不匹配，直接照搬会导致GIS无法识别空间参考。

具体解决方案 #

1. 存储CF标准元数据 #

CF标准是GIS和遥感领域通用的HDF5元数据规范，GDAL、RasterIO均能识别。需为栅格数据集添加以下关键属性：

• crs：存储WKT格式的坐标参考字符串
• geotransform：存储GDAL格式的仿射变换六元组（左上角x、x分辨率、x旋转、左上角y、y旋转、-y分辨率，PRISM栅格y轴从上到下递减，y分辨率为负）
• 为维度标记y（行）和x（列）标签

2. 标准HDF5结构组织示例 #

将每个PRISM栅格作为独立数据集，在数据集层级绑定空间元数据，示例代码如下：

import h5py
import rasterio
from rasterio.crs import CRS

# 读取单张PRISM栅格
prism_tif = "prism_tmin_20240101.tif"
with rasterio.open(prism_tif) as src:
    raster_data = src.read(1)
    crs_wkt = src.crs.to_wkt()
    gdal_gt = src.transform.to_gdal()  # 获取GDAL兼容的仿射变换参数
    rows, cols = src.shape
    nodata_val = src.nodata

# 创建并写入HDF5文件
with h5py.File("prism_h5_collection.h5", "w") as h5_file:
    # 创建栅格数据集
    ds = h5_file.create_dataset("tmin_20240101", data=raster_data, dtype=raster_data.dtype)
    
    # 绑定空间元数据
    ds.attrs["crs"] = crs_wkt
    ds.attrs["geotransform"] = gdal_gt
    ds.attrs["_FillValue"] = nodata_val  # 标记无效值
    ds.attrs["units"] = "degrees Celsius"  # 可选：添加数据单位
    
    # 标记维度含义
    ds.dims[0].label = "y"
    ds.dims[1].label = "x"
    # 可选：存储坐标数组
    h5_file.create_dataset("y_coords", data=src.y)
    h5_file.create_dataset("x_coords", data=src.x)

3. 适配GDAL的额外优化 #

若需更贴合GDAL解析逻辑，可在HDF5根节点添加GDAL专属元数据：

# 在根节点添加子数据集描述（多栅格场景）
h5_file.attrs["GDAL_METADATA"] = (
    f"SUBDATASET_1_NAME=HDF5:{h5_file.filename}://tmin_20240101\n"
    "SUBDATASET_1_DESC=PRISM Daily Minimum Temperature"
)

4. 验证方法 #

• 用GDAL命令行验证：执行gdalinfo HDF5:"prism_h5_collection.h5"://tmin_20240101，若输出包含CRS和地理变换信息，说明GIS可识别。
• 直接在QGIS中导入HDF5文件，选择对应数据集，查看空间位置是否正确加载。

关键注意事项 #

• 不要直接照搬HDF4的元数据结构，HDF5与HDF4的属性存储逻辑、GIS解析规则存在差异。
• 优先使用WKT格式存储CRS，兼容性优于单独的EPSG代码。
• 确保仿射变换参数的符号正确，PRISM栅格y轴从上到下数值递减，y分辨率需设为负数。

内容的提问来源于stack exchange，提问作者Nick C