【GEE】以某个日期为中心镶嵌合成研究区域完整影像

---

前言

一些新手GEE学习记录
本文介绍如何在GEE中以某个日期为中心，向前向后找寻影像镶嵌合成研究区域完整影像，在镶嵌过程中优先显示更靠近该日期的影像像素值。
本文没有实现自动向前向后搜索影像，而是手动修改前后时间范围看合成效果不断扩大或缩小，有兴趣者可以自己写个循环实现。

---

一、影像导入和预处理

1. 预处理函数

//Apply scale and offset factors for Sentinel-2.
//https://blog.csdn.net/m0_71243797/article/details/132417873
exports.applyScaleFactors_S2 = function(image) {
  image = ee.Image(image);
  var opticalBands = image.select('B.').divide(10000);
  return image.addBands(opticalBands, null, true);
};
 
 
// S2去云、去雪
exports.maskS2 = function(image) {
  var scl = image.select('SCL');
  var Cloud_Shadows = 1 << 3;
  var Clouds_Low_Probability = 1 << 7;
  var	Clouds_Medium_Probability = 1 << 8;
  var Clouds_High_Probability = 1 << 9;
  var Cirrus = 1 << 10;
  var	Snow_Ice = 1 << 11;
  
  var qa = image.select('QA60');
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  
  var mask = scl.bitwiseAnd(Cloud_Shadows).eq(0)
    .and(scl.bitwiseAnd(Clouds_Low_Probability).eq(0))
    .and(scl.bitwiseAnd(Clouds_Medium_Probability).eq(0))
    .and(scl.bitwiseAnd(Clouds_High_Probability).eq(0))
    .and(scl.bitwiseAnd(Cirrus).eq(0))
    .and(scl.bitwiseAnd(Snow_Ice).eq(0))
    .and(qa.bitwiseAnd(cloudBitMask).eq(0))
    .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  
  return image.updateMask(mask);
};

// Reverse imagecollection's order
exports.reverseImageColl = function(imageColl){
  //for restore new collection(after reverse)
  var imageColl_new = imageColl;
  //imagecollection to list
  var list= imageColl.toList(100);
  var list_reverse = list.reverse();
  //list to imagecollection
  imageColl_new = ee.ImageCollection.fromImages(list_reverse);
  return imageColl_new;
};

2. 影像导入和预处理

//Load functions
//影像预处理的函数保存在此脚本中，使用时需引入
var pref = require('users/yueeashion/Data_search:Preprocessing_function.js');

//Define research period
//此例以20190901为中心日期，向前向后使用了各半个月的影像合成
var start_1 = '2019-08-15';
var end_1 = '2019-08-31';
var start_2 = '2023-09-01';
var end_2 = '2023-09-15';
//Load region ROI
var table = ee.FeatureCollection("users/yueeashion/A/LosRios_area");
var roi = table; 

//Plot ROI Boundry
var my_style_1={color:'yellow',fillColor:'00000000'};
var my_style_2={color:'red',fillColor:'red'};
Map.addLayer(roi.style(my_style_1),{}, "LosRios");
Map.addLayer(roi_points.style(my_style_2),{}, "Balsa location");


//Filter date and area
var my_filter_1 = ee.Filter.and(
  ee.Filter.date(start_1, end_1),
  ee.Filter.bounds(roi)
);
var my_filter_2 = ee.Filter.and(
  ee.Filter.date(start_2, end_2),
  ee.Filter.bounds(roi)
);
//Define show style of images
var visParams = {bands: 'red,green,blue', min: -0.01, max: 0.2};

//Sentinel-2
// First period images collection(before 09.01)
//First period 和second period影像处理步骤完全相同
var s2_1 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
  .filter(my_filter_1)
  .map(pref.maskS2)  //delete clouds and snow area
  .map(pref.applyScaleFactors_S2)  //Apply scale and offset   

// 识别去云错误影像
  .map(function(image){
    var image_mean = ee.Image(image)
    .select(['B2', 'B3', 'B4', 'B8'])//这里的波段似乎不能乱改
    .reduce(ee.Reducer.mean());
    return image.addBands(image_mean.rename("ref_mean"));
  })
  .map(function(image){
    var mean = image.select('ref_mean').reduceRegion({
      reducer: ee.Reducer.mean(),
      scale: 90,
      });
    return image.set("ref_mean", mean.get('ref_mean'));
  })

  // 波段选择和重命名，波段按自己需要选择保留
  .select(['B2', 'B3', 'B4','B5','B6','B7', 'B8','B8A', 'B11', 'B12'])
  .map(function(image){
    return ee.Image(image).rename(['blue', 'green', 'red', 'nir1','nir2','nir3','nir4','nir5', 'swir1', 'swir2']);
  });
  
