Google Earth Engine——绘制土地覆被变化制图

Google Earth Engine是一个用于行星尺度遥感分析的云计算平台。它提供对多 PB 卫星图像和地理空间数据集目录的访问，包括几乎所有 Landsat、MODIS 和 Sentinel 图像目录，以及 NDVI、土地覆盖分类和森林变化图等衍生产品。此外，Earth Engine 允许用户利用 Google 强大的计算基础设施来处理这些图像。

Earth Engine 最早和最引人注目的应用之一是马里兰大学的马特·汉森 (Matt Hansen) 和他的团队制作了全球森林变化图。在本教程中，我们将通过进行类似的分析来演示 Earth Engine 的使用，但规模要小得多。我们将在巴西的一个 Landsat 场景中绘制 2001 年到 2011 年之间的森林覆盖变化图。

加载影像

在本教程中，我们将制作一张巴西小区域 2001 年至 2011 年间森林变化的地图。为此，我们需要研究区域 2001 年和 2011 年的卫星图像，最好是在一年中的同一时间，以消除季节性作为混杂变量。Landsat 卫星在这里是一个很好的候选者：它们每 16 天重新访问同一区域，分辨率为 30m，并且可以在 Earth Engine 上找到可追溯到 1972 年的整个目录。

尝试在Earth Explorer 数据目录中搜索“landsat”以查看可用的数据集。您会注意到，图像集合可用于每个不同的 Landsat 卫星和不同级别的处理。我们将为Landsat 5和Landsat 7使用正射校正大气顶 (TOA) 反射收集。每个集合由唯一的ImageCollection ID（例如LANDSAT/LT5_L1T_TOA_FMASK）标识，并由该卫星拍摄的每个 Landsat 场景组成，每个场景都由唯一的Image ID（例如LANDSAT/LT5_L1T_TOA_FMASK/LT52240632011210CUB01）标识。

在本教程中，我们分别选择了 2001 年 7 月和 2011 年 7 月的两个 Landsat 场景，作为我们分析的基础。要加载和查看这些场景，请将以下代码粘贴到您的代码编辑器中：

 load landsat scenes for 2001 & 2011 in Brazil
// landsat 5 - 2011-07-29
var ls5 = ee.Image('LANDSAT/LT5_L1T_TOA_FMASK/LT52240632011210CUB01');
print(ls5);

// landsat 7 - 2001-07-09
var ls7 = ee.Image('LANDSAT/LE7_L1T_TOA_FMASK/LE72240632001190EDC00');
print(ls7);

// plot true colour composites
Map.addLayer(ls5, {bands: ['B3', 'B2', 'B1'], gamma: 1.5}, '2011');
Map.addLayer(ls7, {bands: ['B3', 'B2', 'B1'], gamma: 1.5}, '2001');
Map.centerObject(combined, 9);

 

在地图窗格中，您现在应该会看到这些 Landsat 场景被叠加。在地图的左上角，单击图层按钮以切换您看到的是 2001 年还是 2011 年的场景。{bands: ['B3', 'B2', 'B1']}调用中的使用Map.addLayers()告诉地球引擎要绘制哪些波段，在这种情况下会生成真彩色合成。尝试使用{bands: ['B4', 'B3', 'B2']}以获得标准的假色复合高亮植被。除了地图之外，print()代码中的语句还将图像的元数据打印到控制台（右上窗格）。例如，您可以在此处浏览图像以查看可用的频段。

遮罩云彩和阴影

这些图像被特别选择为具有非常有限的云量，但是，有一些我们希望掩盖这些多云像素并将它们从我们的分析中删除。幸运的是，我们正在使用的 Landsat 集合包含一个fmask带整数的额外波段，指示给定像素是阴影 (2)、雪 (3) 还是云 (4)。我们将提取此波段，将其转换为二进制（1 表示要保留的像素，0 表示要因云或阴影而遮盖的像素），然后使用它来遮盖原始图像。我还将重命名波段，以便我们区分年份。

 

 landsat scenes for 2001 & 2011 in Brazil
// landsat 5 - 2011
var ls5 = ee.Image('LANDSAT/LT5_L1T_TOA_FMASK/LT52240632011210CUB01');
// create binary mask band for clouds and shadows
var ls5Mask = ls5.select('fmask').lt(2);
// select bands 1-5 and rename to append year
ls5 = ls5.
  select(['B1', 'B2', 'B3', 'B4', 'B5']).
  rename(['B1_2011', 'B2_2011', 'B3_2011', 'B4_2011', 'B5_2011']);
// mask out clouds
ls5 = ls5.updateMask(ls5Mask);

