[人工智能]【Earth Engine】合成Landsat4/5/7/8/9影像并进行NDVI、NDWI和MNDWI等指数计算

1 简介与摘要

最近在做一个课题，需要逐月还有逐年的NDVI还有一系列衍生数据。
但是逐月的使用Sentinel影像会出现有云，并且掩膜掉容易出现空洞，所以我就想着可以拿Landsat来补充。
同样的思想，在逐年数据上也可以应用于不同的Landsat影像，这样可以丰富影像数量，提升计算质量。比如前一段时间刚上去的Landsat9，就号称可以和Landsat8协同，降低重访周期。
考虑到使用同一系列卫星的场景比较多，所以这篇博文就使用Landsat系列卫星的影像合成。

因此，结合实际应用场景，在这里需要做的工作是：

1. 输入某一个时间（或时间段），得到ImageCollection，并且（以地图和数字的形式）显示兴趣区里的影像数量；
2. 通过ImageCollection数据合成，并且计算一些基础的指数。

2 思路

首先，要获取一个包含我需要的时间段/影像质量/区域范围的Landsat系列影像，也就是一个Landsat的ImageCollection。
在得到ImageCollection之前，需要考虑一个问题——这些影像的相同的波段名可能对应不同的频率，比如说Landsat4/5/7的第四波段是近红外（NIR），而Landsat8/9的第四波段是红光（red），所以需要在这些波段重命名一下，这样合成才不会合错。
然后，就可以通过一个mean()（当然也可以其他的函数），将一个ImageCollection合成一个Image，以方便计算各种指数。
最后，通过预设的指数计算公式，计算指数即可。

总的来说，比较麻烦的是前半部分。后半部分的工作直接套用一些常规的做法就行了。

3 效果预览

惯例，在代码前面先放效果预览。
我在控制台展示了各个卫星匹配到的影像数量，还有合成后总的ImageCollection的影像数量。
在地图中显示合成后的真彩色和假彩色影像，还有计算后的NDVI、NDWI、MNDWI，最后还有各个区域影像数量的可视化。

在接下来的示例里，我时间段选取20200101——20210101，云量低于20%，兴趣区随手勾了一个大概希腊+土耳其（老精罗了）。

我的兴趣区（roi）：

控制台：

上图的意思是Landsat8得到1484幅影像，Landsat7得到1411幅影像，合成后的ImageCollection有2895幅（1484+1411）影像。

合成后的真彩色影像：

合成后的假彩色影像：

NDVI（绿色是高值，红色是低值）：

NDWI（蓝色是高值，白色是低值）：

MNDWI（蓝色是高值，白色是低值）：

影像数量（以像元为基本计算单位，红色是低值，绿色是中值，蓝色是高值）：

4 代码思路

考虑到要调用Landsat4/5/7/8/9，所以先写上这些卫星的预处理函数，包括大气校正和云掩膜；其中4/5/7的函数通用，8/9的函数通用。如下：

function applyScaleFactorsL89(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}

function applyScaleFactorsL457(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBand, null, true);
}

function cloudmaskL89(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 4);
  var cloudsBitMask = (1 << 3);
  // Get the pixel QA band.
  var qa = image.select('QA_PIXEL');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

function cloudMaskL457(image) {
  var qa = image.select('QA_PIXEL');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 3)
                  .and(qa.bitwiseAnd(1 << 9))
                  .or(qa.bitwiseAnd(1 << 4));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
}

写完预处理的函数，就需要改波段名称，以便放在同一个ImageCollection里能自动合成。因此，需要写将波段重命名的函数，如下：

function bandRenameL89(image) {
  var blue = image.select(['SR_B2']).rename('blue');
  var green = image.select(['SR_B3']).rename('green');
  var red = image.select(['SR_B4']).rename('red');
  var nir = image.select(['SR_B5']).rename('nir');
  var swir1 = image.select(['SR_B6']).rename('swir1');
  var swir2 = image.select(['SR_B7']).rename('swir2');
  var new_image = blue.addBands([green, red, nir, swir1, swir2]);
  return new_image;
}