// paint ref_mean条形图, for set the threshold
var refMeanList = s2_1.aggregate_array('ref_mean');
var chart = ui.Chart.array.values(refMeanList, 0)
  .setChartType('ColumnChart')
  .setOptions({
    title: '属性 ref_mean 的条形图',
    hAxis: {title: 'ref_mean'},
    vAxis: {title: '影像序号'},
    colors: ['red']
  });
print(chart);

//删除去云错误影像
s2_1 = s2_1.filter(ee.Filter.lt("ref_mean", 0.2));



// Second period images collection(after 09.01)
var s2_2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
  .filter(my_filter_2)
  .map(pref.maskS2)  //delete clouds and snow area
  .map(pref.applyScaleFactors_S2)  //Apply scale and offset   
// 识别去云错误影像
  .map(function(image){
    var image_mean = ee.Image(image)
    .select(['B2', 'B3', 'B4', 'B8'])
    .reduce(ee.Reducer.mean());
    return image.addBands(image_mean.rename("ref_mean"));
  })
  .map(function(image){
    var mean = image.select('ref_mean').reduceRegion({
      reducer: ee.Reducer.mean(),
      scale: 90,
      });
    return image.set("ref_mean", mean.get('ref_mean'));
  })

  // 波段选择和重命名
  .select(['B2', 'B3', 'B4','B5','B6','B7', 'B8','B8A', 'B11', 'B12'])
  .map(function(image){
    return ee.Image(image).rename(['blue', 'green', 'red', 'nir1','nir2','nir3','nir4','nir5', 'swir1', 'swir2']);
  });
//删除去云错误影像
s2_2 = s2_2.filter(ee.Filter.lt("ref_mean", 0.2));

3. 镶嵌合成

s2_1和s2_2的影像集合中，影像都是按拍摄时间的顺序排列的。mosaic()函数的规则是，将该imageCollection的最后一张影像放在最top的位置，越靠前的影像越在底层，我们若想让9.1日前后的影像都遵循“越靠近9.1日的越优先显示，即越在顶层”，s2_1顺序不变，s2_2要倒转影像顺序。

//---------------------------------Sentinel-2 合成图 --------------------
Define show style of images
var visParams = {bands: 'red,green,blue', min: -0.01, max: 0.2};

//Mosaic function makes the last image on the top, we want to make images which are closer to date 9.1 on the top
//thus s2_1 keep on, s2_2: Reverse its images' order
var s2_2_order = pref.reverseImageColl(s2_2);

//Mosaic two collections
var s2_1_mosaic = s2_1.mosaic();
var s2_2_mosaic = s2_2_order.mosaic();
// Show them
Map.addLayer(s2_1_mosaic.clip(roi), visParams, 's2_1_mosaic');
Map.addLayer(s2_2_mosaic.clip(roi), visParams, 's2_2_mosaic');

//!!!Before this step，test to see which collection is better and put the better one at last
var s2_mosaic = ee.ImageCollection([s2_2_mosaic,s2_1_mosaic]).mosaic().clip(roi);//put the better one at last
print("s2_mosaic",s2_mosaic);
Map.addLayer(s2_mosaic, visParams, 's2_mosaic');

二、像素覆盖度计算

计算当前合成影像占研究区roi的像素比例。
可以用于在此基础上继续写一个循环自动向前向后寻找影像时，终止的条件阈值计算，如像素覆盖度 > 99%时终止。

参考来自：https://blog.csdn.net/weixin_40694662/article/details/124823133【GEE笔记】有效像元（面积、数量）统计

//Pixel proportion
var all_pixel = ee.Image.pixelArea().reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: roi,
    scale: 10,
    maxPixels: 10e16,
}).get("area");
print("all_pixel：", all_pixel);
var real_pixel = ee.Image.pixelArea()
.updateMask(s2_mosaic.select(1).mask())  //select(1)这里其实是只用了s2_mosaic的第2个波段计算非null像素数
.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: roi,
    scale: 10,
    maxPixels: 10e16,
}).get("area");

//Conunt no null pixels in each band
var eachband_pixel = s2_mosaic.reduceRegion({
  reducer: ee.Reducer.count(),
  scale: 10,
  crs: 'EPSG:4326',
  geometry: roi,
  maxPixels: 1e16
});

print("real_pixel:",real_pixel);
print("像素比例", ee.Number(real_pixel).divide(all_pixel));//只用了一个波段的非null像素值除以roi区域内像素总数，得到的是一个数值
print('each band pixels:', eachband_pixel);  //输出每个波段的非null像素数