基于python和google earth engine的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代码进行批量处理。)

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 pd

def 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 flags
            Error_mask = np.char.count(new,'00',start=0,end=2) == 1
            print(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，后续设置为nodata
            Out_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对地理数据进行处理，对地理模块相对熟悉，大家可以联系我,一起学习。

对于待定像元如何处理

名为“Naga Climber”的粉丝把这个图片发给我了。其实处理起来也很简单，谁需要的时候可以加我qq，上面代码里面有我的联系方式。我不是要私聊收费啊，就是我会写，但是我没有应用需求，懒得写，谁需要我就给他写个。

更新的Google Earth Engine 代码

var geometry: FeatureCollection (1 element)

//就是这里，导出的格式写为SHP
Export.table.toDrive({
  collection:geometry,
  description: "geometry_studyArea_parper",
  fileNamePrefix: "geometry",
  fileFormat: "SHP"
});

  

// Helper function to extract the values from specific bits
// The input parameter can be a ee.Number() or ee.Image()
// Code adapted from https://gis.stackexchange.com/a/349401/5160
var bitwiseExtract = function(input, fromBit, toBit) {
  var maskSize = ee.Number(1).add(toBit).subtract(fromBit)
  var mask = ee.Number(1).leftShift(maskSize).subtract(1)
  return input.rightShift(fromBit).bitwiseAnd(mask)
}
 
 
var dataset = ee.ImageCollection('MODIS/061/MYD11A1')
                  .filter(ee.Filter.date('2020-02-01', '2020-02-05'));
                  
var landSurfaceTemperature = dataset.select('LST_Day_1km','QC_Day');
                  

var bitwiseExtract = function(input, fromBit, toBit) {
  var maskSize = ee.Number(1).add(toBit).subtract(fromBit);
  var mask = ee.Number(1).leftShift(maskSize).subtract(1);
  return input.rightShift(fromBit).bitwiseAnd(mask);
};


var qcControlForMod11a1 = function(qcDayLayer, lstDayLayer){
  var lstMasked = function(image){
    var qcDay = image.select(qcDayLayer);
    var lstDay = image.select(lstDayLayer);
    var qaMask = bitwiseExtract(qcDay, 0, 1).lte(1);
    var dataQualityMask = bitwiseExtract(qcDay, 2, 3).eq(0);
    var lstErrorMask = bitwiseExtract(qcDay, 6, 7).eq(0);
    var mask = qaMask.and(dataQualityMask).and(lstErrorMask);
    return lstDay.updateMask(mask);
  };
  return lstMasked; // this is very important!!!
};

var mod11a1WithQc = landSurfaceTemperature.map(qcControlForMod11a1('QC_Day', 'LST_Day_1km'));


var vis = {
  min: 13000.0,
  max: 16500.0,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};

Map.addLayer(landSurfaceTemperature.select('LST_Day_1km').first().clip(geometry), vis, 'Original_1');

Map.addLayer(mod11a1WithQc.select('LST_Day_1km').first().clip(geometry), vis, 'After QC_1');

Map.addLayer(ee.Image(landSurfaceTemperature.select('LST_Day_1km').toList(landSurfaceTemperature.size()).get(1))
  .clip(geometry), vis, 'Original_2');

Map.addLayer(ee.Image(mod11a1WithQc.select('LST_Day_1km').toList(landSurfaceTemperature.size()).get(1))
  .clip(geometry), vis, 'After QC_2');


print(landSurfaceTemperature.size());


var list = mod11a1WithQc.select('LST_Day_1km').toList(landSurfaceTemperature.size()); 


for (var i=0; i<367; i++){
  
  var image = ee.Image(list.get(i));
  var id1 = image.get("system:index");
  
  // var stringALB = ee.String('ERA5_LAND_01Dgre_mm_HB_Daily_Prec_');
  // var string1 = ee.String(id1);
  // var nameOut = (stringALB.cat(string1)).getInfo();
  
  var nameOut = 'MYD11A1_'+  // 导出数据
  Export.image.toDrive({
    image: out.unmask(-9999),
    description: nameOut, //数据名
    folder:'MYD11A1', //存储文件夹
    scale: 1000, //分辨率
    region:geometry,
    maxPixels:1e13,
  });id1.getInfo()+'_LST_Day_1km';
  
  print(nameOut);
  
  // var nameOut = 'ERA5_LAND_01Dgre_K_HB_Daily_Temp_' +  id1.getInfo();
  
  var out =image.clip(geometry)//.float().multiply(0.02);
  
  // // print(out);
  
  // 导出数据
  Export.image.toDrive({
    image: out.unmask(-9999),
    description: nameOut, //数据名
    folder:'MYD11A1', //存储文件夹
    scale: 1000, //分辨率
    region:geometry,
    maxPixels:1e13,
  });
  
}

代码运行结果