3在 Google Earth Engine 中使用 Landsat 8 进行预处理

在 Google Earth Engine 中提供数据之前，拥有（免费！）专门的支持和幕后工作是一个令人难以置信的优势。但是，您可能仍然发现有必要操作您感兴趣的数据集以促进特定的研究应用。在本模块中，我们将使用 Landsat 8 数据，下图详细说明了不同处理级别的几个用例。
 

来自Young 等人的决策工作流程，2017 年，显示了针对不同级别的 Landsat 数据预处理的建议用例。

3.1预处理级别示例。

为了定性地了解不同预处理级别之间的差异，我们可以查看 2018 年夏末美国俄勒冈州南部的几张真彩色图像。在这个时间范围内，来自加利福尼亚州北部的卡尔火灾。为了评估我们的初始图像，我们将加载“原始”Landsat 8 集合。原始数据（也称为“传感器处辐射”）尚未针对任何形式的潜在影响进行校正，并且通常不用于生态研究。然而，为即将到来的水平建立一个基线是有帮助的。运行下面的脚本以生成类似于下面的图像。

注：除以下预处理级别外，Landsat 数据已分为两个质量级别，Tier 1 和 Tier 2。Tier 1 是更高质量的选项。单击此处了解有关各层之间差异的更多信息。

// 设置影像中心和缩放倍数
Map.setCenter(-122.3158, 42.4494, 12);

//导入影像
var raw = ee.Image('LANDSAT/LC08/C01/T1/LC08_045031_20180811');

// 可视化参数
var rawvis = {bands: ['B4','B3','B2'], min: 0.0, max: 30000.0, gamma: 1};

// 加载图层
Map.addLayer(raw, rawvis, 'raw');

来自俄勒冈州南部的原始 Landsat 8 图像。我们的图像以 Mt McLoughlin (2,893 m) 为中心，东边是上克拉马斯湖。如果缩小图像，可以看到 Medford（西）和 Klamath Falls（东）这两个城市。

3.2大气层顶（TOA）

下一级预处理采用我们的“原始”数据，并对太阳活动的影响应用修正，包括太阳辐照度、地球-太阳距离和太阳仰角。对于研究人员来说，大气顶部 (TOA) 通常适用于评估单一日期、单一场景的图像（即相对较小的研究区域中的土地覆盖分类）。这是由于不同程度的太阳效应取决于收集的日期、时间和纬度。将以下代码附加到您的脚本中，生成的图像将类似于下面的图像。

// 对比toa数据
var toa = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_045031_20180811');

// 定义真彩色影像
var toavis = {bands: ['B4','B3','B2'], min: 0.0, max: 0.4};

// 加载toa图层
Map.addLayer(toa, toavis, 'toa');

加载 TOA 集合，由于对太阳效应的影响进行了一些校正，部分图像看起来更清晰。然而，似乎还有工作要做！

3.3表面反射率（SR）

这些数据经过了最高级别的预处理，试图最好地代表地面上的实际情况，其中一定量的太阳能反射回（反射）到空中和星载传感器。然而，即使是表面反射率产品也会受到太阳角度低、云层过多以及覆盖位置超过南北纬 65 度的不利影响（Young 等，2017）。尽管如此，建议对多个日期（例如变化检测）或大地理范围（例如算法预测）的分析使用 Landsat 表面反射率数据。将这最后一段代码添加到您的脚本中，以在您的地图查看器窗格中查看下面的图像。

// 加载一张sr影像
var sr = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_045031_20180811');

// 定义真彩色参数
var srvis = {bands: ['B4','B3','B2'], min: 0, max: 3000, gamma: 1.4};

// 加载图层
Map.addLayer(sr, srvis, 'sr');

应用大气校正似乎大大提高了我们的图像清晰度，尤其是在上克拉马斯湖和麦格洛林山以西的农业区。

3.4 Landsat 影像比较的完整代码

这些示例的代码来自一个特定的 Landsat 图像，但您现在拥有调查任何感兴趣区域（即您的研究区域）以比较不同预处理级别的框架。

// 设置图形展示中心
Map.setCenter(-122.3158, 42.4494, 12);

// 导入一景影像
var raw = ee.Image('LANDSAT/LC08/C01/T1/LC08_045031_20180811');

//定义真彩色影像
var rawvis = {bands: ['B4','B3','B2'], min: 0.0, max: 30000.0, gamma: 1};

// 加载此图层
Map.addLayer(raw, rawvis, 'raw');

//比较，从第1级陆地卫星8号大气层顶部数据加载一张图像。
var toa = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_045031_20180811');

//定义真彩色影像
var toavis = {bands: ['B4','B3','B2'], min: 0.0, max: 0.4};

