Python中gdal读取多波段HDF栅格影像并绘制直方图

  本文详细介绍基于Python语言gdal等模块实现多波段HDF栅格图像文件（即.hdf文件）的读取、处理与像元值可视化等操作。此外，基于gdal等模块读取.tif格式栅格图层文件的方法可以查看Python中gdal实现多幅栅格影像批量绘制直方图，读取单波段.hdf格式栅格图层文件的方法可以查看Python中gdal栅格影像读取计算与写入及质量评估QA波段筛选掩膜。

  本文期望实现的需求为：现有一存放.tif格式的全球LAI产品栅格数据的路径，需将这一路径下的全部LAI产品栅格数据依据另一路径下存放的全球MODIS植被覆盖类型产品栅格数据进行像元分类，并绘制全球每一种植被类型对应的LAI数值直方图。在这里，由于有前述两篇推文作为铺垫，本文对代码的讲解就着重于多波段HDF栅格图像文件的读取部分；其它内容由于前面两篇推文已经详细介绍，这里就不再赘述~

  首先将本文所需代码展示如下：

# -*- coding: utf-8 -*-
"""
Created on Tue Jul 20 11:05:31 2021

@author: fkxxgis
"""

import os
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt

lai_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/h20v09.tif"
mcd_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/MCD12Q1.A2018001.h20v09.006.2019199233851.hdf"
pic_save_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/"

for veg_type in range(9):

    mcd_raster=gdal.Open(mcd_file_path)
    mcd_sub_dataset=mcd_raster.GetSubDatasets()
    hdf_band_num=len(mcd_sub_dataset)
    # for sub_dataset in mcd_sub_dataset:
    #     print(sub_dataset[1])
    # print(mcd_sub_dataset[2][1])
    mcd_sub_type=gdal.Open(mcd_sub_dataset[2][0])
    mcd_raster_array=mcd_sub_type.ReadAsArray()

    lai_raster=gdal.Open(lai_file_path)
    lai_raster_array=lai_raster.ReadAsArray()
    non_veg_type_lai_array=np.where(mcd_raster_array==veg_type+1,lai_raster_array,np.nan)
    plt.hist(non_veg_type_lai_array)
    plt.savefig(pic_save_path+"DRT_"+str(veg_type+1)+".png", dpi=300)
    plt.clf()
    plt.cla()

  我们直接讲解多波段HDF栅格图像文件读取部分的代码：首先，多波段.hdf格式文件的读取在一开始与单波段.hdf格式文件或.tif格式文件的读取一致，即通过gdal.Open()函数实现；但随后，需要额外借助len()函数获取HDF文件对应的波段数量。

  因为我们读取的HDF文件是多波段，因此hdf_band_num肯定是大于1的，那么刚刚读取进来的mcd_sub_dataset其实就是一个列表（List）；其中，这个列表的元素个数就是对应的多波段HDF文件波段数，列表的每一个元素则都是一个元组（tuple）；同时，每一个元组都有两个元素，其每一个元素都是一个字符串；其中第一个元素为当前HDF文件的当前波段对应的文件路径与部分提示信息，第二个元素作为当前HDF文件的当前波段对应的文件像素行列数、名称与数据类型。

  这么说可能不太明白，我们用一个实例来讲解。例如，通过上述代码读取一景具有六个波段的MODIS LAI产品——MCD15A3H产品，其第一个波段为FPAR数据，第二个波段为LAI数据。那么读取其后，得到的mcd_sub_dataset长这个样子：

  可以看到，是一个具有6个元素的列表。

  点开列表，可以看到6个元素每一个都是一个具有2个元素的元组：

  再点开第一个元组，可以看到其具有2个字符串格式的元素：

  其第二个元素包含了该波段对应的数据行数与列数（即[2400×2400]）、数据名称（即Fpar）、数据空间分辨率（即500m）、数据产品简称（即MOD_Grid_MCD15A3H），以及数据格式（即8-bit unsigned integer）；而第一个字符串没有显示完毕，我们可以点击打开看看：

  可以看到第一个元素则包含了该波段对应的数据路径、文件全称，以及部分与第二个元素重复的几个数据信息参数。

  有了上面的分析就比较清楚了，接下来再一次利用gdal.Open()函数读取我们需要的波段，mcd_sub_dataset[2][0]表示第三个波段；其中，第三个波段却用[2]来表示，是因为波段数量（也就是mcd_sub_dataset的Index）是从0开始计算的；而后面的[0]则表示元组中的第一个参数，也就是上面一幅图中显示的该波段对应的数据路径。

  随后，再利用.ReadAsArray()函数将其读取为Array即可。接下来的操作与本文开头提及的那两篇推文就一致，这里不再赘述~

欢迎关注公众号：疯狂学习GIS