function bandRenameL457(image) {
  var blue = image.select(['SR_B1']).rename('blue');
  var green = image.select(['SR_B2']).rename('green');
  var red = image.select(['SR_B3']).rename('red');
  var nir = image.select(['SR_B4']).rename('nir');
  var swir1 = image.select(['SR_B5']).rename('swir1');
  var swir2 = image.select(['SR_B7']).rename('swir2');
  var new_image = blue.addBands([green, red, nir, swir1, swir2]);
  return new_image;
}

在这里，选取常用的光学波段，将波段重命名为blue、green、red、nir、swir1、swir2。
最后，根据波段名称，写指数计算的函数。和其他人的文章不同，这里的输入的波段名称是我们重命名后的（见上面的代码）。这里以NDVI、NDWI和MNDWI为例，如下：

function NDVI(img) {
 var nir = img.select("nir");
 var red = img.select("red");
 var ndvi = img.expression(
   "(nir - red)/(nir + red)",
   {
     "nir": nir,
     "red": red
   }
 );
 return ndvi;
}

function NDWI(img) {
 var green = img.select("green");
 var nir = img.select("nir");
 var ndwi = img.expression(
   "(green - nir)/(green + nir)",
   {
     "green": green,
     "nir": nir
   }
 );
 return ndwi;
}

function MNDWI(img) {
 var green = img.select("green");
 var swir1 = img.select("swir1");
 var mndwi = img.expression(
   "(green - swir1)/(green + swir1)",
   {
     "green": green,
     "swir1": swir1
   }
 );
 return mndwi;
}

好了，到这里，所有预置的函数都写完了。
接下来就是对各个卫星的影像预处理，然后放进同一个ImageCollection里。下面是利用上述的函数进行预处理，可以根据需要自行调整。如下：

// landsat 9
var l9_col = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL89)
                  .map(cloudmaskL89)
                  .map(bandRenameL89)
                  ;
print('landsat9', l9_col.size())
// landsat 8
var l8_col = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL89)
                  .map(cloudmaskL89)
                  .map(bandRenameL89)
                  ;
print('landsat8', l8_col.size())
// landsat 7
var l7_col = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL457)
                  .map(cloudMaskL457)
                  .map(bandRenameL457)
                  ;
print('landsat7', l7_col.size())
// landsat 5
var l5_col = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL457)
                  .map(cloudMaskL457)
                  .map(bandRenameL457)
                  ;
