??本文介紹基于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)
??至此，大功告成，
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/546430.html
標籤：其他
上一篇：AI智能問答助手 AI智能批量文章生成器 網站優化SEO批量內容生成工具 原創文章生成軟體
下一篇：【機器學習】1. 廣義線性模型