本文介紹基于Python語言gdal模塊，實現(xiàn)多波段HDF柵格圖像文件的讀取、處理與像元值可視化（直方圖繪制）等操作。

另外，基于gdal等模塊讀取.tif格式柵格圖層文件的方法可以查看Python批量繪制遙感影像數(shù)據(jù)的直方圖，讀取單波段.hdf格式柵格圖層文件的方法可以查看Python GDAL讀取柵格數(shù)據(jù)并基于質量評估波段QA對指定數(shù)據(jù)加以篩選掩膜。

本文期望實現(xiàn)的需求為：現(xiàn)有一存放.tif格式的全球LAI產品柵格數(shù)據(jù)的路徑，需將這一路徑下的全部LAI產品柵格數(shù)據(jù)依據(jù)另一路徑下存放的全球MODIS植被覆蓋類型產品柵格數(shù)據(jù)進行像元分類，并繪制全球每一種植被類型對應的LAI數(shù)值直方圖。在這里，由于有前述兩篇博客作為鋪墊，本文對代碼的講解就著重于多波段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()函數(shù)實現(xiàn)；但隨后，需要額外借助len()函數(shù)獲取HDF文件對應的波段數(shù)量。

因為我們讀取的HDF文件是多波段，因此hdf_band_num肯定是大于1的，那么剛剛讀取進來的mcd_sub_dataset其實就是一個列表（List）；其中，這個列表的元素個數(shù)就是對應的多波段HDF文件波段數(shù)，列表的每一個元素則都是一個元組（tuple）；同時，每一個元組都有兩個元素，其每一個元素都是一個字符串；其中第一個元素為當前HDF文件的當前波段對應的文件路徑與部分提示信息，第二個元素作為當前HDF文件的當前波段對應的文件像素行列數(shù)、名稱與數(shù)據(jù)類型。

這么說可能不太明白，我們用一個實例來講解。例如，通過上述代碼讀取一景具有六個波段的MODIS LAI產品——MCD15A3H產品，其第一個波段為FPAR數(shù)據(jù)，第二個波段為LAI數(shù)據(jù)。那么讀取其后，得到的mcd_sub_dataset長這個樣子：

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

點開列表，可以看到6個元素每一個都是一個具有2個元素的元組：

再點開第一個元組，可以看到其具有2個字符串格式的元素：

其第二個元素包含了該波段對應的數(shù)據(jù)行數(shù)與列數(shù)（即[2400×2400]）、數(shù)據(jù)名稱（即Fpar）、數(shù)據(jù)空間分辨率（即500m）、數(shù)據(jù)產品簡稱（即MOD_Grid_MCD15A3H），以及數(shù)據(jù)格式（即8-bit unsigned integer）；而第一個字符串沒有顯示完畢，我們可以點擊打開看看：

可以看到第一個元素則包含了該波段對應的數(shù)據(jù)路徑、文件全稱，以及部分與第二個元素重復的幾個數(shù)據(jù)信息參數(shù)。

有了上面的分析就比較清楚了，接下來再一次利用gdal.Open()函數(shù)讀取我們需要的波段，mcd_sub_dataset[2][0]表示第三個波段；其中，第三個波段卻用[2]來表示，是因為波段數(shù)量（也就是mcd_sub_dataset的Index）是從0開始計算的；而后面的[0]則表示元組中的第一個參數(shù)，也就是上面一幅圖中顯示的該波段對應的數(shù)據(jù)路徑。

隨后，再利用.ReadAsArray()函數(shù)將其讀取為Array即可。接下來的操作與本文開頭提及的那兩篇博客就一致，這里不再贅述~

到此這篇關于Python實現(xiàn)批量讀取HDF多波段柵格數(shù)據(jù)并繪制像元直方圖的文章就介紹到這了,更多相關Python讀取HDF多波段柵格數(shù)據(jù)內容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家！

最新評論