GEE:应用遥感影像时空插值技术的实践(插值填补去云空洞)
作者:_养乐多_
本文介绍了几个 Google Earth Engine (GEE) 平台中常用的处理遥感数据中的缺失值代码片段,这些代码可以用于在时间序列中对遥感图像进行线性插值,提供更加连续和完整的时间序列。
第一段代码是一个线性插值函数,它能够在一个时间序列图像集合中对缺失数据进行插值。该函数使用了 GEE 平台中的图像集合操作函数和图像合成函数,并采用了线性插值方法对缺失值进行处理,最后返回一个插值后的图像集合。(该代码可以应用与去云后通过时间序列数据进行插值填补影像空洞)
第二段代码是一个掩模替换函数,它的作用是将原始影像中的某些值(如缺失值)用一个新值替换。它使用了 GEE 平台中的图像操作函数,并使用掩膜对影像进行处理,最后返回替换后的影像。
第三段代码是一个添加时间属性的函数,它会在影像中添加一个时间波段,用于后续时间序列分析。该函数同样使用了 GEE 平台中的图像操作函数和元数据属性获取函数。
最后一段代码是一个提取质量控制 (QA) 位的函数,它能够从一个多波段图像中提取指定位置的 QA 位数据。它使用了位运算函数和移位函数,并使用了 GEE 平台中的图像操作函数。
目录
一、线性插值函数
二、掩模替换函数
三、添加时间属性的函数
四、提取质量控制 (QA) 位函数
一、线性插值函数
//线性插值填补无云像素
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 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);// 待插影像之后的所有影像// 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);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;
}
这段代码实现了对图像集合进行线性插值,生成时间连续的图像序列,从而对时间序列进行进一步分析和处理。下面是代码的详细介绍:
linearInterp(imgcol, frame, nodata, bandname)
:该函数接受四个参数,分别为:
imgcol
:图像集合,即待插值的图像序列。frame
:时间窗口,即选择多少天内的图像插值,默认值为32。nodata
:缺失值,即需要插值的像素点的默认值,默认值为0。bandname
:插值后的波段名称,默认值为原始波段名称。
-
addTimeBand(img)
:该函数用于添加时间属性,即为图像添加一个名为“time”的新波段,该波段的值为图像的起始时间。 -
ee.Filter.maxDifference(frame * (1000*60*60*24), time, null, time)
:该语句用于生成一个时间过滤器,用于对图像集合进行分组,以便进行插值操作。frame * (1000*60*60*24)
表示时间差异的阈值,单位为毫秒。 -
ee.Join.saveAll({matchesKey:'after', ordering:time, ascending:false}).apply(imgcol, imgcol, f1)
:该语句用于将图像集合按照时间降序排序,并为每个图像找到它之后的最近的图像。在该语句中,matchesKey
表示保存关联图像的键值,ordering
表示排序关键字,ascending
表示是否升序排序。 -
ee.Join.saveAll({matchesKey:'before', ordering:time, ascending:true}).apply(c1, imgcol, f2)
:该语句用于将图像集合按照时间升序排序,并为每个图像找到它之前的最近的图像。在该语句中,matchesKey
、ordering
、ascending
的含义与上述相同。 -
interpolated
:该语句是一个图像集合,其中的每个图像都是通过线性插值生成的新图像。具体来说,对于两个相邻的图像before
和after
,首先用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。然后,函数通过循环遍历从start
到end
的所有位数,将每个位的二进制表示中的1向左移动相应的位数,最后将它们相加,得到一个表示QA位的二进制模式。具体而言,这个过程就是将一个二进制的"1"左移若干位(位数由i
决定),然后与变量pattern
相加。
接下来,函数通过image.bitwiseAnd(pattern)
将二进制模式应用到输入的多波段图像image
上,得到一个新的多波段图像,其中仅保留了QA位上的像素值。最后,函数通过.rightShift(start)
将提取出来的QA位向右移动start
位,得到一个单波段图像,即为提取出来的QA位的单波段图像。
该函数通常在遥感影像处理中使用,例如对Landsat影像的QA信息进行处理,以确定哪些像素是有用的(如云、阴影、水等)并排除掉无用的像素。