【GEE】Sentinel-2 最新去云方法总结

第一种方法:使用QA波段去云

这是我们最常用的方法,具体原理就是利用QA60波段标记实现去云,具体代码如下:

代码语言:javascript
复制
var s2 = ee.ImageCollection("COPERNICUS/S2"),
point = /* color: #98ff00 */ee.Geometry.Point([116.20553071765003, 39.404020061278715]);

/**
*

  • remove cloud by QA bands
  • */

function rmCloudByQA(image) {

var qa = image.select('QA60');

var cloudBitMask = 1 << 10;

var cirrusBitMask = 1 << 11;

var mask = qa.bitwiseAnd(cloudBitMask).eq(0)

           .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

return image.updateMask(mask);

}

function main() {

var startDate = "2019-4-4";

var endDate = "2019-4-8";

Map.centerObject(point, 8);

var s2Imgs = s2.filterDate(startDate, endDate)

             .filterBounds(point);

Map.addLayer(s2Imgs.first(), {min:0, max:3000, bands:["B4", "B3", "B2"]}, "raw", false);

s2Imgs = s2Imgs.map(rmCloudByQA);

Map.addLayer(s2Imgs.first(), {min:0, max:3000, bands:["B4", "B3", "B2"]}, "cloud");

}
main();

代码分析:

这段代码中去云的方法是:rmCloudByQA()其他的代码也都是基础代码,我就不在过多解释。

第二种方法:计算云量分数去云

这种方法比较有意思,是一个国外的大牛利用波段组合以及NDSI指数计算得到当前像素是不是云,具体代码如下:

代码语言:javascript
复制
var s2 = ee.ImageCollection("COPERNICUS/S2"),

point = /* color: #98ff00 */ee.Geometry.Point([116.20553071765003, 39.404020061278715]);

var _cloudScore = function(img) {

var rescale = function(img, exp, thresholds) {

return img.expression(exp, {img: img})

    .subtract(thresholds[0]).divide(thresholds[1] - thresholds[0]);

};

var score = ee.Image.constant(1.0);

score = score.min(rescale(img, 'img.blue', [0.1, 0.3]));

score = score.min(rescale(img, 'img.red + img.green + img.blue', [0.2, 0.8]));

score = score.min(rescale(img, 'img.nir + img.swir1 + img.swir2', [0.3, 0.8]));

var ndsi = img.normalizedDifference(['green', 'swir1']);

return score.min(rescale(ndsi, 'img', [0.8, 0.6]));

};

/**
*

  • remove cloud by Score bands
  • */

function rmCloudByScore(image, thread) {

var preBands = ["B2","B3","B4","B8","B11","B12"];

var newBands = ['blue','green','red','nir','swir1','swir2'];

var score = _cloudScore(image.select(preBands, newBands));

score = score.multiply(100).byte().rename('cloud');

return image.addBands(score)

          .updateMask(score.lte(thread));

}

/**
*

  • scale image
  • */

function scaleImage(image) {

var time_start = image.get("system:time_start");

image = image.divide(10000);

image = image.set("system:time_start", time_start);

return image;

}

function main() {

var startDate = "2019-4-4";

var endDate = "2019-4-8";

Map.centerObject(point, 8);

var s2Imgs = s2.filterDate(startDate, endDate)

             .filterBounds(point)

             .map(scaleImage)

             .map(function(image) {

                return rmCloudByScore(image, 30);

             });

Map.addLayer(s2Imgs.first(), {min:0, max:0.3, bands:["B4", "B3", "B2"]}, "cloud");

}

main();

代码分析:

观察这段代码,有一点需要注意的是:调用去云方法rmCloudByScore(image, threshold)需要传入两个参数,一个是影像另外一个是识别云的阈值(也就是大于多大的值就认为是云,范围是0-100)。还有就是调用这个方法之间首先需要将影像的波段做缩放和将波段名称重新命名为指定的名称,这是因为这个算法内部主要逻辑方法_cloudScore(xxx)它内部计算时候取的都是特殊值,具体参考上述代码。

第三种方法:使用s2cloudless算法去云

这是目前最新的一种去云方法,s2cloudless算法是一种机器学习算法,Sentinel-2影像通过它可以计算得到云掩膜文件。这个算法是由SENTINEL Hub开源的一种算法,地址是:

GitHub - sentinel-hub/sentinel2-cloud-detector: Sentinel Hub Cloud Detector for Sentinel-2 images in Python

在GEE上我们不需要自己去用这个算法来计算云掩膜,GEE官方目前正在利用这个算法来生成一个新的数据集,这个数据集就是Sentinel-2云掩膜数据集(“COPERNICUS/S2_CLOUD_PROBABILITY”)。需要注意的是由于目前这个数据集依然在生产中,所以不一定所有的地方和所有时间段都存在,我们在使用的时候一定要注意。

下面以具体的例子说明如何使用这个数据集:

代码语言:javascript
复制
var s2 = ee.ImageCollection("COPERNICUS/S2"),

s2_cloud = ee.ImageCollection(&#34;COPERNICUS/S2_CLOUD_PROBABILITY&#34;),

point = /* color: #98ff00 */ee.Geometry.Point([116.20553071765003, 39.404020061278715]);

/**
*

  • remove cloud by probability bands
  • */

function rmCloudByProbability(image, thread) {

var prob = image.select("probability");

return image.updateMask(prob.lte(thread));

}

/**
*

  • scale image
  • */

function scaleImage(image) {

var time_start = image.get("system:time_start");

image = image.divide(10000);

image = image.set("system:time_start", time_start);

return image;

}

function getMergeImages(primary, secondary) {

var join = ee.Join.inner();

var filter = ee.Filter.equals({

leftField: &#39;system:index&#39;,

rightField: &#39;system:index&#39;

});

var joinCol = join.apply(primary, secondary, filter);

joinCol = joinCol.map(function(image) {

var img1 = ee.Image(image.get(&#34;primary&#34;));

var img2 = ee.Image(image.get(&#34;secondary&#34;));

return img1.addBands(img2);

});

return ee.ImageCollection(joinCol);

}
function main() {

var startDate = "2019-4-4";

var endDate = "2019-4-8";

Map.centerObject(point, 8);

var s2Imgs1 = s2.filterDate(startDate, endDate)

              .filterBounds(point)

              .map(scaleImage);

var s2Imgs2 = s2_cloud.filterDate(startDate, endDate)

                    .filterBounds(point);

var s2Imgs = getMergeImages(s2Imgs1, s2Imgs2);

s2Imgs = s2Imgs.map(function(image) {

return rmCloudByProbability(image, 50);

});

Map.addLayer(s2Imgs.first(), {min:0, max:0.3, bands:["B4", "B3", "B2"]}, "cloud");

}

main();

代码分析:

从上述代码可以看到,由于这是两个数据集,所以我们首先需要将筛选后的两个数据(s2Imgs1和s2Imgs2)利用join合并为一个数据集,具体方法参考getMergeImages(xxx)。合并后的影像中每一幅都包含一个叫做”probability“的波段,这个波段计算的就是云掩膜的信息(值范围是0-100)。程序中我们利用rmCloudByProbability(image, threshold)方法传入影像和云量阈值,就可以实现去云操作。