基于Landsat和MODIS进行多波段融合（反射率、地表温度）

主要任务：

之前帮朋友做一个数据缺失插补的任务，她遇到的主要问题是：landsat去云后，缺失数据比较多，希望用MODIS数据进行填补，更准确的术语是：数据融合（data fusion）。具体来说是用，MCD43A4的反射率、地表温度去插补Landsat 8 的反射率和地表温度，示例的时间选为2017年。

基本原理

下面介绍的是，Nietupski，2021的工作 （Redirecting）。MODIS的时间分辨率高，但是空间分辨率低，Landsat可以捕捉更精细的空间细节，但时间分辨率低。在过去十年中，时空自适应反射融合模型 STARFM被不断开发，用于生成高精度高频的遥感产品。STARFM无法适应空间上突然变化下的融合，为此有一些改进的算法，如ESTARFM（陈晋老师近几年的工作）。Nietupski在这篇文章中，开发了一种能够在GEE上实现类似ESTARFM功能的方法，并主要创新在用于估计物候。

基本理论：两张Landsat图像一般间隔16天，中间缺失的用MODIS进行融合（插补）。在对MODIS进行重采样到30m后，逐像素线性插值Landsat缺失的日期图像，同时考虑临近的空间信息和光谱信息（这部分没太看懂，友友们可以看看原作）。

由于Nietupski，2021 的代码是python API写得（可以在github上搜"Nietupski"），我转为了在 GEE上可运行的代码。

结果：

注意landsat的数据是分path和row的，所以我是根据某个point的经纬度，获取了整块path和row的图像。我后续还在做一个多块Landsat并行融合的代码，bug还是有些，期待友友们贡献智慧，一起交流。

融合前

融合后

[Landsat融合后，数据太大，导致画图超内存。下图展示的是完整融合后的预期结果，空白处被填充，但是不是真实结果，下次学习python API后可能可以解决超内存的问题]

代码

/*-----------                        --------------
                 global variables
-------------                       --------------*/
// Define the region of interest
var point = ee.Geometry.Point([35.93, 15.5]);
var region = point.buffer(250);

// Define temporal bounds
var startDate = '2017-5-01';
var endDate = '2017-6-30';
// Define common band names for fusion
var commonBandNames = ['ndvi', 'blue'];
// Image collections for fusion
var landsatCollection = 'LANDSAT/LC08/C02/T1_L2';
var modisCollection = 'MODIS/006/MCD43A4';
// Band names and indices for Landsat
var bandNamesLandsat = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'QA_PIXEL'];
var landsatBands = [1, 2, 3, 4, 5, 6, 17];
// Band names and indices for MODIS
var bandNamesModis = ['red', 'nir', 'blue', 'green', 'swir1', 'swir2'];
var modisBands = [0, 1, 2, 3, 4, 5];
// Kernel configuration for moving window operations,# radius of moving window
// Note: Generally, larger windows are better but as the window size increases,
// so does the memory requirement and we quickly will surpass the memory
// capacity of a single node (in testing 13 was max size for single band, and 10 was max size for up to 6 bands)
var kernelRadius = 10;
var kernel = ee.Kernel.square(kernelRadius);
var numPixels = Math.pow(kernelRadius * 2 + 1, 2);
// Number of land cover classes
var coverClasses = 7;
// Asset export path
//var path = 'projects/ee-zxy/assets/';


