Google Earth Engine（GEE）——10分钟短文快速了解地球引擎和森林面积损失计算

1. 了解什么是 Google 地球引擎

正如其开发人员所描述的那样，Google Earth Engine 是“世界上最先进的基于云的地理空间处理平台！” 这意味着，通过 Google Earth Engine，您可以访问并高效分析众多开源空间数据库（如 Landsat 和 MODIS 遥感图像、全球森林变化数据集、道路、保护区等）。在进行这些分析时，您使用的是 Google 服务器，因此您可以在您的计算机甚至一台花哨的计算机上进行需要数周甚至数月的分析。

从 Google 地球引擎，您可以将.csv您计算出的任何值的文件和geoTIFF文件（地理参考图像）导出到您的 Google 云端硬盘帐户。

2. 了解您可以使用 GEE 进行哪些类型的分析

使用 GEE，您可以以一种高效的方式回答大规模的研究问题，这在以前是不可能的，非常令人兴奋！您可以使用大型地理空间数据集来解决现代世界中人类面临的大量问题和挑战。稍后我们将了解如何探索 GEE 中可用的数据集，并且还可以导入您自己的地理参考图像（如无人机任务中的照片）。您可以从GEE 开发人员网站上的此页面了解如何导入您自己的栅格数据。

例如，您可以对不同的土地覆盖类型进行分类，您可以计算和提取景观特征的值，例如NDVI（归一化差异植被指数）——针对世界、特定感兴趣区域或世界上许多不同区域。真的，可能性是巨大的，在这里我们只是通过给你一个例子来说明如何使用 GEE 来计算森林覆盖率随时间的变化，从而触及表面。

如果您希望了解更多信息并练习您的 GEE 技能，可以查看Google 地球引擎开发者网站上的教程！

3.熟悉GEE布局

前往Earth Engine，开启你的GEE之旅！

花点时间熟悉一下 Earth Engine 编辑器的布局——就像刚开始学习一门新语言时，一次要吸收的东西似乎很多！使用您的空白脚本，尝试探索不同的选项卡。请注意，如果您绘制多边形或放置点，它们将如何出现在您的脚本中。您可以转到Inspector选项卡，单击地图中的某个地点，然后查看可用于该地点的信息。以下是大多数选项卡功能的概述：

4.学习JavaScript的基本原理

Google 地球引擎使用编程语言JavaScript。

与其他编程语言类似，也有在线支持——您可以谷歌JavaScript和 Earth Engine 教程。一开始看起来很陌生，但是多亏了在线编程社区，您很少会完全从头开始——也就是说，不要因为不能凭JavaScript头脑想出正确的代码而对自己感到难过马上。

JavaScript随着教程的进行，我们将向您介绍有关语法和函数的更多信息，但现在，这里有一些注意事项：

JavaScript以a 结尾的代码行;- 请注意，例如定义变量的代码可以分布在多行中，但您只需要将 a;放在代码块的最后一行的末尾。

要定义新变量，您可以使用：

在我们稍后创建的整个脚本中，您会在多个地方看到此代码的变体。本质上，当您导入数据集、创建新层、计算新值时，所有这些都需要存储为变量，以便您可以映射它们、导出它们等。

要在脚本中添加注释，请使用//. 例如，在您的空白新脚本的开头（如果您在探索时创建了任何多边形或点，您现在可以创建一个新脚本以开始“干净”）。与使用其他编程语言编码时一样，最好留下评论以确保您的脚本概述了您的身份、脚本的目的是什么以及您遵循特定工作流程的原因。以下是一些示例注释 - 您可以在脚本中编写类似的内容：

// Calculating forest cover change in protected areas around the world
// Gergana Daskalova
// 26th Nov 2018

在 JavaScript 中，您必须一次运行整个脚本 - 也就是说，例如，您不能选择脚本的两行并只运行它们，您必须运行整个脚本。Run您可以通过按下按钮来“运行”脚本。这意味着在整个教程中，当您向脚本中添加更多行时，您必须一直按Run才能查看添加的新代码的结果。

5.导入和探索数据——以保护区和森林覆盖变化为例

与任何分析一样，与其说是关于数据，不如说是关于你的研究问题，所以当你开始探索 GEE 时，请记住保留你的研究问题（或科学传播目标，因为 GEE 也非常适合这样做） ） 心里！

研究问题

世界各地不同国家公园的森林覆盖率如何变化？

在 GEE - 保护区导入和探索数据集

要导入保护区数据集（来自世界保护区数据库的世界各地保护区的多边形），protected area请在搜索选项卡中键入并选择数据库的多边形版本（另一个只是点，即一个坐标保护区内的点，而不是它们的轮廓）。