// landsat 7 - 2001
var ls7 = ee.Image('LANDSAT/LE7_L1T_TOA_FMASK/LE72240632001190EDC00');
// create binary mask band for clouds and shadows
var ls7Mask = ls7.select('fmask').lt(2);
// select bands 1-5 and rename to append year
ls7 = ls7.
  select(['B1', 'B2', 'B3', 'B4', 'B5']).
  rename(['B1_2001', 'B2_2001', 'B3_2001', 'B4_2001', 'B5_2001']);
// mask out clouds
ls7 = ls7.updateMask(ls7Mask);

// plot true colour composites
Map.addLayer(ls5, {bands: ['B3_2011', 'B2_2011', 'B1_2011'], gamma: 1.5}, '2011');
Map.addLayer(ls7, {bands: ['B3_2001', 'B2_2001', 'B1_2001'], gamma: 1.5}, '2001');
Map.centerObject(ls7, 9);

请注意，这些 Landsat 图像现在在云层所在的位置出现了洞。

创建多日期堆栈

我们有兴趣识别在 2001 年到 2011 年间从森林变为非森林（反之亦然）的像素。因此我们需要将这两个场景叠加起来，以便我们可以同时分析它们。我通过将 2011 年的红外波段映射到红色并将 2001 年的红外波段映射到绿色和蓝色来可视化这个堆栈。在这个复合图中，砍伐森林的区域将显示为红色。

 landsat scenes for 2001 & 2011 in Brazil
// landsat 5 - 2011
var ls5 = ee.Image('LANDSAT/LT5_L1T_TOA_FMASK/LT52240632011210CUB01');
// create binary mask band for clouds and shadows
var ls5Mask = ls5.select('fmask').lt(2);
// select bands 1-5 and rename to append year
ls5 = ls5.
  select(['B1', 'B2', 'B3', 'B4', 'B5']).
  rename(['B1_2011', 'B2_2011', 'B3_2011', 'B4_2011', 'B5_2011']);
// mask out clouds
ls5 = ls5.updateMask(ls5Mask);

// landsat 7 - 2001
var ls7 = ee.Image('LANDSAT/LE7_L1T_TOA_FMASK/LE72240632001190EDC00');
// create binary mask band for clouds and shadows
var ls7Mask = ls7.select('fmask').lt(2);
// select bands 1-5 and rename to append year
ls7 = ls7.
  select(['B1', 'B2', 'B3', 'B4', 'B5']).
  rename(['B1_2001', 'B2_2001', 'B3_2001', 'B4_2001', 'B5_2001']);
// mask out clouds
ls7 = ls7.updateMask(ls7Mask);

// stack to create multi-date composite
var combined = ls7.addBands(ls5);
// mask clouds
var cloudMask = combined.updateMask(ls7Mask).updateMask(ls5Mask);

// plot true colour composites
Map.addLayer(ls5, {bands: ['B3_2011', 'B2_2011', 'B1_2011'], gamma: 1.5}, '2011');
Map.addLayer(ls7, {bands: ['B3_2001', 'B2_2001', 'B1_2001'], gamma: 1.5}, '2001');
// plot multi-date composite
Map.addLayer(combined, {bands: ['B5_2011', 'B5_2001', 'B5_2001'], gamma: 1.5}, 'combined');
Map.centerObject(combined, 9);

尝试使用多日期图像识别森林砍伐区域，然后在图层管理器中关闭此图像并在 2001 年和 2011 年图像之间滑动以确认确实发生了森林砍伐。

训练多边形

我们将使用随机森林算法执行监督森林变化分类。为此，我们需要为该算法提供训练数据，即我们需要描绘已知命运的区域以表征不同类别的光谱特征。特别是，我们对以下 4 个类感兴趣：

class	name	description
0	forest	forest in 2000 and 2011
1	forestLoss	forest loss
2	nonforest	non-forest in 2000 and 2011
3	forestGain	forest gain

Earth Engine 提供了一个直接在地图上绘制多边形的工具，我们将使用它来为四个类中的每一个描绘多边形。目标是尽可能多地捕获类内的可变性，因此我们将尝试为每个类使用 4-6 个多边形来捕获类的全部多样性。首先，单击地图左上角的多边形工具，这将创建一个新的多边形图层。

现在点击这个新图层旁边的齿轮图标，并在下图中将细节填充为荧光笔。该层将用于第一类：2001 年到 2011 年之间未更改的森林。您需要为其命名（“forest”），将类型设置为“特征”而不是几何，并添加一个新的“类”属性（这将是 0 级）。

