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

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

本文期望實(shí)現(xiàn)的需求為：現(xiàn)有一存放.tif格式的全球LAI產(chǎn)品柵格數(shù)據(jù)的路徑，需將這一路徑下的全部LAI產(chǎn)品柵格數(shù)據(jù)依據(jù)另一路徑下存放的全球MODIS植被覆蓋類(lèi)型產(chǎn)品柵格數(shù)據(jù)進(jìn)行像元分類(lèi)，并繪制全球每一種植被類(lèi)型對(duì)應(yīng)的LAI數(shù)值直方圖。在這里，由于有前述兩篇博客作為鋪墊，本文對(duì)代碼的講解就著重于多波段HDF柵格圖像文件的讀取部分；其它內(nèi)容由于在本文開(kāi)頭提及的前期兩篇博客中已經(jīng)詳細(xì)介紹，這里就不再贅述~

首先將本文所需代碼展示如下：

# -*- 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格式文件的讀取在一開(kāi)始與單波段.hdf格式文件或.tif格式文件的讀取一致，即通過(guò)gdal.Open()函數(shù)實(shí)現(xiàn)；但隨后，需要額外借助len()函數(shù)獲取HDF文件對(duì)應(yīng)的波段數(shù)量。

因?yàn)槲覀冏x取的HDF文件是多波段，因此hdf_band_num肯定是大于1的，那么剛剛讀取進(jìn)來(lái)的mcd_sub_dataset其實(shí)就是一個(gè)列表（List）；其中，這個(gè)列表的元素個(gè)數(shù)就是對(duì)應(yīng)的多波段HDF文件波段數(shù)，列表的每一個(gè)元素則都是一個(gè)元組（tuple）；同時(shí)，每一個(gè)元組都有兩個(gè)元素，其每一個(gè)元素都是一個(gè)字符串；其中第一個(gè)元素為當(dāng)前HDF文件的當(dāng)前波段對(duì)應(yīng)的文件路徑與部分提示信息，第二個(gè)元素作為當(dāng)前HDF文件的當(dāng)前波段對(duì)應(yīng)的文件像素行列數(shù)、名稱(chēng)與數(shù)據(jù)類(lèi)型。

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

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

點(diǎn)開(kāi)列表，可以看到6個(gè)元素每一個(gè)都是一個(gè)具有2個(gè)元素的元組：

再點(diǎn)開(kāi)第一個(gè)元組，可以看到其具有2個(gè)字符串格式的元素：

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

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

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

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

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

最新評(píng)論