本文介紹基於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)
至此，大功告成。
