1 简介与摘要

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

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

输入某一个时间（或时间段），得到ImageCollection，并且（以地图和数字的形式）显示兴趣区里的影像数量；
通过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 bandsvar 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 = 20Map.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 bandsvar 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也很随缘。
如果对你有帮助，还望支持一下~点击此处施舍或扫下图的码。
-----------------------分割线（以下是乞讨内容）-----------------------

【Earth Engine】合成Landsat4/5/7/8/9影像并进行NDVI、NDWI和MNDWI等指数计算相关推荐

Google Earth Engine（GEE）——利用插值方法解决影像去云后的空缺/填充/弥补方法详细讲解（拉萨区域为例）
本篇文章重点要解决的问题就是,很多时候我们无论是在小区域内的单景影像或者是中大尺度的影像,更或是长时间序列的影像研究中,很多情况下我们会因为云量筛选等因素,或多或少的存在影像空白而缺少值,因此如何处 ...
Google Earth Engine (GEE) ——export 导出指定尺寸的影像
本次错误主要是想导指定像素的大小的指定影像,但是再后续处理中缺出现了一些问题.具体问题如下: 原本是想导出一个完整的图像(带有 x 波段,对应于一个分类图,其他的输出到 S2波段) ,像素为64 x ...
Google Earth Engine（GEE）APP——一个监测影像各波段的DN值的app
这是一个非常简单的APP构造,主要用于监测地图上每一个点的波段信息,过程可以分为三步,第一步加载影像,第二步设置面板和小部件等一些信息,第三步,写一个回调函数用于出发和监测过程中的获取点位信息和波段信 ...
GEE_API Docs_Tutorials_1.编程基础和Earth Engine API入门
API Docs_Tutorials_1.编程基础和Earth Engine API入门一.Introduction to JavaScript for Earth Engine(JavaScrip ...
Google Earth Engine（GEE）——MODIS 影像LST地表温度随时间变化的趋势案例分析
该实验室的目标是使用 Google Earth Engine 深入研究气候变量.在本实验结束时,您将能够探索特定感兴趣区域的温度数据的长期趋势. MODIS LST 数据集 MOD11A2 V6 产品 ...
Google Earth Engine（GEE）农作物种植结构提取
目录写在前面 1.构建物候特征 2.构建光谱特征 3.将所有影像合并为一幅影像 4.构建随机森林算法进行分类 5.算法的存储 6.面积统计写在前面前段时间因为考研的原因一直没能更新,已经完成了农 ...
Google Earth Engine（GEE）填补缺失影像
今日分享: Google Earth Engine(GEE)填补缺失影像之前在做月合成NDVI的过程中,发现如果研究区较大时,一个月的影像覆盖不了整个研究区,就会有缺失的地方,还有就是去云之后,有云 ...
Google Earth Engine(GEE)——可视化动态图
代码: var geometry = /* color: #d63000 *//* shown: false *//* displayProperties: [{"type": & ...
Google Earth Engine （GEE） ——全球海岸线数全球海岸线数据集30米分辨率
全球海岸线数据集一个新的30米空间分辨率的全球海岸线矢量(GSV)是由2014年Landsat卫星图像的年度合成物开发的.图像的半自动分类是通过手动选择代表整个全球海岸线上的水和非水类别的训练点来完 ...

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

目录