/*-----------                        --------------
                 Functions group 1
-------------                       --------------*/
// Define the function to mask Landsat images
function maskLandsat(image) {
    // Bits 3 and 5 are cloud shadow and cloud, respectively. 4 is snow
    var cloudShadowBitMask = 1 << 3;
    var cloudsBitMask = 1 << 5;
    var snowBitMask = 1 << 4;
    // Get the pixel QA band
    var qa = image.select('QA_PIXEL');
    // Create the mask
    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
        .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
        .and(qa.bitwiseAnd(snowBitMask).eq(0));
    // Apply the mask
    var maskedImage = image.updateMask(mask);
    // Count the number of non-masked pixels
    var maskedCount = maskedImage.select(['QA_PIXEL'])
        .reduceRegion({
            reducer: ee.Reducer.count(),
            geometry: image.geometry(),
            scale: 30,
            maxPixels: 1e10
        });
    // Count the total number of pixels
    var origCount = image.select(['blue'])
        .reduceRegion({
            reducer: ee.Reducer.count(),
            geometry: image.geometry(),
            scale: 30,
            maxPixels: 1e10
        });
    // Calculate the percent of masked pixels
    var percent = ee.Number(origCount.get('blue'))
        .subtract(maskedCount.get('QA_PIXEL'))
        .divide(origCount.get('blue'))
        .multiply(100)
        .round();
    // Return the masked image with new property and timestamp
    return maskedImage
        .set('CloudSnowMaskedPercent', percent)
        .copyProperties(image, ['system:time_start']);
}

// Define the function to mask MODIS images
function maskMODIS(image) {
    // Calculate the snow water index
    var swi = image.expression(
        '(green * (nir - swir1)) / ((green + nir) * (nir + swir1))', {
            'green': image.select('green'),
            'nir': image.select('nir'),
            'swir1': image.select('swir1')
        }).rename('swi');
    // Mask out values of swi above 0.1
    var mask = swi.lt(0.1);
    return image
        .updateMask(mask)
        .copyProperties(image, ['system:time_start', 'system:id']);
}

// Define the function to add NDVI to images
function addNDVI(image) {
    var ndvi = image.normalizedDifference(['nir', 'red']).rename('ndvi');
    return image.addBands(ndvi);
}

// Define the function to get paired Landsat and MODIS images
function getPaired(startDate, endDate, landsatCollection, landsatBands, bandNamesLandsat, modisCollection, modisBands, bandNamesModis, commonBandNames, region) {
    var landsat = ee.ImageCollection(landsatCollection)
        .filterDate(startDate, endDate)
        .filterBounds(region)
        .select(landsatBands, bandNamesLandsat)
        .map(addNDVI)
        .map(maskLandsat)
        .map(function (image) {
            return image.setMulti({
                'system:time_start': ee.Date(image.date()).millis(),
                'DOY': image.date().format('D'),
                'date':ee.Date(image.date()).format('yyyy-MM-dd')
            });
        })
        .select(commonBandNames);
    //print(landsat)
    
    var modis = ee.ImageCollection(modisCollection)
        .filterDate(startDate, endDate)
        .select(modisBands, bandNamesModis)
        .map(addNDVI)
        .map(maskMODIS)
        .map(function (image) {
            return image.set({'DOY': image.date().format('D'),
            'date':ee.Date(image.date()).format('yyyy-MM-dd')});
        })
        .select(commonBandNames);
    //print(modis)
    
    // Filter the two collections by the date property
    var dayfilter = ee.Filter.equals('date', null, 'date');

    // Define simple join
    var pairedJoin = ee.Join.simple();
    // Define inverted join to find MODIS images without Landsat pair
    var invertedJoin = ee.Join.inverted();

    // Create collections of paired Landsat and MODIS images
    var landsatPaired = pairedJoin.apply(landsat, modis, dayfilter);
    var modisPaired = pairedJoin.apply(modis, landsat, dayfilter);
    var modisUnpaired = invertedJoin.apply(modis, landsat, dayfilter);

    return [landsatPaired, modisPaired, modisUnpaired];
}

// Define the function to get dates from images and append to list
function getDates(image, empty_list) {
    // Get date and update format
    var date = ee.Image(image).date().format('yyyy-MM-dd');

    // Add date to 'empty list'
    var updatelist = ee.List(empty_list).add(date);

    return updatelist;
}

