GEE在2021年8月30日左右甘邀，發(fā)布了Sentinel 2 去云的新教程琅坡，利用Probability產(chǎn)品和 CDI指數(shù)來(lái)去云逗物。官方教程代碼如下：
https://code.earthengine.google.com/a7d6d6defee0ae0d0fdfe4d4c5011306
或者在Examples中也可以找到
以下內(nèi)容是自己的筆記
目錄
- 知識(shí)點(diǎn)
- 官方代碼講解
- 實(shí)際應(yīng)用
一跌前、知識(shí)點(diǎn)
CDI
CDI是David Frantz在2018年創(chuàng)造的指數(shù)钮糖，這個(gè)指數(shù)是用來(lái)探測(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 effectsee.Join函數(shù)
關(guān)于這部分內(nèi)容泊业，推薦去看知乎上的一個(gè)作者，想要了解Join函數(shù)就必須先了解Filter断傲，更加推薦大家去看視頻版脱吱，也就12分鐘。
第10節(jié) GEE的參數(shù)類(lèi)型 (Filter认罩，Join) - 知乎 (zhihu.com)
二、官方代碼講解
1. 調(diào)用三個(gè)數(shù)據(jù)集
// 調(diào)用 Sentinel 2 L1C產(chǎn)品续捂，來(lái)計(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)品的波段選擇沒(méi)有要求垦垂，也可以選擇其他波段，或者全選牙瓢，但注意和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，由用戶(hù)自定義弓候，可以是任何字符串
//為了后面的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'})
}));
// 通過(guò)get函數(shù)獲取右邊數(shù)據(jù)集的影像菇存，再通過(guò)addBands合并成一個(gè)ImageCol
return joined.map(function(image) {
return image.addBands(ee.Image(image.get(propertyName)));
});
}
4. 定義函數(shù)：去云
這一步才是真正的開(kāi)始去云以及去掉陰影夸研，前面都是數(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è)定依鸥，滿(mǎn)足下列條件的都是云亥至，這里的閾值是官方給的
var isCloud = s2c.gt(65).and(cdi.lt(-0.5)).or(cirrus.gt(0.01));
// 以下代碼似乎是通過(guò)面積來(lái)判斷，有知道的可以評(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)用
生成云掩膜的過(guò)程實(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']來(lái)計(jì)算睡榆。如果將L2A和L1C的['B7', 'B8', 'B8A', 'B10']都合并在一起,會(huì)出現(xiàn)波段命名不一樣的情況萍肆，如下圖所示，此時(shí)如果用maskImage,CDI計(jì)算用到的波段可能是L2A的胀屿。
而我們?cè)谑褂眠^(guò)程中塘揣，B8這個(gè)波段對(duì)我們很重要，不能丟掉B8宿崭。因此亲铡，可以先生成掩膜文件，再通過(guò)Map循環(huán)來(lái)掩膜葡兑。
代碼如下：
// 數(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ù)，開(kāi)始掩膜
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)
