Код SLC не заполняет все пробелы Landsat 7 SR

Question

Я новичок в коде Google Earth Engine и пытаюсь применить код SLC-gap к изображениям Landsat 7 Surface Reflectance. Используя ресурсы, доступные в StackOverflow, я сгенерировал код ниже; однако, когда я передаю изображения в QGIS, все еще остаются пробелы. Мой код неверен или я неправильно применил его к изображениям?

Сначала я замаскировал облака на основе полосы pixel_qa данных Landsat SR:

var cloudMaskL7 = function(image) {
var qa = image.select('pixel_qa');
var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));

Затем я удалил краевые пиксели, которые встречаются не во всех полосах:

var mask2 = image.mask().reduce(ee.Reducer.min());
return image.updateMask(cloud.not()).updateMask(mask2);
 };
var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
              .filterDate('2004-09-15', '2004-12-31')
              .map(cloudMaskL7);

var visParams = {
bands: ['B3', 'B2', 'B1'],
min: 0,
max: 3000,
gamma: 1.4,
};
Map.setCenter(36.197, 31.701,7);
Map.addLayer(l7.median(), visParams);

Затем я нанёс на карту функцию за один год данных TOA Landsat 7, взял медиану и отобразил ее для Иордании.

var composite = l7.map(cloudMaskL7)
    .median();
Map.setCenter(36.124, 31.663);
Map.addLayer(composite, {bands: ['B4', 'B3', 'B2'], max: 0.3});

Затем я попытался заполнить пробелы в SLC Landsat 7, применив протокол заполнения промежутка Фаза-2 USGS L7, используя один размер ядра.

var MIN_SCALE = 1/3;
var MAX_SCALE = 3;
var MIN_NEIGHBORS = 144;
var GapFill = function(src, fill, kernelSize) {
var kernel = ee.Kernel.square(kernelSize * 30, 'meters', false);
var common = src.mask().and(fill.mask());
var fc = fill.updateMask(common);
var sc = src.updateMask(common);

Затем я нашел основные коэффициенты масштабирования с помощью регрессии и чередовал полосы для регрессии (предполагается, что полосы имеют одинаковые имена).

   var regress = fc.addBands(sc);
    regress = regress.select(regress.bandNames().sort());
    var fit = regress.reduceNeighborhood(ee.Reducer.linearFit().forEach(src.bandNames()),  kernel, null, false);
    var offset = fit.select('.*_offset');
    var scale = fit.select('.*_scale');

Затем я нашел вторичные коэффициенты масштабирования, используя просто средства и стандартное отклонение.

 var reducer = ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true);
    var src_stats = src.reduceNeighborhood(reducer, kernel, null, false);
    var fill_stats = fill.reduceNeighborhood(reducer, kernel, null, false);
    var scale2 = src_stats.select('.*stdDev').divide(fill_stats.select('.*stdDev'));
    var offset2 = src_stats.select('.*mean').subtract(fill_stats.select('.*mean').multiply(scale2));

    var invalid = scale.lt(MIN_SCALE).or(scale.gt(MAX_SCALE));
    scale = scale.where(invalid, scale2);
    offset = offset.where(invalid, offset2);

Я применил масштабирование и маскировку пикселей, у которых не было достаточно соседей.

var count = common.reduceNeighborhood(ee.Reducer.count(), kernel, null, true, 'boxcar');
var scaled = fill.multiply(scale).add(offset)
      .updateMask(count.gte(MIN_NEIGHBORS));

  return src.unmask(scaled, true);
};

var source = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var fill = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');

Я загрузил таблицу границ и фильтр.

var Jordan = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
    .filter(ee.Filter.or(
        ee.Filter.eq('country_co', 'JO')));
        var clippedJordan = composite.clipToCollection(Jordan);

Я отобразил результаты для Джордана; однако пробелы в SLC, по-видимому, не заполняются. Я продолжаю вычислять значения MSAVI2, используя эти изображения, поэтому оставшиеся пробелы влияют на результаты.

var mc = Map.setCenter(36.274, 31.682, 6);
var visParams = {bands: ['B3', 'B2', 'B1']};
Map.addLayer(clippedJordan, visParams, 'clipped composite');

Любой совет будет с благодарностью!

Answer 1

Я не знаком с новейшей технологией Landsat 7, заполняющей пробелы для изображений SLC-off, но вот значительно упрощенная версия того, что вы пытались сделать.Я удалил много (ненужных?) Вещей, сильно увеличил размер ядра и увеличил период времени, в течение которого генерируется медианная замена.Это может быть близко к тому, что вы хотите:

var cloudMaskL7 = function(image) {
  var qa = image.select('pixel_qa');
  var cloud = qa.bitwiseAnd(1 << 5)
                    .and(qa.bitwiseAnd(1 << 7))
                    .or(qa.bitwiseAnd(1 << 3));
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
};

var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
    .map(cloudMaskL7);

var kernelSize = 10;
var kernel = ee.Kernel.square(kernelSize * 30, 'meters', false);

var GapFill = function(image) {
  var start = image.date().advance(-1, 'year');
  var end = image.date().advance(1, 'year');
  var fill = l7.filterDate(start, end).median();
  var regress = fill.addBands(image); 
  regress = regress.select(regress.bandNames().sort());
  var fit = regress.reduceNeighborhood(ee.Reducer.linearFit().forEach(image.bandNames()), kernel, null, false);
  var offset = fit.select('.*_offset');
  var scale = fit.select('.*_scale');
  var scaled = fill.multiply(scale).add(offset);
  return image.unmask(scaled, true);
};

// TESTING CODE
var point = ee.Geometry.Point(36.124, 31.663);
Map.centerObject(point, 11);

var check = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
    .filterBounds(point)
    .filterDate('2004-09-15', '2004-12-31');
var checkImage = ee.Image(check.first());
var visParams = {bands: ['B4', 'B3', 'B2'], min: 200, max: 5500};
Map.addLayer(checkImage, visParams, 'source');

// Test composite.
var checkStart = checkImage.date().advance(-1, 'year');
var checkEnd = checkImage.date().advance(1, 'year');
var composite = l7.filterDate(checkStart, checkEnd).median();
Map.addLayer(composite, visParams, 'median');

// Rough implementation for comparison.
var replaced = checkImage.unmask(composite);
Map.addLayer(replaced, visParams, 'simple');

// Fancy implementation.
var filled = ee.Image(check.map(GapFill).first());
Map.addLayer(filled, visParams, 'filled');

РЕДАКТИРОВАТЬ: Ответ теперь показывает, как отобразить это на коллекции.Берегись, потому что я не знаю, насколько хорошо это будет масштабироваться.Если вы решите нанести на карту большую область или длинный временной ряд, вас предупредят.