Python gdal讀取MODIS遙感影像並結合質量控制QC波段掩膜資料
本文介紹基於Python中GDAL模組，實現MODIS遙感影像資料的讀取、計算，並基於質量控制QC波段進行影象掩膜的方法。
前期的文章Python GDAL讀取柵格資料並基於質量評估波段QA對指定資料加以篩選掩膜詳細介紹了基於Python語言gdal等模組實現遙感影像柵格資料的讀取，以及基於質量評估（QA）波段實現柵格像元篩選與掩膜的全部操作。而在本文，我們依據前述這一篇部落格的程式碼，結合大家更為熟悉的MODIS系列遙感影像產品，基於其質量評估波段進行具體的對照講解。也就是說，本文重點不在於程式碼的講解（具體程式碼在前述這一篇部落格中已經很詳細地介紹了），而是將上述程式碼在更為具體的一個實踐中加以應用，告訴大家該如何選擇波段、處理質量評估QA波段並進行篩選操作等。同時，這裡還有一點需要注意：在MODIS系列遙感影像中，質量評估波段更應該稱為質量控制波段，因為其官方手冊中將其寫作Quality Control，因此後文就寫作質量控制波段或QC波段。
首先，需要下載好對應的MODIS資料，大家可以依據文章批次下載MODIS遙感影像：基於LAADS DAAC的方法中的方法進行下載。本文就以一景MODIS的LAI產品——MCD15A3H產品為例進行操作。
下載後，開啟HDF檔案可以看到，其具有很多波段，同時包括質量控制QC波段；且在FPAR與LAI波段中，像元數值方面還具有精度較低的像元值、填充值等無效數值。上述這些都需要我們在讀取資料時加以識別、處理與篩選。
由於MODIS系列遙感影像產品種類較多，不同產品之間的屬性差異較大；因此建議大家每次使用一種MODIS產品時，都到官網檢視其基本資訊，有需要的話還可以在官網下載對應產品的使用者手冊。前面提到，本文所用產品為MCD15A3H，因此可以在其官網查閱其基本資訊或下載使用者手冊檢視更為詳細的產品屬性。
例如，下圖所示即為使用者手冊中關於這一產品一景影像中波段分佈情況與每一個波段具體資訊的介紹表格；其中包括了波段含義、資料型別、填充值範圍、有效值範圍與縮放係數等關鍵引數，這些對於後期我們用gdal讀取.hdf格式柵格檔案而言具有重要意義。
接下來，質量控制QC波段同樣是執行柵格讀取操作前有必要了解的資訊。下圖所示即為使用者手冊中關於這一產品一景影像中質量控制QC波段具體資訊介紹的表格，其中包含了當前一景影像中FPAR與LAI產品的每一個像元所對應的演演算法、感測器、雲覆蓋等資訊。這裡需要注意的是：在MCD15A3H產品中是有2個質量控制QC波段的，這個是第一個QC，而第二個QC主要包括水陸區域、冰雪區域、氣溶膠等資訊，本文中暫且不涉及第二個QC。
其中，由上表可知，QC波段的資訊一共是由0至7共8個位元位（即Bit No.）組成，其中，由若干個位元位又可以組成Bit-word，每一個Bit-word就代表某一種QC波段資訊。結合上圖，我們可以對照下圖這樣一個範例進行理解：
結合以上基本資訊，我們已經對MCD15A3H產品的基本資訊有了一定了解。接下來就可以進行柵格資料的讀取與處理、篩選了。
在這裡需要注意的是，之前的兩篇文章Python GDAL讀取柵格資料並基於質量評估波段QA對指定資料加以篩選掩膜以及Python批次讀取HDF多波段柵格資料並繪製像元直方圖已經對本次所要用到的大部分需求與程式碼加以實現並進行了詳細講解，這裡就不再贅述。本文程式碼所實現功能與上述第一篇部落格中的需求一致，唯一不同的是將GLASS產品更改為了MCD15A3H產品，且僅需對MCD15A3H產品的主演演算法像元加以做差計算（也就是篩選出MCD15A3H產品中第一個QC波段對應二進位制數的第一位為0的像元，其它像元就不用參與差值計算了）。
具體程式碼如下：
# -*- coding: utf-8 -*-
"""
Created on Sun Jul 25 14:57:45 2021
@author: fkxxgis
"""
import os
import copy
import numpy as np
from osgeo import gdal
rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/RT_LAI/"
mcd15_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/mcd15A3H/"
out_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/"
rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
rt_file_split=rt_file.split("_")
rt_hv=rt_file_split[3][:-4]
mcd15_file_list=os.listdir(mcd15_file_path)
for mcd15_file in mcd15_file_list:
if rt_hv in mcd15_file:
rt_file_tif_path=rt_file_path+rt_file
mcd15_file_tif_path=mcd15_file_path+mcd15_file
drt_out_file_path=out_file_path+"drt/"
if not os.path.exists(drt_out_file_path):
os.makedirs(drt_out_file_path)
drt_out_file_tif_path=drt_out_file_path+rt_hv+".tif"
eco_out_file_path=out_file_path+"eco/"
if not os.path.exists(eco_out_file_path):
os.makedirs(eco_out_file_path)
eco_out_file_tif_path=eco_out_file_path+rt_hv+".tif"
wat_out_file_path=out_file_path+"wat/"
if not os.path.exists(wat_out_file_path):
os.makedirs(wat_out_file_path)
wat_out_file_tif_path=wat_out_file_path+rt_hv+".tif"
tim_out_file_path=out_file_path+"tim/"
if not os.path.exists(tim_out_file_path):
os.makedirs(tim_out_file_path)
tim_out_file_tif_path=tim_out_file_path+rt_hv+".tif"
rt_raster=gdal.Open(rt_file_tif_path)
rt_raster_array=rt_raster.ReadAsArray()
rt_lai_array=rt_raster_array[0]
rt_qa_array=rt_raster_array[1]
rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
rt_lai_array_fin=rt_lai_array_mask*0.001
mcd15_raster=gdal.Open(mcd15_file_tif_path)
mcd15_sub_dataset=mcd15_raster.GetSubDatasets()
# for sub_dataset in mcd15_sub_dataset:
# print(sub_dataset[1])
# print(mcd15_sub_dataset[1][1])
# print(mcd15_sub_dataset[2][1])
mcd15_sub_lai=gdal.Open(mcd15_sub_dataset[1][0])
mcd15_sub_qc=gdal.Open(mcd15_sub_dataset[2][0])
mcd15_lai_array=mcd15_sub_lai.ReadAsArray()
mcd15_qc_array=mcd15_sub_qc.ReadAsArray()
mcd15_lai_array_mask=np.where(mcd15_lai_array>248,np.nan,mcd15_lai_array)
mcd15_lai_array_fin=mcd15_lai_array_mask*0.1
rt_qa_array_bin=copy.copy(rt_qa_array)
rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
for i in range(rt_qa_array_row):
for j in range(rt_qa_array_col):
rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
mcd15_qc_array_bin=copy.copy(mcd15_qc_array)
mcd15_qc_array_row,mcd15_qc_array_col=mcd15_qc_array.shape
for i in range(mcd15_qc_array_row):
for j in range(mcd15_qc_array_col):
mcd15_qc_array_bin[i][j]="{:08b}".format(mcd15_qc_array[i][j])[-1:]
mcd15_lai_main_array=np.where(mcd15_qc_array_bin==1,np.nan,mcd15_lai_array_fin)
lai_dif=rt_lai_array_fin-mcd15_lai_main_array
lai_dif=lai_dif*1000
drt_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
np.nan,lai_dif)
eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
np.nan,lai_dif)
wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
np.nan,lai_dif)
tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
np.nan,lai_dif)
row=rt_raster.RasterYSize
col=rt_raster.RasterXSize
geotransform=rt_raster.GetGeoTransform()
projection=rt_raster.GetProjection()
# 輸出為int格式後，所得結果中0就是NoData
driver=gdal.GetDriverByName("Gtiff")
out_drt_lai=driver.Create(drt_out_file_tif_path,row,col,1,gdal.GDT_Int16)
out_drt_lai.SetGeoTransform(geotransform)
out_drt_lai.SetProjection(projection)
out_drt_lai.GetRasterBand(1).WriteArray(drt_lai_dif_array)
out_drt_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Int16)
out_eco_lai.SetGeoTransform(geotransform)
out_eco_lai.SetProjection(projection)
out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
out_eco_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Int16)
out_wat_lai.SetGeoTransform(geotransform)
out_wat_lai.SetProjection(projection)
out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
out_wat_lai=None
driver=gdal.GetDriverByName("Gtiff")
out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Int16)
out_tim_lai.SetGeoTransform(geotransform)
out_tim_lai.SetProjection(projection)
out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
out_tim_lai=None
print(rt_hv)
至此，大功告成。