print('landsat5', l5_col.size())
// landsat 4
var l4_col = ee.ImageCollection('LANDSAT/LT04/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL457)
                  .map(cloudMaskL457)
                  .map(bandRenameL457)
                  ;
print('landsat4', l4_col.size())

需要注意的是，roi是我们划定的研究区，sart_date和end_date是我们需要的影像时间范围，cloudCover是我们预设的最大云量，为了方便修改参数，一些影像要求（日期和云量）在这里我使用变量代替。另外，在每个卫星预处理函数后面，用size()函数print出了匹配到的影像数量。效果图见前一章。
接下来的工作就是把这些影像，合成一个ImageCollection。这里需要用到的函数是merge()，具体如下：

var image = l9_col
            .merge(l8_col)
            .merge(l7_col)
            .merge(l5_col)
            .merge(l4_col)
            ;
print("final image count", image.size(), image)

最后，调用指数计算的函数对最终合成的ImageCollection进行计算，不过在此之前需要先mean合成一下，如下：

print("final image count", image.size(), image)
var final_image = image.mean().clip(roi);
var image_ndvi = NDVI(final_image)
var image_ndwi= NDWI(final_image)
var image_mndwi = MNDWI(final_image)

最后简单写个可视化，包括NDVI、NDWI、MNDWI：

Map.addLayer(final_image, {bands: ["red", "green", "blue"], min:0.0, max:0.25}, "image")
Map.addLayer(final_image, {bands: ["nir", "red", "green"], min:0.0, max:0.25}, "image2")
var ndvi_palettes = ["#e700d5", "#e60000", "#e69f00", "#dfe200", "#7ebe00", "#00a10c", "#008110"];
Map.addLayer(image_ndvi.clip(roi), {min:-0.5, max:0.7, palette:ndvi_palettes}, "ndvi");
var ndwi_palettes = ["ffffff","#f9f9f9","#d8fdf4","#7dd5e9","3d7ede","243ad4","#1c00b8", "#250081"];
Map.addLayer(image_ndwi.clip(roi), {min:-0.5, max:0.7, palette:ndwi_palettes}, "ndwi");
Map.addLayer(image_mndwi.clip(roi), {min:-0.5, max:0.9, palette:ndwi_palettes}, "mndwi");

然后，顺手把影像数量的可视化写了，这里用到count()，如下：

var images_count = image.select('red').count();
Map.addLayer(images_count.rename('count').clip(roi), {min:10, max:80, palette:["#ff0000", "#fbff00", "#1dff02", "#02adff"]}, "images count");

好了，完成了。最后导出就行。

Export.image.toDrive(
  {
    image: image_ndvi.clip(roi),
    folder: "ndvi",
    description: "ndvi" + year_name,
    scale: 30,
    region: roi,
    maxPixels: 1e13
  })
Export.image.toDrive(
  {
    image: image_ndvi.clip(roi),
    folder: "ndwi",
    description: "ndwi" + year_name,
    scale: 30,
    region: roi,
    maxPixels: 1e13
  })
Export.image.toDrive(
  {
    image: image_ndvi.clip(roi),
    folder: "mndwi",
    description: "mndwi" + year_name,
    scale: 30,
    region: roi,
    maxPixels: 1e13
  })

对了，要记得在最前面补上需要更改的参数。

var year_name = 2020;
var start_date = (year_name) + '-01-01';
var end_date   = (year_name + 1) + '-01-01';
var cloudCover = 20

5 完整代码

我的代码链接在这，可以直接使用。
完整代码如下（和链接中相同）：

var year_name = 2020;
var start_date = (year_name) + '-01-01';
var end_date   = (year_name + 1) + '-01-01';
var cloudCover = 20

Map.addLayer(roi)

// indices
function NDVI(img) {
 var nir = img.select("nir");
 var red = img.select("red");
 var ndvi = img.expression(
   "(nir - red)/(nir + red)",
   {
     "nir": nir,
     "red": red
   }
 );
 return ndvi;
}

function NDWI(img) {
 var green = img.select("green");
 var nir = img.select("nir");
 var ndwi = img.expression(
   "(green - nir)/(green + nir)",
   {
     "green": green,
     "nir": nir
   }
 );
 return ndwi;
}

function MNDWI(img) {
 var green = img.select("green");
 var swir1 = img.select("swir1");
 var mndwi = img.expression(
   "(green - swir1)/(green + swir1)",
   {
     "green": green,
     "swir1": swir1
   }
 );
 return mndwi;
}


function bandRenameL89(image) {
  var blue = image.select(['SR_B2']).rename('blue');
  var green = image.select(['SR_B3']).rename('green');
  var red = image.select(['SR_B4']).rename('red');
  var nir = image.select(['SR_B5']).rename('nir');
  var swir1 = image.select(['SR_B6']).rename('swir1');
  var swir2 = image.select(['SR_B7']).rename('swir2');
  var new_image = blue.addBands([green, red, nir, swir1, swir2]);
  return new_image;
}

function bandRenameL457(image) {
  var blue = image.select(['SR_B1']).rename('blue');
  var green = image.select(['SR_B2']).rename('green');
  var red = image.select(['SR_B3']).rename('red');
  var nir = image.select(['SR_B4']).rename('nir');
  var swir1 = image.select(['SR_B5']).rename('swir1');
  var swir2 = image.select(['SR_B7']).rename('swir2');
  var new_image = blue.addBands([green, red, nir, swir1, swir2]);
  return new_image;
}

function applyScaleFactorsL89(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}

function applyScaleFactorsL457(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBand, null, true);
}

function cloudmaskL89(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 4);
  var cloudsBitMask = (1 << 3);
  // Get the pixel QA band.
  var qa = image.select('QA_PIXEL');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