// Define the function to make subcollections from paired images
function makeSubcollections(paired) {
    // Empty list to store dates
    var empty_list = ee.List([]);

    // Fill empty list with dates
    var dateList = paired[0].iterate(getDates, empty_list);

    // Function to create individual subcollection
    function getSub(ind) {
        // Get Landsat images
        var lan_01 = paired[0]
            .filterDate(ee.List(dateList).get(ind),
                ee.Date(ee.List(dateList).get(ee.Number(ind).add(1)))
                    .advance(1, 'day'))
            .toList(2);

        // Get MODIS paired images
        var mod_01 = paired[1]
            .filterDate(ee.List(dateList).get(ind),
                ee.Date(ee.List(dateList).get(ee.Number(ind).add(1)))
                    .advance(1, 'day'))
            .toList(2);

        // Get MODIS images between these two dates
        var mod_p = paired[2]
            .filterDate(ee.List(dateList).get(ind),
                ee.Date(ee.List(dateList).get(ee.Number(ind).add(1)))
                    .advance(1, 'day'));

        mod_p = mod_p.toList(mod_p.size());

        // Combine collections to one object
        var subcollection = ee.List([lan_01, mod_01, mod_p]);

        return subcollection;
    }

    // Filter out sub collections from paired and unpaired collections
    var subcols = ee.List.sequence(0, ee.List(dateList).length().subtract(2))
        .map(getSub);

    return subcols;
}

/*-----------                        --------------
                 Functions group 2 prepare
-------------                       --------------*/

// 定义函数来对MODIS图像进行重采样和配准到Landsat参考图像
function registerImages(landsat_t01, modis_t01, modis_tp) {
    // 使用bicubic插值对图像进行重采样
    landsat_t01 = landsat_t01.map(function(image) {
        return ee.Image(image).resample('bicubic');
    });

    modis_t01 = modis_t01.map(function(image) {
        return ee.Image(image).resample('bicubic');
    });

    modis_tp = modis_tp.map(function(image) {
        return ee.Image(image).resample('bicubic');
    });

    // 将MODIS图像配准到Landsat t0图像
    modis_t01 = modis_t01.map(function(image) {
        return ee.Image(image).register({
            referenceImage: ee.Image(landsat_t01.get(0)),
            maxOffset: 150.0,
            stiffness: 7.0
        });
    });

    modis_tp = modis_tp.map(function(image) {
        return ee.Image(image).register({
            referenceImage: ee.Image(landsat_t01.get(0)),
            maxOffset: 150.0,
            stiffness: 7.0
        });
    });

    return [landsat_t01, modis_t01, modis_tp];
}

// 定义函数来确定每个Landsat图像的相似性阈值
function threshold(landsat, coverClasses) {
    var getThresh = function(image) {
        // 计算图像中每个波段的标准差
        var stddev = ee.Image(image).reduceRegion({
            reducer: ee.Reducer.stdDev(),
            bestEffort: true,
            maxPixels: 1e6
        });

        // 将标准差字典转换为多波段图像
        stddev = ee.Image(stddev.toImage());

        // 获取波段名称并重命名以用于阈值
        var names = stddev.bandNames().map(function(bn) {
            return ee.String(bn).cat('_thresh');
        });

        // 根据标准差和地物分类数计算阈值
        var thresh = stddev.multiply(ee.Image.constant(2).divide(coverClasses)).rename(names);

        return thresh;
    };

    // 对每个Landsat图像计算阈值
    var threshs = ee.List(landsat).map(getThresh);

    return threshs;
}

// 使用计算出的阈值掩蔽不相似的像素
function threshMask(neighLandsat_t01, thresh, commonBandNames) {
    var masks = ee.List([0, 1]).map(function(i) {
        return commonBandNames.map(function(name) {
            // 获取t0或t1邻域图像，计算与中心像素的距离
            var dist = ee.Image(neighLandsat_t01.get(i))
                .select([ee.String(name).cat('_(.+)')])
                .select([ee.String(name).cat('_0_0')])
                .subtract(ee.Image(neighLandsat_t01.get(i)).select([ee.String(name).cat('_(.+)')]))
                .abs();

            // 掩蔽超过阈值的像素
            return dist.lte(ee.Image(thresh.get(i)).select([ee.String(name).cat('_(.+)')]));
        });
    });

    return masks;
}