您导入的数据集出现在脚本的顶部 - 它目前被称为table不是特别有用的信息，因此您可以将其重命名为其他名称，例如parks.

请记住保存您的脚本并经常保存它！保存后，您会看到该文件出现在脚本选项卡下方的左侧。

可视化世界各地的保护区

接下来，我们将使用该Map函数来映射数据集，并添加一个图层。然后，您可以从地图窗口右上角的图层选项卡打开和关闭该图层。您还可以更改不透明度。

转到Inspector选项卡，单击地图上某处的一个点并查看features该点的名称 - 保护区的名称、面积、建立时间等。

环游世界，找到一个国家公园并“检查”它——你能找到名称、面积等吗——所有这些信息都在Inspector选项卡下。

在 GEE 中导入和探索数据集 - 森林覆盖变化 

与导入保护区数据集的方式类似，转到搜索选项卡，输入global forest change并选择Hansen 等人。数据集。

查看此数据集中包含的不同类型的信息 - 这将帮助您熟悉我们稍后的分析结果。

调用 objectgfc或任何您想要的其他名称，但请记住，如果您将其命名为其他名称，则必须gfc在接下来的所有代码中更改为您的新名称！接下来，我们将再次映射我们的数据集。

// Add the Global Forest Change dataset
Map.addLayer(gfc);

目前，我们只有一张黑色和红色的地图——黑色代表没有森林的地方，红色代表有森林覆盖的地方。这不是非常有用，在本教程的过程中，我们将努力使这张地图变得更好！

再次转到Inspector选项卡，单击红色部分地图上某处的一个点，然后查看features森林覆盖变化图层。如果它显示loss: 0, gain: 0，则意味着在这个特定像素中，没有发生森林损失或增加。

您还可以打开和关闭层，如果您不希望每次重新运行脚本时都执行该操作，则可以“注释掉”代码的某些部分。例如，映射保护区数据集需要相当长的时间，所以如果你不想多次这样做，你可以//在那行代码前面添加。//当您希望再次映射这些数据时，您可以随时删除。像这样：

// If you want to visualise the PAs around the world, you can use:
// Map.addLayer(parks);

如果您想将大量代码行转换为注释或将大量注释重新转换为代码，您可以在计算机Cmd + /上使用键盘Mac快捷键。Ctrl + /Windows

我们现在准备好改进我们的地图并得出森林损失和增加的量化值！

6.可视化森林覆盖变化

首先，最好定义分析的规模 - 在我们的例子中，它是 30 m，这是全球森林变化数据集的分辨率。如果给定像素经历了森林损失，这意味着在该 30 mx 30 m 正方形的某个地方，森林覆盖率有所下降。

您还可以将比例设置为自动检测数据集的分辨率并将其用作您的比例。

在脚本中键入以下代码：

// Set the scale for our calculations to the scale of the Hansen dataset
// which is 30m
var scale = gfc.projection().nominalScale();

下一步是为 2000 年（数据库启动时）的树木覆盖率、2016 年之前的损失和 2016 年森林覆盖率的增加创建变量。在栅格数据中，图像通常有不同的“波段”（例如，红色、绿色、UV），我们可以选择要使用的波段。在这种情况下，gfc对象的不同波段代表森林覆盖率、森林损失和森林增加，因此我们将为每个波段创建一个变量。

为此，我们将使用该select()功能。请注意，与其他编程语言（如）不同R，在JavaScript您将要应用函数的对象放在第一位，然后才是实际函数。

// Create a variable for the original tree cover in 2000
var treeCover = gfc.select(['treecover2000']);
 
// Convert the tree cover layer because the treeCover by default is in
// hundreds of hectares, but the loss and gain layers are just in hectares!
treeCover = treeCover.divide(100);
 
// Create a variable for forest loss
var loss = gfc.select(['loss']);
 
// Create a variable for forest gain
var gain = gfc.select(['gain']);

绘制全球森林覆盖、森林流失和森林增加地图

现在我们有了三个变量，我们可以为每个变量创建一个图层，然后我们可以使用我们选择的颜色绘制它们。我们将使用与Map.addLayer之前相同的功能，但除了添加对象名称外，我们还将指定颜色以及我们想要调用特定图层的内容。

请注意，我们还引入了一个新功能updateMask()。它所做的是掩盖 2000 年没有森林覆盖的区域——它们变得透明，所以我们不仅可以看到黑色，还可以看到海洋、河流、大陆轮廓等。

// Add the tree cover layer in light grey
Map.addLayer(treeCover.updateMask(treeCover),
    {palette: ['D0D0D0', '00FF00'], max: 100}, 'Forest Cover');
 