function cloudMaskL457(image) {
  var qa = image.select('QA_PIXEL');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 3)
                  .and(qa.bitwiseAnd(1 << 9))
                  .or(qa.bitwiseAnd(1 << 4));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
}



// get image collection
// landsat 9
var l9_col = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL89)
                  .map(cloudmaskL89)
                  .map(bandRenameL89)
                  ;
print('landsat9', l9_col.size())
// landsat 8
var l8_col = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL89)
                  .map(cloudmaskL89)
                  .map(bandRenameL89)
                  ;
print('landsat8', l8_col.size())
// landsat 7
var l7_col = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL457)
                  .map(cloudMaskL457)
                  .map(bandRenameL457)
                  ;
print('landsat7', l7_col.size())
// landsat 5
var l5_col = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL457)
                  .map(cloudMaskL457)
                  .map(bandRenameL457)
                  ;
print('landsat5', l5_col.size())
// landsat 4
var l4_col = ee.ImageCollection('LANDSAT/LT04/C02/T1_L2')
                  .filterBounds(roi)
                  .filterDate(start_date, end_date)
                  .filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
                  .map(applyScaleFactorsL457)
                  .map(cloudMaskL457)
                  .map(bandRenameL457)
                  ;
print('landsat4', l4_col.size())

// combine, mean and calculate
var image = l9_col
            .merge(l8_col)
            .merge(l7_col)
            .merge(l5_col)
            .merge(l4_col)
            ;
print("final image count", image.size(), image)
var final_image = image.mean().clip(roi);
var image_ndvi = NDVI(final_image)
var image_ndwi= NDWI(final_image)
var image_mndwi = MNDWI(final_image)

Map.addLayer(final_image, {bands: ["red", "green", "blue"], min:0.0, max:0.25}, "image")
Map.addLayer(final_image, {bands: ["nir", "red", "green"], min:0.0, max:0.25}, "image2")
var ndvi_palettes = ["#e700d5", "#e60000", "#e69f00", "#dfe200", "#7ebe00", "#00a10c", "#008110"];
Map.addLayer(image_ndvi.clip(roi), {min:-0.5, max:0.7, palette:ndvi_palettes}, "ndvi");
var ndwi_palettes = ["ffffff","#f9f9f9","#d8fdf4","#7dd5e9","3d7ede","243ad4","#1c00b8", "#250081"];
Map.addLayer(image_ndwi.clip(roi), {min:-0.5, max:0.7, palette:ndwi_palettes}, "ndwi");
Map.addLayer(image_mndwi.clip(roi), {min:-0.5, max:0.9, palette:ndwi_palettes}, "mndwi");


var images_count = image.select('red').count();
Map.addLayer(images_count.rename('count').clip(roi), {min:10, max:80, palette:["#ff0000", "#fbff00", "#1dff02", "#02adff"]}, "images count");



// export to drive
Export.image.toDrive(
  {
    image: image_ndvi.clip(roi),
    folder: "ndvi",
    description: "ndvi" + year_name,
    scale: 30,
    region: roi,
    maxPixels: 1e13
  })
Export.image.toDrive(
  {
    image: image_ndvi.clip(roi),
    folder: "ndwi",
    description: "ndwi" + year_name,
    scale: 30,
    region: roi,
    maxPixels: 1e13
  })
Export.image.toDrive(
  {
    image: image_ndvi.clip(roi),
    folder: "mndwi",
    description: "mndwi" + year_name,
    scale: 30,
    region: roi,
    maxPixels: 1e13
  })

6 后记

刚开始学习GEE，可能有些地方不太专业不太科学，希望诸位同行前辈不吝赐教或者有什么奇妙的想法可以和我共同探讨一下。可以在csdn私信我或者联系我邮箱（chinshuuichi@qq.com），不过还是希望大家邮箱联系我，csdn私信很糟糕而且我上csdn也很随缘。

如果觉得有帮助，希望大家能打赏支持一下（我的乞讨码在下面）~
-----------------------分割线（以下是乞讨内容）-----------------------