MODIS質量控制文件，對MODIS產品進行提取

MODIS數據簡介

我們拿到的MODIS數據，多數人認為只要有值的地方，就是準確數據，我們直接就可以拿來使用，只有空值的區域，數據才會異常（多數本科生是這樣認為的）；然而并非如此，往往一個MODIS產品一個像元處，只有當所有輸入的反演參數都為異常值時，這個像元才會被設置為異常，即設置為空值。 因此，我們所能看到的擁有像元值的地方，就會因為輸入的反演參數都為異常程度，會有不同的質量。MODIS數據的生產商，也考慮到了數據生產過程中的數據異常情況，為了讓客戶能夠更好的使用數據，為此提供了質量空值文件（Ｑc，Ｑa）。這些信息的進一步了解,可以查看官方提供的pdf文檔，如 MOD11_User_Guide_V6.pdf. 質量空值文件多以二進制形式進行編碼，并且將多個數據圖層的質量控制參數，編寫在HDF文件的一個數據集中．本文章以MOD11A1陸表溫度日產品為例，使用python讀取二進制文件數據，以掩膜的形式將滿足要求的柵格值，提取到一個新的TIF文件中，供后續進一步使用。下圖分別為一個HDF數據的圖層(數據集)分布和QC_Day白天質量控制文件。

每一個產品像元對應一個質量控制圖層的像元.每個質量控制像元包含一個8位的整型數值,我們需要將其轉化位二進制數值,才能進行讀取\解譯.下圖為解譯的示意圖:

如圖所示,從左開始0\1位代表Mandatory QA flags，2\3位代表Data quality flag，4\5位為Emis Error flag，6\7位代表LST LST Error flag。(在python代碼中,0位在最左邊)。基于上面的理論，我們使用python讀取QC_Day的tif文件(由于前期涉及到鑲嵌\投影等步驟,所以使用MRT軟件，將HDF數據圖層,轉化為TIF文件，然后再使用python代碼進行批量處理。)

代碼

下面貼上代碼: #!/usr/bin/python # -*- coding: utf-8 -*- """本代碼，將質量控制文件，進行提取,將滿足要求的保存為tif。質量 LST error flag <= 1k@version: Anaconda @author: LeYongkang @contact: 1363989042@qq.com @software: PyCharm @file: 12_CombineLst_PredictedAndModis.py @time: 2020/2/5 0005 下午 11:10 """ import numpy as np from osgeo import gdal import os import pandas as pddef opentif(filepath):"""輸入文件名，返回數組、寬度、高度、仿射矩陣信息、投影信息:param filepath: 文件的完整路徑:return: im_data,im_width,im_height,im_geotrans,im_proj"""dataset = gdal.Open(filepath)im_width = dataset.RasterXSize #柵格矩陣的列數im_height = dataset.RasterYSize #柵格矩陣的行數data = dataset.ReadAsArray(0,0,im_width,im_height)#獲取數據# data = dataset.ReadAsArray() # 獲取數據im_data = np.array(data)print("opentif的shape")print(im_data.shape)im_geotrans = dataset.GetGeoTransform()#獲取仿射矩陣信息im_proj = dataset.GetProjection()#獲取投影信息return(im_data,im_width,im_height,im_geotrans,im_proj)def savetif(dataset,path,im_width,im_height,im_geotrans,im_proj):"""將數組保存為tif文件:param dataset: 需要保存的數組:param path: 需要保存出去的路徑，包含文件名:param im_width: 數組寬度:param im_height: 數組寬度:param im_geotrans: 仿射矩陣信息:param im_proj: 投影信息:return:"""print(dataset)driver = gdal.GetDriverByName("GTiff")outdataset = driver.Create(path, im_width, im_height, 1, gdal.GDT_Float32)print(path)outdataset.SetGeoTransform(im_geotrans) # 寫入仿射變換參數outdataset.SetProjection(im_proj) # 寫入投影outdataset.GetRasterBand(1).WriteArray(dataset)outdataset.GetRasterBand(1).SetNoDataValue(0)print("yes")if __name__ == "__main__":inDir = r"I:\2018_Parper2\MYD_2018\MYD11A1\.GeoTif_Mosaic\.Tif"# inFile = "MOD11A1.A2018001.QC_Day.tif"Out_Dir = r"I:\2018_Parper2\MYD_2018\MYD11A1\.GeoTif_Mosaic\.Tif_Masked"# 獲取LST質量控制文件的列表InList_Qc = [infile for infile in os.listdir(inDir) if infile.endswith(".QC_Day.tif")]for InFile in InList_Qc:##################################################################################################if not os.path.exists(Out_Dir +os.sep + InFile[:-10] + "LST_Day_1km.tif"):################################################################################################### 獲取完整路徑in_Full_Dir = inDir + os.sep + InFile# 打開TIF文件，獲取TIF文件的信息InData = opentif(in_Full_Dir)in_Array = InData[0]in_Array= np.array(in_Array,dtype = np.uint8)print(in_Array)print(" ")# 將十進制轉回到二進制binary_repr_v = np.vectorize(np.binary_repr)new = binary_repr_v(in_Array, 8)print(new)# 6-7位，是控制LST質量的字段，‘00’代表 LST error flag <= 1k# start=0,end=2：代表 LST LST Error flag# start=2,end=4：代表 Emis Error flag# start=4,end=6：代表Data quality flag# start=6,end=8：代表Mandatory QA flagsError_mask = np.char.count(new,'00',start=0,end=2) == 1print(Error_mask)# 打開LST文件，獲取文件名# G:\2018\MODIS\MOD11A1\.GeoTif_Mosaic_10000\Tif\MOD11A1.A2018002.QC_Day.tif 質量控制文件# G:\2018\MODIS\MOD11A1\.GeoTif_Mosaic_10000\Tif\MOD11A1.A2018002.LST_Day_1km.tif LST文件in_Full_Dir_Lst = in_Full_Dir[:-10] + "LST_Day_1km.tif"Lst_Array = opentif(in_Full_Dir_Lst)[0]# 將滿足質量條件的提取出來，不滿足條件的設置為0，后續設置為nodataOut_Lst_Array = np.where(Error_mask,Lst_Array,0)print(Out_Lst_Array)print(in_Full_Dir_Lst.split("\\")[-1])# 將masked后的LST保存，將 0 設置為SetNoDataValue()if not os.path.exists(Out_Dir + os.sep + in_Full_Dir_Lst.split("\\")[-1]):print(os.path.exists(Out_Dir + os.sep + in_Full_Dir_Lst.split("\\")[-1]))savetif(Out_Lst_Array,Out_Dir + os.sep + in_Full_Dir_Lst.split("\\")[-1],InData[1],InData[2],InData[3],InData[4])

代碼運行結果

上圖為運行結果展示，彩色為 Lst error flag <= 1k,底圖為未使用代碼提取的所有LST像元點。

官方也提供了基于arcgis的python工具箱方法(arcgis-modis-python-toolbox-v1.0)\LDOPE-1.7軟件,但是本人短時間也沒搞明白批量處理. 此博客,為本人第一次編寫，若有錯誤不妥之處，還望批評指正。此外,本人較多使用python對地理數據進行處理，對地理模塊相對熟悉，大家可以聯系我,一起學習。

總結

以上是生活随笔為你收集整理的基于python的MODIS数据质量控制------以MOD11A1为例的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。