现在使用所有三层卫星图像来识别在整个研究期间仍然是森林的区域，并创建一些多边形来描绘这些区域。保持你的多边形小，记住捕捉这个类中的多样性。一旦你完成了这个类，移动到其他类，确保为每个类创建一个新层，并在层属性中填写正确的信息。调用您的图层：forest、  forestLoss、nonforest和forestGain。

完成后，向上滚动到代码编辑器的顶部，您将看到导入这些多边形图层的新部分。此时请务必保存您的代码，以免丢失这些多边形！

提取单元格值

最后，我们需要将这四个训练层合二为一，并从多边形内提取图像单元格值。这将生成一个表格，将每个类别的像素与这些像素中的光谱带值相关联。很可能当您选择多边形时，某些类（例如未变化的森林）很容易找到示例，因此这些类的训练多边形覆盖了更大的区域。理想情况下，我们希望每个班级有相同数量的训练单元。此外，Earth Engine 施加了使用限制，如果训练多边形包含太多单元格，则会超出这些限制并返回错误。为了解决这个问题，我们将在多边形内进行二次采样。将以下代码添加到您已有的脚本的末尾。

// subsample training polygons with random points
// this ensures all classes have same sample size
// also EE can't handle too many cells at once
var trainingLayers = [forest, forestLoss, nonforest, forestGain];
var n = 500;
// loop over training layers
for (var i = 0; i < trainingLayers.length; i++) { 
  // sample points within training polygons
  var pts = ee.FeatureCollection
    .randomPoints(trainingLayers[i].geometry(), n);
  // add class
  var thisClass = trainingLayers[i].get('class');
  pts = pts.map(function(f) {
    return f.set({class: thisClass});
  });
  // extract raster cell values
  var training = combined.sampleRegions(pts, ['class'], 30);
  // combine trainging regions together
  if (i === 0) {
    var trainingData = training;
  } else {
    trainingData = trainingData.merge(training);
  }
}

随机森林

我们现在有了图像和训练数据，是时候运行随机森林分类了。将以下代码添加到您的脚本中以拟合随机森林模型并绘制生成的森林变化图。

 classify with random forests
// use bands 1-5 from each time period
var bands = ['B1_2001', 'B1_2011', 'B2_2001', 'B2_2011', 'B3_2001', 'B3_2011',
             'B4_2001', 'B4_2011', 'B5_2001', 'B5_2011'];
// fit a random forests model
var classifier = ee.Classifier.randomForest(30)
  .train(trainingData, 'class', bands);
// produce the forest change map
var classified = combined.classify(classifier);
var p = ['00ff00', 'ff0000', '000000', '0000ff'];
// display
Map.addLayer(classified, {palette: p, min: 0, max: 3}, 'classification');

下面是一个示例森林变化图，您的可能会略有不同，因为您可能选择了不同的训练区域。在这张地图中，forest为绿色，nonforest为黑色，forestLoss为红色，forestGain为蓝色。

准确度评估

在我们使用刚刚创建的地图之前，了解它的准确度很重要。例如，如果分类地图显示给定区域发生了森林损失，那么我们对该区域实际经历森林损失的信心有多大？

混淆矩阵是评估分类算法性能的标准方法。它采用已知类别的案例（例如训练数据或独立验证数据集）并将它们与预测类别进行比较。矩阵的行是实际类的实例，而列是预测类的实例。矩阵的对角线给出了正确分类的数量，而非对角线给出了错误分类的数量。例如，如果我们只有两个类，矩阵可能如下所示：

在这个例子中，12 类的 10 例被正确分类，而 2 类的 8 例中有 5 例被正确分类 查看非对角线分量，在 2 例中，1 类被错误地分配到 2 类，3案例类别 2 被错误地分配给类别 1。总体准确度是正确分类的总数占案例总数的比例，在这种情况下为15/20=75.

要计算森林变化图的混淆矩阵和整体准确度，请将以下代码添加到脚本末尾：

// accuracy assessement
var confMat = classifier.confusionMatrix();
print('Confusion matrix: ', confMat);
print('Overall accuracy: ', confMat.accuracy());

您现在应该看到控制台的混淆矩阵和整体准确性。

在这种情况下，我们在所有类中的准确性都非常好，但是请注意，我们使用训练数据来执行验证。在实践中，最好收集一个独立的验证数据集，或者将训练数据集划分为训练和验证子集，以避免准确性评估中的偏差。

结论

在本教程中，我们使用监督分类为巴西的单个 Landsat 场景构建了森林变化图。然而，这只是触及了 Earth Engine 所能实现的功能的表面。我们可以扩展我们的分析以包括更大的区域或研究不同地理位置或生物群落的土地变化。而且，应用不仅限于保护背景下的土地覆盖变化，Earth Engine 广泛适用于任何需要分析地球表面时空趋势的任务。