// 准备MODIS图像以便进行分析
function prepMODIS(modis_t01, modis_tp, kernel, numPixels, commonBandNames, pixelBandNames) {
    // 转换图像为邻域图像
    var neighMod_t01 = modis_t01.map(function(image) {
        return ee.Image(image).neighborhoodToBands(kernel);
    });

    var neighMod_tp = modis_tp.map(function(image) {
        return ee.Image(image).neighborhoodToBands(kernel);
    });

    // 转换为数组图像
    var modArr = ee.Image(neighMod_t01.get(0)).toArray()
        .arrayCat(ee.Image(neighMod_t01.get(1)).toArray(), 0);

    // 创建按像素位置切片的数组列表
    var modPixArrays_t01 = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
        return modArr.arraySlice(0, ee.Number(i).int(),
           ee.Number(numPixels).multiply(ee.List(commonBandNames).length().multiply(2)).int(), numPixels);
    });

    var modPixArrays_tp = neighMod_tp.map(function(image) {
        return ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
            return ee.Image(image).toArray().arraySlice(0, ee.Number(i).int(),
                ee.Number(numPixels).multiply(ee.List(commonBandNames).length()).int(), numPixels);
        });
    });

    // 展平数组并根据doy、波段和像素位置命名
    var modSorted_t01 = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
        return ee.Image(modPixArrays_t01.get(i)).arrayFlatten([pixelBandNames.get(i)]);
    });

    var modSorted_tp = ee.List.sequence(0, modPixArrays_tp.length().subtract(1)).map(function(i) {
        return ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(x) {
            return ee.Image(ee.List(modPixArrays_tp.get(i)).get(x)).arrayFlatten([commonBandNames])
                .set('DOY', ee.Image(modis_tp.get(i)).get('DOY'));
        });
    });

    return [modSorted_t01, modSorted_tp];
}

// 准备Landsat图像以便进行分析
function prepLandsat(landsat_t01, kernel, numPixels, commonBandNames, doys, coverClasses) {
    // 转换图像为邻域图像
    var neighLandsat_t01 = landsat_t01.map(function(image) {
        return ee.Image(image).neighborhoodToBands(kernel);
    });

    // 创建像素位置列表
    var pixPositions = ee.Image(neighLandsat_t01.get(0)).bandNames().map(function(bn) {
        return ee.String(bn).replace('[a-z]+_', '_');
    }).slice(0, numPixels);

    // 创建波段名称列表以重命名输出数组
    var pixelBandNames = pixPositions.map(function(position) {
        return doys.map(function(doy) {
            return commonBandNames.map(function(bn) {
                return ee.String(doy).cat('_').cat(ee.String(bn)).cat(ee.String(position));
            });
        }).flatten();
    });

    // 转换为数组并排序
    var lanArr_t01 = ee.Image(neighLandsat_t01.get(0)).toArray()
        .arrayCat(ee.Image(neighLandsat_t01.get(1)).toArray(), 0);

    var pixArrays = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
        return lanArr_t01.arraySlice(0, ee.Number(i).int(),
            ee.Number(numPixels).multiply(ee.List(commonBandNames).length().multiply(2)).int(), numPixels);
    });

    // 展平数组并根据doy、波段和像素位置命名
    var lanSorted = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
        return ee.Image(pixArrays.get(i)).arrayFlatten([pixelBandNames.get(i)]);
    });

    // 为图像确定阈值
    var thresh = threshold(landsat_t01, coverClasses);

    // 使用阈值掩蔽窗口图像
    var mask_t01 = threshMask(neighLandsat_t01, thresh, commonBandNames);

    // 将掩蔽的列表转换为图像然后转为数组
    var maskArr_t01 = ee.ImageCollection(ee.List(mask_t01.flatten())).toBands().toArray();

    var maskArrays = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
        return maskArr_t01.arraySlice(0, ee.Number(i).int(),
            ee.Number(numPixels).multiply(ee.List(commonBandNames).length().multiply(2)).int(), numPixels);
    });

    // 展平掩蔽数组并根据doy、波段和像素位置命名
    var masksSorted = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(i) {
        return ee.Image(maskArrays.get(i)).arrayFlatten([pixelBandNames.get(i)]);
    });

    // 掩蔽Landsat图像
    var maskedLandsat = ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(index) {
        return ee.Image(lanSorted.get(index)).updateMask(ee.Image(masksSorted.get(index)));
    });

    return [maskedLandsat, pixPositions, pixelBandNames];
}

