GEE合成无云的Landsat-8和哨兵-2数据

GIS与遥感开发平台

发布于 2022-04-29 08:30:17

5.7K0

发布于 2022-04-29 08:30:17

Google Earth Engine（合成的影像）

Landsat-8（30m分辨率）和哨兵-2（10m分辨率）都是我们常用的遥感影像数据。

但是，遥感影像在很多情况下都会有云，我们就需要对影像进行去云处理，这样会使影像出现数据的缺失。我们可以利用Google earth engine对多景影像进行合成，得到无云的影像，而且也没有数据的缺失（类似于MODIS8天反射率产品）。

还有一种情况就是，我们的研究区域一景影像覆盖不了，我们通常都是通过下载影像利用ENVI镶嵌裁剪。而用Google earth engine，我们就可以直接下载镶嵌好的数据，十分的方便。

今天，我们就以Landsat-8和哨兵-2为例，演示如何合成无云的影像并且下载。（GEE数据下载方法上期已经讲过，就不再写在代码里了。）

哨兵-2数据合成

哨兵-2数据的空间分辨率为10m，时间分辨率为5天。

本次我们利用ee.Geometry.Rectangle函数制作一个矩形的区域来选择我们的哨兵-2的数据。

选择的数据为2018-06-01到2018-10-30之间的数据，我们把云量大于百分20的数据筛选掉。并且通过哨兵-2的质量检查波段进行去云处理。然后，对所有的数据进行平均值合成，就可以得到一景无云的哨兵数据。

代码：

var bound=ee.Geometry.Rectangle([115.3814, 41.0928,117.5020,  39.4227]);
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11是云，我们要把它mask掉
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
//哨兵的像元值是反射率的10000倍，要除以10000
  return image.updateMask(mask).divide(10000);
}
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2018-06-01', '2018-10-30')
                 .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .filterBounds(bound)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B8', 'B4', 'B3'],
};

Map.setCenter(116, 40, 12);
Map.addLayer(dataset.median(), rgbVis, 'RGB');

结果（GEE显示）：

本地显示（ENVI）：

Landsat-8合成无云的影像

下面我们来讲一下如果利用Landsat-8来进行合成无云的影像，并且我们还要计算一下Landsat-8的NDVI。然后对NDVI进行合成。最后我们选择了一个研究区域（以北京为例），对研究区域进行裁剪（需要自己在Assets里面上传北京的shp文件）。最后，我们就会得到北京区域在2017-06-01到2017-10-01期间内的最大值合成的NDVI。

代码：

var table=ee.FeatureCollection("users/boyxiaozheng/feature/beijing");
var landsat8=ee.ImageCollection("LANDSAT/LC08/C01/T1_RT");

function cal_ndvi(image){
  var mask=ee.Algorithms.Landsat.simpleCloudScore(image).select(['cloud']).lte(20);
  var ndvi=image.normalizedDifference(['B5', 'B4']).updateMask(mask);
  return ndvi;
}
var landsat_images=landsat8.filterDate('2017-06-01', '2017-10-01')
                .filterBounds(table.geometry())
                .filter(ee.Filter.lt('CLOUD_COVER',10))
                .map(cal_ndvi);


var visParam = {
 min: 0,
 max: 1,
 palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
   '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};

Map.setCenter(115,39);
Map.addLayer(table, {}, 'default display');
Map.addLayer(landsat_images.max().clip(table.geometry()),visParam);

运行结果（GEE显示）：