作者: _养乐多_

本文介绍了几个 Google Earth Engine (GEE) 平台中常用的处理遥感数据中的缺失值代码片段,这些代码可以用于在时间序列中对遥感图像进行线性插值,提供更加连续和完整的时间序列。

第一段代码是一个线性插值函数,它能够在一个时间序列图像集合中对缺失数据进行插值。该函数使用了 GEE 平台中的图像集合操作函数和图像合成函数,并采用了线性插值方法对缺失值进行处理,最后返回一个插值后的图像集合。(该代码可以应用与去云后通过时间序列数据进行插值填补影像空洞)

第二段代码是一个掩模替换函数,它的作用是将原始影像中的某些值(如缺失值)用一个新值替换。它使用了 GEE 平台中的图像操作函数,并使用掩膜对影像进行处理,最后返回替换后的影像。

第三段代码是一个添加时间属性的函数,它会在影像中添加一个时间波段,用于后续时间序列分析。该函数同样使用了 GEE 平台中的图像操作函数和元数据属性获取函数。

最后一段代码是一个提取质量控制 (QA) 位的函数,它能够从一个多波段图像中提取指定位置的 QA 位数据。它使用了位运算函数和移位函数,并使用了 GEE 平台中的图像操作函数。


目录

 一、线性插值函数

二、掩模替换函数

三、添加时间属性的函数

四、提取质量控制 (QA) 位函数

五、应用案例


 一、线性插值函数

使用的时候,直接调用要插值的影像集合YourImageCollection,注意这里需要单波段影像集合,需要使用select('InterpBandName')函数筛选一下。

调用的时候直接

var nodata = -9999;
var frame = 8*4;
var bandname = 'InterpBandName';
var InterpedCollection = linearInterp(YourImageCollection, frame, nodata, bandname);

print('interped image collection:', InterpedCollection );
//线性插值填补无云像素
function linearInterp(imgcol, frame, nodata, bandname){
    frame  = frame  || 32;
    nodata = nodata || 0;

    var time   = 'system:time_start';
    imgcol = imgcol.map(addTimeBand);

    var maxDiff = ee.Filter.maxDifference(frame * (1000*60*60*24), time, null, time);
    var cond    = {leftField:time, rightField:time};

    // 待插影像之后的所有影像
    // Images after 按时间降序排列 (最近的是排序后最后一个)
    var f1 = ee.Filter.and(maxDiff, ee.Filter.lessThanOrEquals(cond));
    var c1 = ee.Join.saveAll({matchesKey:'after', ordering:time, ascending:false})
        .apply(imgcol, imgcol, f1);
    
    // 待插影像之前的所有影像
    // Images before, 按时间升序排序 (最近的是排序后最后一个)
    var f2 = ee.Filter.and(maxDiff, ee.Filter.greaterThanOrEquals(cond));
    var c2 = ee.Join.saveAll({matchesKey:'before', ordering:time, ascending:true})
        .apply(c1, imgcol, f2);
    
    var interpolated = ee.ImageCollection(c2.map(function(img) {
        img = ee.Image(img);
        var before = ee.ImageCollection.fromImages(ee.List(img.get('before'))).mosaic();
        var after  = ee.ImageCollection.fromImages(ee.List(img.get('after'))).mosaic();
        img = img.set('before', null).set('after', null);
        //替换空值,确认线性Interp有结果
        before = replace_mask(before, after, nodata);
        after  = replace_mask(after , before, nodata);
        //计算图像的时间之间的比率
        var x1 = before.select('time').double();
        var x2 = after.select('time').double();
        var now = ee.Image.constant(img.date().millis()).double();
        //待插影像时间减去前面的合成时间,除以后面的合成影像时间和前面的影像时间的差
        var ratio = now.subtract(x1).divide(x2.subtract(x1));  //x1 = x2时,为零
        //计算插值图像
        before = before.select(0); //移除时间波段
        after  = after.select(0);
        img    = img.select(0); 
        var interp = after.subtract(before).multiply(ratio).add(before).toFloat();
        var qc = img.mask().not().rename('qc');
        interp = replace_mask(img, interp, nodata).rename(bandname);
        //return interp.addBands(qc).copyProperties(img, img.propertyNames());
        return interp.copyProperties(img, img.propertyNames());
    }));
    return interpolated;
}

这段代码实现了对图像集合进行线性插值,生成时间连续的图像序列,从而对时间序列进行进一步分析和处理。下面是代码的详细介绍:

  1. linearInterp(imgcol, frame, nodata, bandname):该函数接受四个参数,分别为:
  • imgcol:图像集合,即待插值的图像序列。
  • frame:时间窗口,即选择多少天内的图像插值,默认值为32。
  • nodata:缺失值,即需要插值的像素点的默认值,默认值为0。
  • bandname:插值后的波段名称,默认值为原始波段名称。
  1. addTimeBand(img):该函数用于添加时间属性,即为图像添加一个名为“time”的新波段,该波段的值为图像的起始时间。
  2. ee.Filter.maxDifference(frame * (1000*60*60*24), time, null, time):该语句用于生成一个时间过滤器,用于对图像集合进行分组,以便进行插值操作。frame * (1000*60*60*24) 表示时间差异的阈值,单位为毫秒。
  3. ee.Join.saveAll({matchesKey:'after', ordering:time, ascending:false}).apply(imgcol, imgcol, f1):该语句用于将图像集合按照时间降序排序,并为每个图像找到它之后的最近的图像。在该语句中,matchesKey 表示保存关联图像的键值,ordering 表示排序关键字,ascending 表示是否升序排序。
  4. ee.Join.saveAll({matchesKey:'before', ordering:time, ascending:true}).apply(c1, imgcol, f2):该语句用于将图像集合按照时间升序排序,并为每个图像找到它之前的最近的图像。在该语句中,matchesKeyorderingascending 的含义与上述相同。
  5. interpolated:该语句是一个图像集合,其中的每个图像都是通过线性插值生成的新图像。具体来说,对于两个相邻的图像 beforeafter,首先用 replace_mask 函数替换缺失值,并计算出它们之间的时间比率,然后根据时间比率,通过线性插值生成一张新图像 interp。最后,将新图像的波段名称重命名为 bandname,并复制原始图像的属性。