/*-----------                        --------------
                 Functions group 3 fusion
-------------                       --------------*/
// 计算Landsat和MODIS像素之间的平均绝对光谱距离
function calcSpecDist(maskedLandsat, modSorted_t01, numPixels, pixPositions) {
    var sDist = ee.List.sequence(0, ee.Number(numPixels).subtract(1))
        .map(function(index) {
            return ee.Image(maskedLandsat.get(index))
                .subtract(ee.Image(modSorted_t01.get(index)))
                .abs()
                .reduce(ee.Reducer.mean());
        });

    return ee.ImageCollection(sDist)
        .toBands()
        .rename(pixPositions.map(function(name) {
            return ee.String('sDist').cat(name);
        }));
}

// 计算窗口中每个像素与中心像素之间的空间距离
function calcSpatDist(positions) {
    // 窗口宽度
    var w2 = positions.length().sqrt().subtract(1).divide(2);

    // 窗口中每个像素的距离
    var dist = positions.map(function(position) {
        return ee.Image.constant(1)
            .add(ee.Number.parse(ee.String(position).match('(-?[0-9]+)', 'g').get(0)).pow(2))
            .add(ee.Number.parse(ee.String(position).match('(-?[0-9]+)', 'g').get(1)).pow(2))
            .sqrt()
            .divide(w2);
    });

    return ee.ImageCollection(dist)
        .toBands()
        .rename(positions.map(function(bn) {
            return ee.String('corr').cat(bn);
        }));
}

// 从空间和光谱距离图像创建对角权重矩阵
function calcWeight(spatDist, specDist) {
    var disIndex = specDist.multiply(spatDist);
    var num = ee.Image.constant(1).divide(disIndex);
    var sumNum = ee.ImageCollection(num.bandNames().map(function(bn) {
        return num.select([bn]).rename('num');
    })).sum();
    var W = num.divide(sumNum);

    return W.unmask().toArray().toArray(1).matrixToDiag();
}

// 计算线性回归的转换系数
function calcConversionCoeff(maskedLandsat, modSorted_t01, doys, numPixels, commonBandNames) {
    // 重新格式化Landsat和MODIS数据
    var lanMod = doys.map(function(doy) {
        return ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(index) {
            return ee.Image.constant(1).rename(['intercept'])
                .addBands(ee.Image(modSorted_t01.get(index))
                    .select(ee.String(doy).cat('.+'))
                    .rename(commonBandNames.map(function(bn) {
                        return ee.String(bn).cat('_modis');
                    })))
                .addBands(ee.Image(maskedLandsat.get(index))
                    .select(ee.String(doy).cat('.+'))
                    .rename(commonBandNames.map(function(bn) {
                        return ee.String(bn).cat('_landsat');
                    })));
        });
    });

    // 将集合转换为数组并解决线性回归系数
    lanMod = ee.ImageCollection(lanMod.flatten());
    var coeffs = lanMod.reduce(ee.Reducer.linearRegression(ee.List(commonBandNames).length().add(1), ee.List(commonBandNames).length()))
        .select([0], ['coefficients'])
        .arraySlice(0, 1, ee.List(commonBandNames).length().add(1));

    return coeffs;
}

