??本文介紹基于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即可，接下來的操作與本文開頭提及的那兩篇博客就一致，這里不再贅述~
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/545397.html
標籤：Python