综上所述,该代码实现了对图像序列进行线性插值,并返回一系列时间连续的图像序列,为后续的时间序列分析提供了重要的基础。

二、掩模替换函数

//掩模替换               
function replace_mask(img, newimg, nodata) {
    nodata   = nodata || 0;
    var mask = img.mask();
    img = img.unmask(nodata);
    img = img.where(mask.not(), newimg);
    img = img.updateMask(img.neq(nodata));
    return img;
}

这段代码实现的是掩模替换的功能。该函数的输入参数包括一个原始图像img、一个用于替换无效数据的新图像newimg和一个可选的无效数据值nodata。输出为经过掩模替换后的图像。

  • 函数中的第一行nodata = nodata || 0表示如果未传入无效数据值nodata,则默认值为0。
  • 接下来,该函数通过img.mask()获取了原始图像的掩模(mask),然后通过img.unmask(nodata)将原始图像的无效数据(即值等于nodata的像素)转换为NaN,以便后续处理。
  • 然后,函数调用img.where(mask.not(), newimg),该函数会将掩模中值为false的像素替换为新图像newimg中对应像素的值。这个操作就是所谓的掩模替换。
  • 最后,函数调用img.updateMask(img.neq(nodata)),该函数会将原始图像中值不等于nodata的像素重新赋上掩模,即原始图像中的有效数据区域。

最终函数返回经过掩模替换和更新掩模后的图像。

三、添加时间属性的函数

//添加时间属性
function addTimeBand(img) {
    /** 确保掩模一致 */
    var mask = img.mask();
    var time = img.metadata('system:time_start').rename("time").mask(mask);
    return img.addBands(time);
}

这段代码实现的是为图像添加时间属性的功能。该函数的输入参数为一个图像img,输出为添加时间属性后的图像。

  • 函数中的第一步操作是获取图像的掩模,以确保添加的时间属性的掩模和原始图像一致。具体而言,函数通过img.mask()获取了原始图像的掩模(mask)。
  • 接下来,函数通过img.metadata('system:time_start')获取原始图像的起始时间,即该图像的时间属性。获取的时间属性被重命名为time,然后通过.mask(mask)将其掩模与原始图像的掩模一致。这样就保证了时间属性的掩模和原始图像的掩模一致。
  • 最后,函数通过img.addBands(time)将时间属性作为一维带(band)添加到原始图像中,从而得到了一个带有时间属性的新图像。

这个函数通常在处理遥感影像时使用,由于遥感影像一般都包含时间属性,因此为遥感影像添加时间属性可以更好地进行时间序列分析,例如通过时间属性对影像进行分类、监测变化等。

四、提取质量控制 (QA) 位函数

//返回提取的QA位的单波段图像
var getQABits = function(image, start, end) {
    var pattern = 0;
    for (var i = start; i <= end; i++) {
      pattern += 1 << i;
    }
    return image.bitwiseAnd(pattern).rightShift(start);
};

这段代码实现的是从多波段图像中提取QA(Quality Assurance)位的单波段图像的功能。该函数的输入参数包括一个多波段图像image,一个表示QA位的起始位置start,和一个表示QA位的终止位置end。输出为提取出来的QA位的单波段图像。

函数中的核心操作是通过位运算提取QA位。

  • 首先,函数定义了一个变量pattern,并将其初始化为0。然后,函数通过循环遍历从startend的所有位数,将每个位的二进制表示中的1向左移动相应的位数,最后将它们相加,得到一个表示QA位的二进制模式。具体而言,这个过程就是将一个二进制的"1"左移若干位(位数由i决定),然后与变量pattern相加。
  • 接下来,函数通过image.bitwiseAnd(pattern)将二进制模式应用到输入的多波段图像image上,得到一个新的多波段图像,其中仅保留了QA位上的像素值。
  • 最后,函数通过.rightShift(start)将提取出来的QA位向右移动start位,得到一个单波段图像,即为提取出来的QA位的单波段图像。

该函数通常在遥感影像处理中使用,例如对Landsat影像的QA信息进行处理,以确定哪些像素是有用的(如云、阴影、水等)并排除掉无用的像素。

五、应用案例

1.《GEE:线性插值方法填补去云空洞——以 Landsat 和 NDVI 为例》

2.《GEE:线性插值方法填补去云空洞——以 MCD43A4 和 NDVI 时间序列为例》