// Add the loss layer in pink
Map.addLayer(loss.updateMask(loss),
            {palette: ['#BF619D']}, 'Loss');
 
// Add the gain layer in yellow
Map.addLayer(gain.updateMask(gain),
            {palette: ['#CE9E5D']}, 'Gain');

请记住单击Run以便您看到新绘制的地图。如果您关闭绘制的前两层（保护区和通用 GFC 层），或者您可以保持保护区层打开，但通过拖动下面的栏来降低不透明度，森林层可能更容易看到层。

您可以使用十六进制代码指定颜色，这些代码#CE9E5D是上面代码中的数字和字母组合，例如黄色。您可以在线找到这些示例，例如本网站。

您还可以在地图视图和卫星视图之间切换。如果你足够放大并转到卫星视图，你实际上可以开始发现一些模式，比如亚马逊道路沿线的森林流失。

7. 计算特定位置的总森林覆盖增益和损失

至此我们可以看到森林发生损失和增加的地方，所以我们知道森林变化的程度，但是我们不知道森林变化的幅度，所以我们下一步是转换经历过的像素数增益或损失（请记住它们只是 0 或 1 值，0 表示否，1 表示是）到面积，例如平方公里。

对于我们之前创建的每个变量（森林覆盖、森林损失和森林增加），我们现在将创建代表森林覆盖、损失和增加面积的新变量。为了实现这一点，我们将使用该ee.Image.pixelArea()函数，并且我们必须使用multiply我们的原始变量（例如，treeCover），类似于当您将米转换为厘米时，您将乘以 100。这里我们希望面积以平方公里为单位，所以从平方米到平方公里，我们还要除以 1 000 000。最后，我们从新变量中选择第一个波段 - 分别是森林覆盖面积、损失面积和增益面积。

// The units of the variables are numbers of pixels
// Here we are converting the pixels into actual area
// Dividing by 1 000 000 so that the final result is in km2
var areaCover = treeCover.multiply(ee.Image.pixelArea())
                .divide(1000000).select([0],["areacover"]);
 
var areaLoss = loss.gt(0).multiply(ee.Image.pixelArea()).multiply(treeCover)
              .divide(1000000).select([0],["arealoss"]);
 
var areaGain = gain.gt(0).multiply(ee.Image.pixelArea()).multiply(treeCover)
              .divide(1000000).select([0],["areagain"]);

计算特定区域的森林损失和收益

我们通常对从世界各地特定地点的地理空间数据中提取价值感兴趣。在这里，我们的问题是关于国家公园森林覆盖率的变化，因此要回答这个问题，我们需要计算我们选择的国家公园而不是整个世界发生了多少森林覆盖率变化。

第一步是创建一个过滤变量，其中包含我们感兴趣的领域。在这里，我们将过滤parks包括世界上所有保护区的原始变量，减少到只有四个保护区。我们将使用ee.Filter.or()添加多个过滤条件。

// Create a variable that has the polygons for just a few
// national parks and nature reserves
var parks = parks.filter(ee.Filter.or(
    ee.Filter.eq("NAME", "Yellowstone"),
    ee.Filter.eq("NAME", "Sankuru"),
    ee.Filter.eq("NAME", "Cairngorms"),
    ee.Filter.eq("NAME", "Redwood")));

现在我们准备好计算森林损失和增加的面积，激动人心的时刻！我们将使用 GEE 术语中称为“reducer”的东西——一个汇总函数。我们将把它应用于我们的parks变量，我们将使用我们之前定义的比例（30m，数据集的分辨率）。结果将存储在两个新变量中，statsLoss和statsGain。

// Sum the values of loss pixels.
var statsLoss = areaLoss.reduceRegions({
  reducer: ee.Reducer.sum(),
  collection: parks,
  scale: scale
});
 
// Sum the values of gain pixels.
var statsGain = areaGain.reduceRegions({
  reducer: ee.Reducer.sum(),
  collection: parks,
  scale: scale
});

8.导出结果-汇总表

在这个阶段，我们已经计算了我们选择的保护区的森林损失和增加面积，但我们还没有真正看到或可视化这些数字。

我们可以导出.csv结果文件，在这种情况下，它们将转到您的 Google Drive 帐户。将下面的代码添加到您的脚本中，然后Run再次按下。你会看到Task选项卡亮起，去看看吧。你将有两个任务，你必须按下Run它们旁边的按钮（否则任务已经为你准备好了，但你还没有真正启动它们的完成），然后你会开始看到一个计时器 - 它反映了多少时间自您开始任务以来已通过。根据您的任务，它可能需要几秒到几小时。在我们的例子中应该是秒！

我们使用花括号来指定我们要导出的对象以及我们想要调用文件的内容，例如NP_forest_loss.

