GEE, Sentinel 2 cloud and shadow

GEE在2021年8月30日左右，發(fā)布了Sentinel 2 去云的新教程，利用Probability產(chǎn)品和 CDI指數(shù)來去云。官方教程代碼如下：
https://code.earthengine.google.com/a7d6d6defee0ae0d0fdfe4d4c5011306
或者在Examples中也可以找到

image.png

以下內(nèi)容是自己的筆記

目錄

  • 知識(shí)點(diǎn)
  • 官方代碼講解
  • 實(shí)際應(yīng)用

一、知識(shí)點(diǎn)

  • CDI
    CDI是David Frantz在2018年創(chuàng)造的指數(shù)，這個(gè)指數(shù)是用來探測(cè)Sentinel 2中的云，并且還可以區(qū)分高亮云和建筑物，在GEE中的調(diào)用是ee.Algorithms.Sentinel2.CDI 。由于CDI指數(shù)是使用Sentinel 2 Level 1C 產(chǎn)品生成的，所以用這個(gè)算法計(jì)算CDI時(shí)，要用到Sentinel 2 Level 1C 。想要更加了解這個(gè)指數(shù)的可以看David Frantz在RSE上發(fā)表的文章Improvement of the Fmask algorithm for Sentinel-2 images: Separating clouds from bright surfaces based on parallax effects

  • ee.Join函數(shù)
    關(guān)于這部分內(nèi)容，推薦去看知乎上的一個(gè)作者，想要了解Join函數(shù)就必須先了解Filter，更加推薦大家去看視頻版，也就12分鐘。
    第10節(jié) GEE的參數(shù)類型 (Filter，Join) - 知乎 (zhihu.com)

二、官方代碼講解

1. 調(diào)用三個(gè)數(shù)據(jù)集

// 調(diào)用 Sentinel 2 L1C產(chǎn)品，來計(jì)算CDI
var s2 = ee.ImageCollection('COPERNICUS/S2');
// 調(diào)用 Probability產(chǎn)品
var s2c = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY');
// 調(diào)用 Sentine 2 L2A產(chǎn)品
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR');

2. 篩選

// ROI
var roi = ee.Geometry.Point([-122.4431, 37.7498]);
Map.centerObject(roi, 11);

// 定義時(shí)間范圍
var start = ee.Date('2019-03-01');
var end = ee.Date('2019-09-01');

// 選出CDI計(jì)算中要用到的波段
s2 = s2.filterBounds(roi).filterDate(start, end)
    .select(['B7', 'B8', 'B8A', 'B10']);

s2c = s2c.filterDate(start, end).filterBounds(roi);
// 這里對(duì)Sentinel 2 L2A產(chǎn)品的波段選擇沒有要求，也可以選擇其他波段，或者全選，但注意和L1C區(qū)分
s2Sr = s2Sr.filterDate(start, end).filterBounds(roi)
    .select(['B2', 'B3', 'B4', 'B5']);

3. 定義函數(shù)：將兩個(gè)數(shù)據(jù)集合并
下面這個(gè)函數(shù)，實(shí)際上是將collectionB 按照時(shí)間屬性，加入到collectionA中。比如：有2個(gè)數(shù)據(jù)集，一個(gè)是EVI，一個(gè)是NDVI，它們都有一年的時(shí)間序列，我想按照時(shí)間整合兩個(gè)數(shù)據(jù)集到一個(gè)數(shù)據(jù)集當(dāng)中，以便于分析。
注意：ee.Join.saveFirst 返回的是左邊數(shù)據(jù)集collectionA，只是將 collectionB當(dāng)作一個(gè)屬性添加到collectionA中的屬性當(dāng)中而已

// propertyName 是一個(gè)MatchKey，由用戶自定義，可以是任何字符串
//為了后面的map函數(shù)，ee.ImageCollection 一定要加，因?yàn)閑e.Join函數(shù)返回的是Join對(duì)象
function indexJoin(collectionA, collectionB, propertyName) {
  var joined = ee.ImageCollection(ee.Join.saveFirst(propertyName).apply({
    primary: collectionA,
    secondary: collectionB,
    condition: ee.Filter.equals({
      leftField: 'system:index',
      rightField: 'system:index'})
  }));
  // 通過get函數(shù)獲取右邊數(shù)據(jù)集的影像，再通過addBands合并成一個(gè)ImageCol
  return joined.map(function(image) {
    return image.addBands(ee.Image(image.get(propertyName)));
  });
}

4. 定義函數(shù)：去云
這一步才是真正的開始去云以及去掉陰影，前面都是數(shù)據(jù)集的準(zhǔn)備

function maskImage(image) {
  // 計(jì)算CDI，下面的B10是Sentinel2 L1C的波段
  var cdi = ee.Algorithms.Sentinel2.CDI(image);
  var s2c = image.select('probability');
  var cirrus = image.select('B10').multiply(0.0001);

 // 閾值設(shè)定，滿足下列條件的都是云，這里的閾值是官方給的
  var isCloud = s2c.gt(65).and(cdi.lt(-0.5)).or(cirrus.gt(0.01));

  // 以下代碼似乎是通過面積來判斷，有知道的可以評(píng)論一下
  // Reproject is required to perform spatial operations at 20m scale.
  // 20m scale is for speed, and assumes clouds don't require 10m precision.
  isCloud = isCloud.focal_min(3).focal_max(16);
  isCloud = isCloud.reproject({crs: cdi.projection(), scale: 20});

  // Project shadows from clouds we found in the last step. This assumes we're working in
  // a UTM projection.
  var shadowAzimuth = ee.Number(90)
      .subtract(ee.Number(image.get('MEAN_SOLAR_AZIMUTH_ANGLE')));

  // With the following reproject, the shadows are projected 5km.
  isCloud = isCloud.directionalDistanceTransform(shadowAzimuth, 50);
  isCloud = isCloud.reproject({crs: cdi.projection(), scale: 100});

  isCloud = isCloud.select('distance').mask();
  return image.select('B2', 'B3', 'B4').updateMask(isCloud.not());
}

