【Python实例】hdf文件简介及基于Python导入hdf文件

【Python实例】hdf文件简介及基于Python导入hdf文件

• .hdf文件概述
• hdf工具-HDFView
• 基于Python导入hdf文件
• • 实例1：判断是否为hdf文件
  • 实例2：打开并读取hdf4文件（地面反照率数据）
  • 实例3：绘制图形（地面反照率数据）
• 参考

HDF（Hierarchical Data Format）是一种用于存储和组织大量数据的文件格式。HDF 文件可用于多种应用，特别是在科学和工程领域。

.hdf文件概述

一个HDF5文件是一种存放两类对象的容器：dataset和group。 Dataset是类似于数组的数据集，而group是类似文件夹一样的容器，存放dataset和其他group。
在使用h5py的时候需要牢记一句话：groups类比词典，dataset类比Numpy中的数组。

• HDF4 文件通常以 .hdf 或 .h4 结尾，而 HDF5 文件通常以 .h5 或 .hdf5 结尾。
  需要专门的软件（如HDF View）才能打开预览文件的内容。HDF5 文件结构中有 2 primary objects: Groups 和 Datasets。

• 每个 dataset 可以分成两部分: 原始数据 (raw) data values 和 元数据 metadata (a set of data that describes and gives information about other data => raw data)。
  对于每一个dataset 而言，除了数据本身之外，这个数据集还会有很多的属性 attribute。在hdf5中，还同时支持存储数据集对应的属性信息，所有的属性信息的集合就叫做metadata.

hdf工具-HDFView

HDFView的下载及安装可参见另一博客-【hdf文件工具】HDFView安装及使用教程。

基于Python导入hdf文件

常用工具和库

• h5py：Python 中用于读取和写入 HDF5 文件的库。
• PyTables：另一个 Python 库，专注于高性能的 HDF5 文件操作。

实例1：判断是否为hdf文件

具体Python代码如下：

import os

def is_hdf5(file_path):
    """检查文件是否为 HDF5 格式"""
    with open(file_path, 'rb') as f:
        signature = f.read(8)
    return signature == b'\x89HDF\r\n\x1a\n'

def is_hdf4(file_path):
    """检查文件是否为 HDF4 格式"""
    with open(file_path, 'rb') as f:
        # HDF4 文件的开头通常包含特定的标志
        signature = f.read(4)
    return signature.startswith(b'HDF')

def is_hdf(file_path):
    """检查文件是否为 HDF 文件（HDF4 或 HDF5）"""
    return is_hdf4(file_path) or is_hdf5(file_path)

hdf_file_path = 'D:/path/to/your/file.hdf'  # 确保路径正确

if not os.path.exists(hdf_file_path):
    print("文件不存在，请检查路径。")
else:
    try:
        if is_hdf(hdf_file_path):
            if is_hdf5(hdf_file_path):
                import h5py
                with h5py.File(hdf_file_path, 'r') as hdf:
                    print("打开 HDF5 文件成功！")
                    print("Keys in the HDF file:", list(hdf.keys()))
            elif is_hdf4(hdf_file_path):
                from pyhdf import HDF
                hdf = HDF(hdf_file_path, HDF.READ)
                print("打开 HDF4 文件成功！")
                print("可用数据集:", hdf.datasets())
        else:
            print("文件格式未知，不支持该文件类型。")
    
    except Exception as e:
        print("发生错误:", e)

实例2：打开并读取hdf4文件（地面反照率数据）

以下载的2020年地面反照率数据为例，数据详细介绍可参见另一博客-【数据集】Global Land Surface Satellites (GLASS)：LAI、Albedo、LST、FVC。

Python代码如下：

import h5py
from pyhdf.SD import SD, SDC

# 打开 HDF 文件
hdf_file_path = 'D:/0 DataBase/8 GLASS02B03.V50_Albedo/GLASS02B03.V50.A2020001.2021300.hdf'  # 替换为你的文件路径

try:
    # 直接打开 HDF 文件
    data = SD(hdf_file_path, SDC.READ)
    print("成功打开 HDF 文件！")

    # 打印文件中的所有数据集
    datasets = data.datasets()
    print("文件中的数据集:")

    for dataset_name in datasets.keys():
        # 获取数据集信息
        dataset_info = datasets[dataset_name]
        data_type = dataset_info[0]  # 数据类型
        data_shape = dataset_info[1]  # 数据形状

        print(f"数据集名称: {dataset_name}, 数据类型: {data_type}, 数据形状: {data_shape}")

        # 读取数据集
        dataset = data.select(dataset_name)
        data_values = dataset[:]

        # 打印数据集的部分内容（前10个值）
        print(f"数据集 '{dataset_name}' 的前10个值: {data_values.flatten()[:10]}")

    # 关闭 HDF 文件
    data.end()

except OSError as e:
    print("打开 HDF 文件失败:", e)
except Exception as e:
    print("发生错误:", e)

输出结果如下：（数据空间分辨率：0.05°；空间范围：全球）