Export.table.toDrive({
  collection: statsLoss,
  description: 'NP_forest_loss'});
  
Export.table.toDrive({
  collection: statsGain,
  description: 'NP_forest_gain'});

去检查你在你的谷歌驱动器中的文件。一直向右滚动以查看该sum列，该列显示森林损失或增加的面积（以平方公里为单位）（取决于您正在查看的文件）。

9. R 中的进一步可视化——两全其美！

我们热衷于在我们的分析中纳入不同的平台和语言，发挥各自的优势。R并且R像这样的软件包ggplot2在您如何可视化您的发现方面提供了更大的灵活性，因此我们现在将切换R到制作我们研究的四个保护区的森林损失和收益的条形图。

注意：您还可以在 Earth Engine 中制作图表，因此这取决于个人喜好以及最适合您自己的工作流程的图表。您可以在开发者网站上找到有关如何在 Earth Engine 中创建图形的教程。

打开RStudio（或仅R取决于您的喜好）并通过转到开始一个新脚本File / New file / R Script。如果您以前从未使用R过，可以在此处找到我们的教程介绍。R

# Load libraries ----
library(ggplot2)
devtools::install_github('Mikata-Project/ggthemr') # to install the ggthemr package
# if you don't have it already
library(ggthemr)  # to set a custom theme but non essential!
library(forcats)  # to reorder categorical variables

我们可以使用ggthemr包为我们的情节设置一个主题（更改颜色和背景）。您可以在此处探索不同的颜色选项。

# Set theme for the plot
ggthemr('dust', type = "outer", layout = "minimal")
 
# This theme will now be applied to all plots you make, if you wanted to
# get rid of it, use:
# ggthemr_reset()

接下来，将您的工作目录设置为您保存我们导出到 Google 云端硬盘并读入文件的数据的位置。

# Read in the data ----
NP_forest_gain <- read.csv("NP_forest_gain.csv")
NP_forest_loss <- read.csv("NP_forest_loss.csv")

我们将合并这两个对象（一个用于森林损失，一个用于森林增加），以便我们可以在同一图中可视化它们。我们可以创建一个“标识符”列，以便我们知道哪些值指的是森林覆盖率的增加，哪些值指的是森林覆盖率的损失。

# Create identifier column for gain vs loss
NP_forest_gain$type <- "Gain"
NP_forest_loss$type <- "Loss"
 
# Bind the objects together
forest_change <- rbind(NP_forest_gain, NP_forest_loss)

我们可以制作一个条形图来可视化我们四个研究地点在 2000 年至 2016 年间失去和增加的森林覆盖量。因为一个更大的国家公园可以仅仅因为它更大而失去更多的森林（即，有更多的森林可以失去），我们可以将森林变化可视化为公园总面积的百分比。我们通过指定在下面的代码中执行此操作y = sum/GIS_AREA（或者您可以在数据框中创建一个新列，如果您愿意，可以在其中计算这些百分比）。

我们之前选择的ggthemr主题使图表更具信息图表的感觉。如果您需要更多标准格式，您可以添加+ theme_bw或添加+ theme_classic()到您的条形图代码。

(forest_barplot <- ggplot(forest_change, aes(x = NAME, y = sum/GIS_AREA, 
                                             fill = fct_rev(type))) +
    geom_bar(stat = "identity", position = "dodge") +
    labs(x = NULL, y = "Forest change (% of park area)\n") +
    # Expanding the scale removes the emtpy space below the bars
    scale_y_continuous(expand = c(0, 0)) +
    theme(text = element_text(size = 16),  # makes font size larger
          legend.position = c(0.1, 0.85),  # changes the placement of the legend
          legend.title = element_blank(),  # gets rid of the legend title
          legend.background = element_rect(color = "black", 
                                           fill = "transparent",   # removes the white background behind the legend
                                           linetype = "blank")))

请注意，将整个 ggplot 代码放在方括号 () 中会创建绘图，然后在绘图查看器中显示它。如果您没有括号，则您只是创建了对象，但还没有将其可视化。然后，您必须调用该对象，以便forest_barplot在您创建“forest_barplot”对象后只需键入即可显示该对象。

我们可以使用该ggsave函数来保存我们的图形。该文件将保存到您的工作目录所在的任何位置，您可以通过getwd()在控制台中运行来检查。

ggsave(forest_barplot, filename = "forest_barplot.png",
       height = 5, width = 7)

现在我们可以看到在我们感兴趣的保护区内获得和损失了多少森林，我们可以回到我们最初的研究问题，即森林变化在保护区之间有何不同，我们可以看看我们是否可以发现任何模式 -是否有任何类型的保护区更有可能失去森林？