在GEE使用线性插值和平均插值插补缺少的影像

GEE缺失影像插补平滑完整方案1：使用线性平均插值

KeyPoints

• 遥感数据缺失和平滑是一个问题。
• 数据缺失存在两种情况，一种是局部数据缺失，一种是整张影像缺失，影像后续分析。
• 数据平滑有助于去除噪声，使得数据更为清晰、准确。此外，遥感数据中常常包含着地形、云层等导致数据的波动较大，数据平滑可以改善图像质量。
• 对于数据处理，采取“先插补，再平滑”的策略。本期推文解决插补问题，下次解决平滑问题。
• 本文利用线性平均和历史平均插值处理两种数据缺失的情况。

线性插值

首先需要导入线性插值函数库，可以参考https://github.com/gee-hydro/gee_docs

var pkg_smooth = require('users/kongdd/public:Math/pkg_smooth.js');

线性插值函数pkg_smooth.linearInterp用法如下：

usage: pkg_smooth.linearInterp(imgcol, frame, nodata)

- imgcol        The ImageCollection to be interpolated. 
- frame         The window size (unit: day) to seach the nearenearest valid good values
                before and after the current point. Then used the nearest valid good
                values to linterpolate missing value.
- nodata        The value of `NODATA`. The default is 0. Caution about band's data 
                type.

提供一个实例：

var imgcol_interp = pkg_smooth.linearInterp(imgcol_input, 50);

即对imgcol_input进行线性插值，考虑前后50天的情况。

历史平均插值

在介绍pkg_smooth.historyInterp历史平均插值之前，先介绍pkg_trend.aggregate_prop函数，该函数类似R中的aggregate。

想象一种情况，在使用MODIS 8-day NDVI产品时，我们想计算每个月平均NDVI值，就可以用这种语句进行。

var NDVI_month    = pkg_trend.aggregate_prop(NDVI_d8.select(0), 'month'   , 'mean')

但在这之前需要对NDVI的ImageCollection每张图像添加month属性，也有代码很方便地进行实现：

var NDVI_d8 = NDVI_d8.map(pkg_trend.add_dn(false, 8))

效果如图，会添加一系列时间的属性，方便后续aggregate

image-20230610193610817

用这两个函数还能把小时尺度降水合成日尺度等一系列操作：

    var dataset = ee.ImageCollection('NASA/GPM_L3/IMERG_V05')
    var startdate = ee.Date.fromYMD(2014,3,1)
    var enddate = ee.Date.fromYMD(2014,4,1)
    var precipitation = dataset.filter(ee.Filter.date(startdate,enddate)).select('precipitationCal')
    print(precipitation)
    
    var difdate = enddate.difference(startdate, 'day')
    
    // Time lapse
    var lapse = ee.List.sequence(0, difdate.subtract(1))
    var startdate = ee.Date('2014-01-01')
    var listdates = lapse.map(function(day){
      return startdate.advance(day, 'day')
    })
    
    var pts = ee.FeatureCollection(ee.List([]))
    var newft = ee.FeatureCollection(listdates.iterate(function(img, ft) {
      // Cast
      ft = ee.FeatureCollection(ft)
      var day = ee.Date(img)
      // Filter the collection in one day
      var day_collection = precipitation.filterDate(day, day.advance(1, 'day'))
      // Get the sum of all 24 images into one Image
      var sum = ee.Image(day_collection.sum())
      // Return the FeatureCollection with the new properties set
      return sum
    }, listdates))

https://stackoverflow.com/questions/56306009/aggregate-gpm-hourly-data-to-daily-in-gee

历史平均插值有几种情况：

• 利用其他年份的同一个dn值进行插值，如每一年的这一天平均值，最精确
• 利用其他年份的同一个月的平均值进行插值，较为精确
• 利用其它年份的年平均值进行插值，最不精确

下面是三种历史平均插值代码：

var imgcol_interp = pkg_smooth.linearInterp(Emiss_raw, nday);
var imgcol_hisavg_d8    = pkg_trend.aggregate_prop(Emiss_d8.select(0), 'dn'   , 'median')
var imgcol_hisavg_month = pkg_trend.aggregate_prop(Emiss_d8.select(0), 'month', 'median')
var imgcol_hisavg_year  = pkg_trend.aggregate_prop(Emiss_d8.select(0), 'year' , 'median');