// 加载此图层
Map.addLayer(toa, toavis, 'toa');

// 加载地表反射率数据
var sr = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_045031_20180811');

// 定义真彩色影像
var srvis = {bands: ['B4','B3','B2'], min: 0, max: 3000, gamma: 1.4};

// 加载sr影像
Map.addLayer(sr, srvis, 'sr');

3.5云遮蔽

正如我们发现的那样，当我们到达表面反射产品时，大气、地形和几何校正的预处理工作已经完成。在 Google Earth Engine 中提供 Landsat 数据之前，处理 Landsat 数据未包含的一个非常重要的步骤是去除近地天气现象。这通常以云的形式出现。云在热带（即茂密的雨林）和水体上空特别普遍。在本节中，我们将在明尼苏达州东北部和安大略省西南部的 Quetico-Superior 地区看看后者，那里的地形在夏季难以捕捉，因为大面积的云层覆盖，部分是由数百个小到中等大小的湖泊。

夏季云彩聚集在奎蒂科省立公园的让湖上空。图片来源：泰伊舒尔克。

3.5.1单图像屏蔽：第 1 部分

让我们从加载我们知道是多云的图像开始。在这个地区有很多选择可供选择，但这张八月下旬的图片向我们展示了多种形式的云。需要考虑的另一件事是，云具有在下面的土地上投下阴影的额外效果，这进一步扩大了我们最终必须移除的地理范围。启动一个新脚本，运行下面的代码以生成如下图所示的图像。

// 加载一个初始图像来测试查看多云的图像。
var sr = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_027026_20180829');

// 为真彩色图像定义visparams。
var srvis = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.4};

// 把你的多云图像添加到地图上。
Map.addLayer(sr, srvis, 'single_scene');
Map.setCenter(-91.8859, 48.8936, 8.81);

Quetico-Superior 国家上空的多云图像。

3.5.2单图像屏蔽：第 2 部分

现在来处理我们的云问题。Landsat 提供了一个 pixel_qa 波段，简而言之，它根据先前量化的特征（例如云和雾的可能性）分配不同的值。您会发现用于构建云遮罩的更复杂的代码，但这是从图像中删除那些讨厌的白色斑点的简单、保守的方法。将以下代码附加到您现有的脚本并重新运行以查看类似于下面的图像。请记住取消选中图层控件中的“single_scene”！

// 定义简单的云层掩码，基于 "pixel_qa "波段的数值。
//本质上，322=土地，324=水的值。
var qa = sr.select('pixel_qa');
var mask = qa.eq(322).or(qa.eq(324));
var sr_cm = sr.updateMask(mask);

// 加载图层

Map.addLayer(sr_cm, srvis, 'single_scene_masked');

成功去除了云和云阴影，但生成的图像不会留下很多可用像素。

3.5.3跨多个日期的掩码

在带有云遮罩的单个 Landsat 场景中，我们失去了相当多的地理覆盖范围。但还有另一种方式！我们还可以在一个日期范围内应用掩码。为此，我们需要创建一个函数，我们将在模块 9 中详细介绍。现在，使用下面的函数（和其余代码）来继续你的脚本。将代码附加到您现有的脚本中。

虽然谷歌地球引擎中确实有来自 Landsat 的预制 NDVI 图像集，但这些数据集仅在 2017 年可用。因此，我们还将计算 NDVI 并将其添加到我们的图像集。这将使我们能够为 2018 年整个生长季节的每个像素生成一个中值，测量研究区域的植被健康状况。重新运行您的代码，生成的图像应如下所示。

// 定义去云函数，这里用QA波段
var cloudMask = function(image) {
  // var mask = image.select('pixel_qa').eq(322).or(image.select('pixel_qa').eq(324));
  var mask = image.select('pixel_qa').eq(322);
  return image.mask(mask);
};

// 使用NDVI
var collection = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate('2018-06-01','2018-09-30')
.filter(ee.Filter.calendarRange(140, 270)) // roughly, our potential growing season
.filterMetadata('WRS_PATH', 'equals', 27)
.filterMetadata('WRS_ROW', 'equals', 26)
.map(cloudMask)
.map(function NDVI(i){
  return i.addBands(i.normalizedDifference(['B5','B4']).rename('NDVI'))  // generates NDVI band
});

var median = collection.median();

// 红色=低值，深绿色=高值
var collectionVis = {bands: 'NDVI', min: 0.5, max: 0.95, palette: ['red','yellow','green','003300']};

Map.addLayer(median, collectionVis, 'NDVI');

可视化 2018 年整个生长季节的 NDVI 中值。随意切换不同的背景层。您的图像可能与此处显示的图像不完全相同。

3.5.4可视化图像计数