// 使用权重和系数从MODIS预测Landsat
function predictLandsat(landsat_t01, modSorted_t01, doys, modSorted_tp, weights, coeffs, commonBandNames, numPixels) {
    // 根据doy分离MODIS图像
    var modSplit_t01 = doys.map(function(doy) {
        return modSorted_t01.map(function(image) {
            return ee.Image(image).select(ee.String(doy).cat('.+')).rename(commonBandNames);
        });
    });

    // 计算t0和t1的MODIS与预测时间的MODIS的差异
    var diffMod = modSplit_t01.map(function(imagelist) {
        return ee.ImageCollection(ee.List.sequence(0, ee.Number(numPixels).subtract(1)).map(function(index) {
            return ee.Image(modSorted_tp.get(index)).subtract(ee.Image(ee.List(imagelist).get(index)));
        })).toArray();
    });

    // 基于系数和MODIS反射差异计算y_hat，并进行加权
    var sumPixel = diffMod.map(function(arrayImage) {
        return weights.matrixMultiply(ee.Image(arrayImage).matrixMultiply(coeffs))
            .arrayReduce(ee.Reducer.sum(), [0])
            .arrayProject([1])
            .arrayFlatten([ee.List(commonBandNames)]);
    });

    // 创建基于MODIS反射差异的权重
    var sumDiff_t0 = ee.Image.constant(1).divide(ee.Image(diffMod.get(0)).arrayReduce(ee.Reducer.sum(), [0]).abs());
    var sumDiff_t1 = ee.Image.constant(1).divide(ee.Image(diffMod.get(1)).arrayReduce(ee.Reducer.sum(), [0]).abs());
    var weight_t0 = sumDiff_t0.divide(sumDiff_t0.add(sumDiff_t1)).arrayProject([1]).arrayFlatten([ee.List(commonBandNames)]);
    var weight_t1 = sumDiff_t1.divide(sumDiff_t0.add(sumDiff_t1)).arrayProject([1]).arrayFlatten([ee.List(commonBandNames)]);
    var temporalWeight = ee.List([weight_t0, weight_t1]);

    // 根据像素差异和原始图像预测新图像
    var predictions = ee.List([0, 1]).map(function(time) {
        return ee.Image(landsat_t01.get(time))
            .add(ee.Image(sumPixel.get(time)))
            .multiply(ee.Image(temporalWeight.get(time)));
    });

    var mergedPrediction = ee.Image(predictions.get(0)).add(ee.Image(predictions.get(1)));

    return mergedPrediction.setMulti({
        'DOY': ee.Image(modSorted_tp.get(0)).get('DOY'),
        'system:time_start': ee.Image(modSorted_tp.get(0)).get('system:time_start')
    });
}


/*-----------                        --------------
                 Main 
-------------                       --------------*/

// sorted, filtered, paired image retrieval
var pairedCollections = getPaired(startDate, endDate, landsatCollection, landsatBands, bandNamesLandsat, modisCollection, modisBands, bandNamesModis, commonBandNames, region);
print('pairedCollections:', pairedCollections);
var subs = makeSubcollections(pairedCollections);
print('Subcollections:', subs);


// 创建一个空的ImageCollection用于存储所有的预测结果
var allPredictions = ee.ImageCollection([]);

// 获取配对图像列表的数量
var numLists = subs.length();