5. 合并與去云
官方這里給的是無(wú)云中值合成案例

// 將Probability 加入到 Sentinel2 L2A中
var withCloudProbability = indexJoin(s2Sr, s2c, 'cloud_probability');

// 將Sentinel2 L1C 加入到 withCloudProbability 中
var withS2L1C = indexJoin(withCloudProbability, s2, 'l1c');

// 對(duì)最后合并的數(shù)據(jù)集進(jìn)行去云
var masked = ee.ImageCollection(withS2L1C.map(maskImage));

// 中值合成，下面設(shè)置8是為了避免memory errors
var median = masked.reduce(ee.Reducer.median(), 8);

三、實(shí)際應(yīng)用

生成云掩膜的過程實(shí)際上只用到L1C產(chǎn)品和Probability產(chǎn)品，如果按照官方的步驟，把L2A產(chǎn)品和Probability產(chǎn)品先合并，那么L2A產(chǎn)品中的波段不能含有 ['B7', 'B8', 'B8A', 'B10']，因?yàn)楹竺孢€會(huì)和L1C產(chǎn)品合并，CDI只用到了L1C產(chǎn)品中的 ['B7', 'B8', 'B8A', 'B10']來計(jì)算。如果將L2A和L1C的['B7', 'B8', 'B8A', 'B10']都合并在一起,會(huì)出現(xiàn)波段命名不一樣的情況，如下圖所示，此時(shí)如果用maskImage,CDI計(jì)算用到的波段可能是L2A的。

image.png

而我們?cè)谑褂眠^程中，B8這個(gè)波段對(duì)我們很重要，不能丟掉B8。因此，可以先生成掩膜文件，再通過Map循環(huán)來掩膜。
代碼如下：

// 數(shù)據(jù)集
var s2 = ee.ImageCollection('COPERNICUS/S2');
var s2c = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY');
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR');

var roi = ee.Geometry.Point([-122.4431, 37.7498]);
Map.centerObject(roi, 11);

var start = ee.Date('2020-01-01');
var end = ee.Date('2021-01-01');

//篩選
s2 = s2.filterBounds(roi).filterDate(start, end)
    .select(['B7', 'B8', 'B8A', 'B10']);
s2c = s2c.filterDate(start, end).filterBounds(roi);

var s2Sr = s2Sr
            .filterBounds(roi)
            .filterDate(start,end)
            .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE',50))

// Join two collections on their 'system:index' property.
function indexJoin(collectionA, collectionB, propertyName) {
  var joined = ee.ImageCollection(ee.Join.saveFirst(propertyName).apply({
    primary: collectionA,
    secondary: collectionB,
    condition: ee.Filter.equals({
      leftField: 'system:index',
      rightField: 'system:index'})
  }));
  // Merge the bands of the joined image.
  return joined.map(function(image) {
    return image.addBands(ee.Image(image.get(propertyName)));
  });
}

// 生成一個(gè)云掩膜波段，并將這個(gè)掩膜波段命名為 CloudMask
function maskImage(image) {
  // Compute the cloud displacement index from the L1C bands.
  var cdi = ee.Algorithms.Sentinel2.CDI(image);
  var prob = image.select('probability');
  var cirrus = image.select('B10').multiply(0.0001);
  var isCloud = prob.gt(65).and(cdi.lt(-0.5)).or(cirrus.gt(0.01));

  isCloud = isCloud.focal_min(3).focal_max(16);
  isCloud = isCloud.reproject({crs: cdi.projection(), scale: 20});

  var shadowAzimuth = ee.Number(90)
      .subtract(ee.Number(image.get('MEAN_SOLAR_AZIMUTH_ANGLE')));

  isCloud = isCloud.directionalDistanceTransform(shadowAzimuth, 50);
  isCloud = isCloud.reproject({crs: cdi.projection(), scale: 100});

  isCloud = isCloud.select('distance').mask().rename('CloudMask');
  return isCloud.not();
}

// 將probability產(chǎn)品合并到s2中
var withCloudProbability = indexJoin(s2, s2c, 'cloud_probability');
// 生成掩膜文件
var masked = ee.ImageCollection(withCloudProbability.map(maskImage));
//將掩膜文件合并到L2A產(chǎn)品中
var s2SrWithMask = indexJoin(s2Sr,masked,'could_mask');

// 使用Map函數(shù)，開始掩膜
var S2_cloudMask = s2SrWithMask.map(function(img){
    var mask = img.select('CloudMask')
    return img.updateMask(mask)
              .divide(10000)
              .select('B.*')
              .toFloat()
              .copyProperties(img, ["system:time_start"])
})

print('S2_cloudMask:',S2_cloudMask)
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容