成功打开 HDF 文件！
文件中的数据集:
数据集名称: BSA_VIS, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'BSA_VIS' 的前10个值: [9360 9360 9360 9360 9360 9360 9360 9360 9360 9360]
数据集名称: BSA_NIR, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'BSA_NIR' 的前10个值: [5560 5560 5560 5560 5560 5560 5560 5560 5560 5560]
数据集名称: BSA_shortwave, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'BSA_shortwave' 的前10个值: [7960 7960 7960 7960 7960 7960 7960 7960 7960 7960]
数据集名称: WSA_VIS, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'WSA_VIS' 的前10个值: [9160 9160 9160 9160 9160 9160 9160 9160 9160 9160]
数据集名称: WSA_NIR, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'WSA_NIR' 的前10个值: [5360 5360 5360 5360 5360 5360 5360 5360 5360 5360]
数据集名称: WSA_shortwave, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'WSA_shortwave' 的前10个值: [7760 7760 7760 7760 7760 7760 7760 7760 7760 7760]
数据集名称: QC, 数据类型: ('YDim:GLASS02B03', 'XDim:GLASS02B03'), 数据形状: (3600, 7200)
数据集 'QC' 的前10个值: [32747 32747 32747 32747 32747 32747 32747 32747 32747 32747]

实例3：绘制图形（地面反照率数据）

绘制图形如下：

相关代码如下：

import h5py
from pyhdf.SD import SD, SDC
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap

# 设置字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'

# 打开 HDF 文件
hdf_file_path = 'D:/0 DataBase/8 GLASS02B03.V50_Albedo/GLASS02B03.V50.A2020001.2021300.hdf'  # 替换为你的文件路径

try:
    # 直接打开 HDF 文件
    data = SD(hdf_file_path, SDC.READ)
    print("成功打开 HDF 文件！")

    # 打印文件中的所有数据集
    datasets = data.datasets()
    print("文件中的数据集:")

    for dataset_name in datasets.keys():
        # 获取数据集信息
        dataset_info = datasets[dataset_name]
        data_type = dataset_info[0]  # 数据类型
        data_shape = dataset_info[1]  # 数据形状

        print(f"数据集名称: {dataset_name}, 数据类型: {data_type}, 数据形状: {data_shape}")

        # 提取数据集 'BSA_VIS'
        dataset_name = 'BSA_VIS'
        dataset = data.select(dataset_name)
        data_values = dataset[:].astype(np.float32)

        # 生成经纬度网格
        latitudes = np.linspace(-90, 90, data_values.shape[0])  # 3600 行
        longitudes = np.linspace(-180, 180, data_values.shape[1])  # 7200 列

        # 找出最大值和最小值
        bsa_max = np.nanmax(data_values[data_values != -1])  # 使用 nanmax 处理无效数据
        bsa_min = np.nanmin(data_values[data_values != -1])

        # 打印最大值和最小值
        print(f"BSA_VIS 最小值: {bsa_min}, 最大值: {bsa_max}")

        # 绘制数据集图形
        plt.figure(figsize=(10, 6))

        # 创建 Basemap 实例
        m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90,
                    llcrnrlon=-180, urcrnrlon=180, resolution='c')

        # 将 -1 设置为缺测值
        data_values[data_values == -1] = np.nan  # 用 NaN 替换 -1

        # 绘制陆地边界
        m.drawcoastlines()
        m.drawcountries()
        m.drawmapboundary(fill_color='lightblue')

        # 绘制经纬度网格和刻度标签
        m.drawparallels(np.arange(-90., 91., 30.), labels=[1, 0, 0, 0], fontsize=10)  # 纬度
        m.drawmeridians(np.arange(-180., 181., 60.), labels=[0, 0, 0, 1], fontsize=10)  # 经度

        # 反转纬度数组
        #latitudes = latitudes[::-1]  # 反转纬度顺序

        # 将数据投影到 Basemap 中
        x = longitudes
        y = latitudes
        data_values = np.ma.masked_invalid(data_values)  # 避免无效数据

        # 使用 contourf 绘制数据
        lon_mesh, lat_mesh = np.meshgrid(x, y)
        cmap = plt.get_cmap('viridis')
        cmap.set_bad(color='white')  # 将缺测值设为白色

        cs = m.contourf(lon_mesh, lat_mesh, data_values[::-1, :], cmap=cmap,
                        vmin=bsa_min, vmax=bsa_max)  # 设置颜色范围

        # 添加颜色条
        cbar = plt.colorbar(cs, orientation='vertical', pad=0.02)
        cbar.set_label('BSA_VIS')  # 设置颜色条标签
        cbar.ax.tick_params(labelsize=10)  # 设置颜色条刻度标签大小

        # 设置经纬度刻度标签
        m.drawmeridians(np.arange(-180, 181, 30), labels=[1, 0, 0, 0])  # 经度
        m.drawparallels(np.arange(-90, 91, 30), labels=[0, 0, 0, 1])  # 纬度

        plt.xlabel('Longitude (°)', fontsize=14, labelpad=15)
        plt.ylabel('Latitude (°)', fontsize=14, labelpad=30)
        plt.title('Surface Albedo in 2020', fontsize=16)

        # 设置坐标轴刻度值字体大小
        plt.tick_params(axis='both', labelsize=14)

        plt.show()

        # 关闭 HDF 文件
        data.end()


    # 关闭 HDF 文件
    data.end()

except OSError as e:
    print("打开 HDF 文件失败:", e)
except Exception as e:
    print("发生错误:", e)

参考