// Process each pair by iterating through the sequence
ee.List.sequence(0, numLists.subtract(1)).getInfo().forEach(function(i){
  // Determine the number of MODIS images in the pair
  var sub = ee.List(subs.get(i));
  var numImgs = ee.List(sub.get(2)).length();
  //print('Number of MODIS images in pair', numImgs);

  // Determine the remaining number of images if not a multiple of 10
  var remaining = numImgs.mod(10);

  // Create the starting index sequence for MODIS images
  var indexSeq = ee.List.sequence(0, numImgs.subtract(remaining), 10);

  // Images for grouping and prediction
  var subList = ee.List(sub.get(2));

  // Process each batch using map (simulate forEach by processing side effects)
  var processBatch=indexSeq.getInfo().forEach(function(x){
    // Starting index
    var start = ee.Number(x);
    // Ending index
    var end = ee.Algorithms.If(start.add(10).gt(numImgs), numImgs, start.add(10));

    // Group of images to predict
    var predGroup = subList.slice(start, end);

    var landsatT01 = ee.List(sub.get(0));
    var modisT01 = ee.List(sub.get(1));
    var modisTp = predGroup;

    // Get starting and ending day of year values and year
    var startDay = ee.Number.parse(ee.ImageCollection(predGroup).first().get('DOY'));
    var endDay = ee.Number.parse(ee.ImageCollection(predGroup).sort('system:time_start', false).first().get('DOY'));
    var year = ee.Date(ee.ImageCollection(predGroup).sort('system:time_start', false).first().get('system:time_start')).format('Y');

    // Get DOY list
    var doys = landsatT01.map(function(img) {
      return ee.String(ee.Image(img).get('DOY')).cat('_');
    });

    // Register images
    var registeredImages = registerImages(landsatT01, modisT01, modisTp);
    landsatT01 = registeredImages[0];
    modisT01 = registeredImages[1];
    modisTp = registeredImages[2];

    // Prepare Landsat images (masking and formatting)
    var landsatPrep = prepLandsat(landsatT01, kernel, numPixels, commonBandNames, doys, coverClasses);
    var maskedLandsat = landsatPrep[0];
    var pixPositions = landsatPrep[1];
    var pixBN = landsatPrep[2];

    // Prepare MODIS images (masking and formatting)
    var modisPrep = prepMODIS(modisT01, modisTp, kernel, numPixels, commonBandNames, pixBN);
    var modSortedT01 = modisPrep[0];
    var modSortedTp = modisPrep[1];
    //print('modSortedTp', modSortedTp)
    
    // Calculate spectral distance
    var specDist = calcSpecDist(maskedLandsat, modSortedT01, numPixels, pixPositions);
    // Calculate spatial distance
    var spatDist = calcSpatDist(pixPositions);
    // Calculate weights
    var weights = calcWeight(spatDist, specDist);
    // Calculate conversion coefficients
    var coeffs = calcConversionCoeff(maskedLandsat, modSortedT01, doys, numPixels, commonBandNames);
    //print(weights,coeffs)
    
    // Predict MODIS images
    var predictions = ee.List(modSortedTp).map(function(image) {
      return predictLandsat(landsatT01, modSortedT01, doys, ee.List(image), weights, coeffs, commonBandNames, numPixels);
    });
    //print(predictions.getInfo())
    
    // Create new band names for multi-band NDVI image
    var dates = ee.List(modisTp).map(function(img) {
      return ee.Image(img).get('system:time_start');
    });

    var predNames = ee.List.sequence(0, predictions.length().subtract(1))
      .map(function(i) {
        return commonBandNames.map(function(name) {
          return ee.Date(dates.get(i)).format('yyyy-MM-dd').cat('-').cat(name);
        });
      }).flatten();
    //print(predNames)
    
    // Convert predictions to multi-band image and apply new band names
    var preds = ee.ImageCollection(predictions).toBands().rename(predNames);
    //print('multi-band image', preds)
    
    // Split multi-band image into single-band images and add to image collection
    var imageCollection = ee.ImageCollection(predNames.map(function(bandname) {
      var bandnameLen = ee.String(bandname).length()
      var date = ee.String(bandname).slice(0, 4).cat('-')
          .cat(ee.String(bandname).slice(5, 7)).cat('-')
          .cat(ee.String(bandname).slice(8, 10));
      return preds.select(ee.String(bandname)).rename(ee.String(bandname).slice(11, bandnameLen))
        .set('system:time_start', ee.Date(date).millis())
        .set('date', date);
    }));
    //print('one imagecollection, but date of band twice', imageCollection)
    
    // ImageCollection merged by date
    var mergedImageCollection = ee.ImageCollection(
      imageCollection.distinct('date').map(function(image) {
      var date = ee.String(image.get('date'));
      var imagesOnDate = imageCollection.filter(ee.Filter.equals('date', date));
      var mergedImage = imagesOnDate.toBands(); // 0_ndvi; 1_blue
      var bandNames = mergedImage.bandNames().map(function(name) {
        return ee.String(name).split('_').slice(1).join('_');
      });
      mergedImage = mergedImage.rename(bandNames);
      return mergedImage.set('system:time_start', image.get('system:time_start')).set('date', date);
    }));
    //print('ImageCollection,', mergedImageCollection);
    
    // merge to the null imagecollection
    allPredictions = allPredictions.merge(ee.ImageCollection(mergedImageCollection));
  });
});