我们的 NDVI 图像看起来很不错。但我们对这些价值观有多自信？具体来说，重要的是评估有多少图像实际上构成了我们的中值。我们可以通过可视化每个像素位置使用的图像数量的总和，快速检查我们得到的 NDVI 中位数是否在整个研究区域具有代表性。将以下代码附加到您的脚本中，然后单击“运行”。您应该会看到类似于下图的图像。

对于计数层，如果我们在 NDVI 值中发现空间异常，则有多种选择。我们可以扩大我们的季节性日期范围或选择包括多年的数据。最终，我们可能会接受失败并决定天气多云而无法使用并决定探索不同的数据集——这是完全可以接受的！

// 统计ndvi波段数量
var counts = collection.select('NDVI').count();

//获取影像数量
print(collection.size());

// 在这个调色板中，红色=低值，蓝色=高值。
var countVis = {min: 0, max: 6, palette: ['b2182b','ef8a62','fddbc7', 
                                          'd1e5f0','67a9cf','2166ac']};

Map.addLayer(counts, countVis, 'counts');

可视化用于计算每个像素的中值的图像数量。较深的红色值较低，较深的蓝色值较高。

3.5.5生成和导出直方图

为了了解计数值的定量分布，我们可以构建计数数据的直方图。如果我们的研究区域多云以至于平均值仅代表一两个值，我们可能需要重新考虑我们的源数据或收集年份。将此最终代码块添加到您的脚本中，您将能够在“控制台”选项卡中看到直方图。您将制作自己的几何形状（如果您需要复习如何做到这一点，请参阅模块 1）。您的直方图可能会根据您的形状略有不同，但它足够大以包含图像，它应该类似于下图中的分布。

//画出你自己的矩形，并尝试覆盖大部分的Landsat场景 
// 然后，我们将制作一个分布的直方图。
//由于场景面积大，GEE要求我们进行一些聚合。 
// 让我们对我们的分布有一个大致的了解。
var histogram = ui.Chart.image.histogram(counts, geometry, 120);

// 展示直方图
print(histogram);

使用直方图函数量化来自计数层的值的分布数量。

// Load an initial image to test view a cloudy image.
var sr = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_027026_20180829');

// Define visparams for true-color image.
var srvis = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.4};

// Add your cloudy image to the map.
Map.addLayer(sr, srvis, 'single_scene')
Map.setCenter(-91.8859, 48.8936, 8.81);

// Define simple cloud mask, based in values from the 'pixel_qa' band.
// Essentially, values of 322 = land and 324 = water.
var qa = sr.select('pixel_qa');
var mask = qa.eq(322).or(qa.eq(324));
var sr_cm = sr.updateMask(mask);

// Add your new map for comparison to the cloudy image.
Map.addLayer(sr_cm, srvis, 'single_scene_masked');

// Define how your cloudMask function should work.
var cloudMask = function(image) {
  // var mask = image.select('pixel_qa').eq(322).or(image.select('pixel_qa').eq(324));
  var mask = image.select('pixel_qa').eq(322);
  return image.mask(mask);
};

// This time we'll look at a more ecologically-focused example, using NDVI.
var collection = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate('2018-06-01','2018-09-30')
.filter(ee.Filter.calendarRange(140, 270)) // roughly, our potential growing season
.filterMetadata('WRS_PATH', 'equals', 27)
.filterMetadata('WRS_ROW', 'equals', 26)
.map(cloudMask)
.map(function NDVI(i){
  return i.addBands(i.normalizedDifference(['B5','B4']).rename('NDVI'))  // generates NDVI band
});

var median = collection.median();

// Red = lower values and darker green = higher values.
var collectionVis = {bands: 'NDVI', min: 0.5, max: 0.95, palette: ['red','yellow','green','003300']};

Map.addLayer(median, collectionVis, 'NDVI');

// Create counts band.
var counts = collection.select('NDVI').count();

// Find our potential maximum count by getting the size of the initial image collection.
print(collection.size());

// In this palette, red = lower values and blue = higher values.
var countVis = {min: 0, max: 6, palette: ['b2182b','ef8a62','fddbc7', 
                                          'd1e5f0','67a9cf','2166ac']};

Map.addLayer(counts, countVis, 'counts');

// Draw your own rectangle and try to cover most of the Landsat scene and 
// then we'll make a histogram of the distribution. Because the scene area 
// is so large  GEE requires us to perform some aggregating. But this will 
// give us a general sense of our distribution.
var histogram = ui.Chart.image.histogram(counts, geometry, 120);

// Display the histogram.
print(histogram);

6参考文献

Young, NE, Anderson, RS, Chignell, SM, Vorster, AG, Lawrence, R., & Evangelista, PH (2017)。Landsat 预处理的生存指南。生态学，98（4），920-932。