基于python的MODIS数据质量控制------以MOD11A1为例

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对地理数据进行处理，对地理模块相对熟悉，大家可以联系我,一起学习。