Python批次讀取HDF多波段柵格資料並繪製像元直方圖

瘋狂學習GIS發表於2023-03-01
Python直方圖

  本文介紹基於Python語言gdal模組，實現多波段HDF柵格影像檔案的讀取、處理與像元值視覺化（直方圖繪製）等操作。

  另外，基於gdal等模組讀取.tif格式柵格圖層檔案的方法可以檢視Python批次繪製遙感影像資料的直方圖，讀取單波段.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即可。接下來的操作與本文開頭提及的那兩篇部落格就一致，這裡不再贅述~

相關文章