var imgcol_his_d8 = pkg_smooth.historyInterp(imgcol_interp, imgcol_hisavg_d8   , 'dn')
var imgcol_his_1m = pkg_smooth.historyInterp(imgcol_his_d8, imgcol_hisavg_month, 'month');
var imgcol_his_1y = pkg_smooth.historyInterp(imgcol_his_1m, imgcol_hisavg_year , 'year');

先通过aggregate_prop计算平均值（dn，月或年），再通过historyInterp函数插值，第一个参数是待插值ImgCol，中间的参数是计算出的平均值，最后参数是在哪个尺度下进行插值。

两种数据缺失的情况

接下来展示两种数据缺失的情况，并提供解决思路

数据局部缺失

如图所示，MODIS的某年反照率数据，可能由于云层或其他原因，有大面积缺失

image-20230610190036519

先线性插值，再历史平均插值，代码如下，用dn插值最精确：

var imgcol_interp = pkg_smooth.linearInterp(Emiss_raw, nday);
var imgcol_hisavg_d8    = pkg_trend.aggregate_prop(Emiss_d8.select(0), 'dn'   , 'median')
var imgcol_hisavg_month = pkg_trend.aggregate_prop(Emiss_d8.select(0), 'month', 'median')
var imgcol_hisavg_year  = pkg_trend.aggregate_prop(Emiss_d8.select(0), 'year' , 'median');

var imgcol_his_d8 = pkg_smooth.historyInterp(imgcol_interp, imgcol_hisavg_d8   , 'dn')

效果如图：

整张影像缺失

一些遥感数据会出现整张影像缺失的情况，如下图的情况，整个6月是没有影像的。这种情况不能直接插补，需要用其他方法。

image-20230610202508181

以GRACE为例，首先需要手动找到缺失影像的日期，建立一个列表。

根据日期列表建立空影像，再把数据补入进来，按时间排序：

var missing_tiles = ee.List(['2003-06-01','2004-01-01', '2011-01-01', '2011-06-01', 
'2011-11-01', '2012-04-01', '2012-05-01', '2012-10-01', 
'2013-03-01', '2013-08-01', '2013-09-01', '2014-02-01',
'2014-07-01', '2014-12-01', '2015-05-01', '2015-07-01',
'2015-10-01', '2015-11-01', '2016-02-01', '2016-04-01',
'2016-09-01', '2016-10-01'])

/** Select GRACE data ----------------------------------- */
var GRACE_month  = ee.ImageCollection('NASA/GRACE/MASS_GRIDS/LAND')
            .filter(ee.Filter.date('2003-01-01', '2017-01-01'))
            .select('lwe_thickness_csr').map(pkg_trend.add_dn(false)).map(function(img){return img.toFloat()});
print(GRACE_month, 'GRACE_ori')

        
var prj = GRACE_month.first().projection();    
var GRACE_miss = ee.ImageCollection(missing_tiles.map(function(date){

  date = ee.Date.parse('YYYY-MM-dd', date)
  var id = ee.String('NASA/GRACE/MASS_GRIDS/LAND/').cat(date.format('YYYYMMdd')).cat('_').cat(date.format('YYYYMMdd')).cat('_linear')
  var img = ee.Image().select(0).rename('lwe_thickness_csr').reproject(prj)
  return(img.set('system:time_start', date.millis())
            .set('system:id', id)
            .set('system:index', date.format('YYYY-MM-dd'))
  )
  
}))
  
GRACE_month = ee.ImageCollection(GRACE_month.merge(GRACE_miss))
GRACE_month = GRACE_month.sort('system:time_start')
print(GRACE_month, 'GRACE_month+miss')

对合并的ImageCollection进行线性插值

var imgcol_interp = pkg_smooth.linearInterp(imgcol_input, 50);

结果如图，补全了缺失的6月影像

image-20230610203021233

完整代码后台回复【GEE插补】

本文利用线性平均和历史平均插值处理两种数据缺失的情况，下期利用Savitzky-Golay滤波和Whittaker滤波再对缺失的影像进行平滑并分享代码。

https://code.earthengine.google.com/d23dcc273921e4f78a5aba2e2fd1f663

本文由 mdnice 多平台发布