Google Earth Engine (GEE) 实现对MODIS产品批量质量控制

前言

最近处在学位论文初稿完成和申博面试之间的空档期，所以想学点新东西。GEE火了很久了，一直没有真正使用过，还是坚持下载数据到本地，然后用Python处理。主要是怕一旦接触GEE，就会花费很多心思在上面。再者，我本身是做算法的，不仅仅是数据分析，需要使用本地数据来应用算法。现在觉得下载到本地然后处理数据有很多不确定性，不如直接在GEE上完成预处理然后导出，可以省多很多精力给算法研究本身。所以，今后主要关注于使用GEE进行数据预处理和导出的相关实现，看了一下午官方文档和一些资料，发现GEE对于新手还是很友好的。如果有编程基础，很快就能入门。

言归正传，https://spatialthoughts.com/2021/08/19/qa-bands-bitmasks-gee/实现了利用QC影像对单幅image质量控制。本文将以MODIS地表温度产品MOD11A1为例，实现对一个时间段内的imageCollection进行批量质量控制。内容涉及到GEE中多参数嵌套map的原理和实现。

提取指定位置qc二进制码函数

对于qc码不再赘述，上面的链接包括我之前的博客都有介绍。MOD11A1的qc码共4对，8位，每一对包含不同方面的产品质量。下面这个函数就负责提取指定位置的qc码的十进制数字，来实现对产品不同方面的质量控制。

// 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
// the decmial to binary: 0-00, 1-01, 2-10, 3-11, 4-100
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);
};

对单幅image进行质量控制

想要实现对整个image进行质量控制。需要包含两个步骤，提取4对质量控制码的十进制数字并建立掩膜；将掩膜应用于image。
下面的代码来自于前言中的链接，运行的时候别忘了把上面的提取指定位置qc二进制码函数粘贴进来。

var modisLST = ee.ImageCollection("MODIS/006/MOD11A1")
var lsib = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017")
var australia = lsib.filter(ee.Filter.eq('country_na', 'Australia'))
var geometry = australia.geometry()
var terra = modisLST
  .filter(ee.Filter.date('2001-01-01', '2010-01-01'))
  .select('LST_Day_1km','QC_Day');
   
// Get a single image for testing
var image = ee.Image(terra.first())
var lstDay = image.select('LST_Day_1km')
var qcDay = image.select('QC_Day')
// Let's extract all pixels from the input image where
// Bits 0-1 <= 1 (LST produced of both good and other quality)
// Bits 2-3 = 0 (Good data quality)
// Bits 4-5 Ignore, any value is ok
// Bits 6-7 = 0 (Average LST error ≤ 1K)
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)
var lstDayMasked = lstDay.updateMask(mask)  
var visParams = {min:13000, max:16000, palette: ['green', 'yellow', 'red']}
Map.addLayer(lstDay.clip(geometry), visParams, 'Original LST Image');
Map.addLayer(lstDayMasked.clip(geometry), visParams, 'LST Masked');

对imageCollection质量控制的嵌套函数

上面实现了对单个image的qc影像提取的功能，想要实现对整个imageCollection进行质量控制。需要将对单幅image进行质量控制的过程封装成函数，然后应用map函数将其应用于imageCollection中每一幅image。
这里的难点在于，需要质量控制的波段和qc波段的名称不一致，需要向map函数中传递两个参数来指定imageCollection中的波段，所以要用到嵌套函数。下面是函数的代码，我写了比较详细的英文注释来解释嵌套函数原理。这里再用中文解释一下，外部函数负责将波段名称传递到内部函数，内部函数包含了我们真正想要执行的功能，并返回处理好的image，外部函数还负责返回内部函数。

// the principle is that
// external function play the role to transform parameters, such as the qcDayLayer and lstDayLayer
// external function also play the role to return the internal function, such as the lstMasked
// after execute the external function, 
// the internal function will executed the same as the normal function
// if you want achieve the same goal, please put the final results that you want to the internal function,
// and return it 
// the external function only play the role to transform parameters to the internal function
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!!!
};

下面这一句是imageCollection执行qc控制的语句，mod11a1是MOD11A1的imageCollection。

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

qcControlForMod11a1('QC_Day', 'LST_Day_1km')负责将qc和白天地表温度波段的名称传入了qcControlForMod11a1函数，qcControlForMod11a1接收到这两个参数后，会返回lstMasked。而lstMasked是一个函数，所以接下来的过程就和普通map函数执行的过程一致了，即对imageCollection中每一幅image执行lstMasked，然后返回经过掩膜后的image。

完整代码

// 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
// the decmial to binary: 0-00, 1-01, 2-10, 3-11, 4-100
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);
};

// the principle is that
// external function play the role to transform parameters, such as the qcDayLayer and lstDayLayer
// external function also play the role to return the internal function, such as the lstMasked
// after execute the external function, 
// the internal function will executed the same as the normal function
// if you want achieve the same goal, please put the final results that you want to the internal function,
// and return it 
// the external function only play the role to transform parameters to the internal function
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 imageCollection = ee.ImageCollection("MODIS/006/MOD11A1")

var lsib = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017");
var china = lsib.filter(ee.Filter.eq('country_na', 'China'));
var geometry = china.geometry();

var mod11a1 = imageCollection.filterDate('2021-01-01', '2021-01-03').select('LST_Day_1km', 'QC_Day');
print(mod11a1);

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

var vis = {
  min: 12500,
  max: 15500,
  palette: [
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef', '3ae237',
    'b5e22e', 'd6e21f', 'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08',
    'ff500d', 'ff0000', 'de0101', 'c21301'
  ],
};

Map.addLayer(mod11a1.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(mod11a1.select('LST_Day_1km').toList(mod11a1.size()).get(1))
  .clip(geometry), vis, 'Original_2');

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

应用结果

可以看到01和02日的image都成功实现了质量控制，后续就可以导出数据。
之后会再更新一些如何导出WGS84地理坐标系下指定空间分辨率影像的教程，有需要的同学可以先关注O(∩_∩)O。