// Output all predictions
print('All Predictions:', allPredictions);


// 创建一个新的列表，手动设置system:index; 直接在imagecollection上不能修改system:index
var imageList = allPredictions.toList(allPredictions.size());
var updatedImageList = imageList.map(function(img) {
  return ee.Image(img).set('system:index', ee.Image(img).get('date'));
});
// 将列表转换回ImageCollection
allPredictions = ee.ImageCollection(updatedImageList);
// 输出所有预测结果
print('allPredictions, after revising system:index', allPredictions);
print('dimension: ', allPredictions.first().select(0).getInfo()['bands'][0]['dimensions']) // 不清楚，为什么生成的imagecollection是7631*7831，也就是2.07*2.11°左右

/*-----------                        --------------
                 proprocess, draw, save 
-------------                       --------------*/
// before data fusion, draw ndvi, also for blue band 
var beforeNDVIImageCollection = ee.ImageCollection(pairedCollections[0]).select('ndvi');
print(beforeNDVIImageCollection);
var beforeNDVIChart = ui.Chart.image.series({
  imageCollection: beforeNDVIImageCollection,
  region: point,
  reducer: ee.Reducer.mean(),
  scale: 100,
  xProperty: 'date'
}).setOptions({
  title: 'NDVI before fusion',
  hAxis: {title: 'Date', format: 'yyyy-MM-dd', gridlines: {count: 10}},
  vAxis: {title: 'NDVI', minValue: 0, maxValue: 0.2},
  lineWidth: 2,
  pointSize: 4
});
print(beforeNDVIChart);

var beforeNDVIImage = ee.ImageCollection(pairedCollections[0]).select('ndvi').first();
print('beforeNDVI image', beforeNDVIImage);
Map.addLayer(beforeNDVIImage, {min: 0, max: 0.2, palette: ['blue', 'white', 'green']}, 'NDVIbef');

// after data fusion
var afterNDVIImageCollection = allPredictions.select('ndvi');
print(afterNDVIImageCollection);
var afterNDVIChart = ui.Chart.image.series({
  imageCollection: afterNDVIImageCollection,
  region: point,
  reducer: ee.Reducer.mean(),
  scale: 100,
  xProperty: 'date'
}).setOptions({
  title: 'NDVI after fusion',
  hAxis: {title: 'Date', format: 'yyyy-MM-dd', gridlines: {count: 10}},
  vAxis: {title: 'NDVI', minValue: 0, maxValue: 0.2},
  lineWidth: 2,
  pointSize: 4
});
print(afterNDVIChart);

var afterNDVIImage = allPredictions.select('ndvi').first();
print('afterNDVI image', afterNDVIImage)
Map.addLayer(afterNDVIImage, {min: 0, max: 0.2, palette: ['blue', 'white', 'green']}, 'NDVIaft');

Map.centerObject(point, 7);
Map.addLayer(ee.Image().paint(point.buffer(3000).bounds(), 1, 2).visualize({palette: ['blue']}), {}, 'point在的buffer');