谷歌地球引擎GEE在各地表类型中分别采样

本文介绍在 谷歌地球引擎 （Google Earth Engine， GEE ）中，按照给定的 地表分类数据 ，对 每一种不同的地物类型 ，分别加以全球范围内 随机抽样点 自动批量选取的方法。 本文是谷歌地球引擎（Google Earth Engine， GEE ）系列教学文章的第 22 篇 。

首先，来具体明确一下本文的需求。我们现在有一个 MODIS 的 地表分类 数据产品 MCD12Q1 ，其 LC_Type3 波段可以提供全球陆地范围内，每一年的 地表植被分类 数据，空间分辨率为 500 m 。如下图所示，就是其对于地表植被类型的分类；紫色框内的 8 种地表类型，就是 8 种不同的植被类型。

我们希望实现的是，在全球陆地范围内，针对上述 8 种植被类型，对于 每一种植被类型 ，分别筛选 3000 个左右的随机采样点；因为一共是 8 种植被类型，所以相当于我们最终希望得到的全部采样点个数大约为 8 * 3000 = 24000 。其中，每一种植被类型对应的 3000 个左右的采样点，都单独作为一个 FeatureCollection ，然后保存到 Assets 中。

明确了需求，即可开始代码的撰写。本文用到的代码如下。

var modis_type = ee.Image("projects/ee-ucanuse6/assets/vegetation_type"),
    rectangle = 
    /* color: #98ff00 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[-164.97163926414936, 74.9796307938761],
          [-164.97163926414936, -57.6624255131542],
          [176.74711073585064, -57.6624255131542],
          [176.74711073585064, 74.9796307938761]]], null, false);

// Map.addLayer(rectangle);
var landcoverTypes = [1, 2, 3, 4, 5, 6, 7, 8];

var modis_type_vis = {
  min: 1.0,
  max: 17.0,
  palette: [
    '05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff'
  ],
};
Map.addLayer(modis_type, modis_type_vis, "MODIS_Type");

landcoverTypes.forEach(function(type) {
  var type_mask = modis_type.eq(type);
  var type_image = modis_type.updateMask(type_mask);
  var type_point = type_image.sample({
    scale: 500,
    region: rectangle,
    numPixels: 140000,
    seed: 7,
    dropNulls: true,
    geometries: true
  });
  print(type_point);
  Map.addLayer(type_point, {}, "Point");
  
  Export.table.toAsset({
    collection: type_point,
    description: 'point_' + type,
    assetId: 'point_' + type,
  });
});

其中， // Map.addLayer(rectangle); 这句代码以上的内容，是获取我这里处理好的 MODIS 地表分类数据，以及我手动绘制的一个包含了全球主要陆地部分的 Polygon 。这一部分的代码，大家复制之后，可以转换成 Imports 的格式，如下图红色框内所示。

代码接下来的部分的含义如下。

首先， var landcoverTypes = [1, 2, 3, 4, 5, 6, 7, 8]; 定义了一个包含地物类型的数组 landcoverTypes ，这里的值 1 到 8 分别对应着我们所需的 8 种土地植被覆盖类别。

随后， var modis_type_vis = { ... }; 这部分代码定义了 MODIS 地物类型数据的可视化参数，其中包含了最小值、最大值和调色板信息，用于根据类型值将地物类型数据映射到不同颜色，从而可以更好地可视化 MODIS 地表分类数据产品。随后， Map.addLayer(modis_type, modis_type_vis, "MODIS_Type"); 将 MODIS 地表分类数据产品添加到地图中，并使用上述定义的可视化参数进行显示。

接下来， landcoverTypes.forEach(function(type) { ... }); 表示对 landcoverTypes 数组中的每个值进行循环，也就是对 8 种地表类型数据分别加以循环。

随后， var type_mask = modis_type.eq(type); 表示根据当前循环的类型值，创建一个布尔类型的掩膜 type_mask ，用于筛选出与该类型相符的像素。其次， var type_image = modis_type.updateMask(type_mask); 表示使用掩膜 type_mask 对地表类型数据进行屏蔽（掩膜），保留与当前类型相符的像素，并将其他像素设置为 null 。

其次， var type_point = type_image.sample({ ... }); 表示对掩膜后的地表类型数据进行抽样，生成样本点；在这里，我们使用给定的参数（比例尺、区域、像素数等）指定抽样的方式。

最后，使用 print(type_point); 打印该地表类型的样本点信息，以查看抽样结果；通过 Map.addLayer(type_point, {}, "Point"); 将样本点图层添加到地图中进行可视化。

完成上述操作后，通过 Export.table.toAsset({ ... }); 将样本点导出，将其存储为 GEE 中的 Asset 。

运行上述代码，首先我们看一下 Console 一栏的结果——也就是 print(type_point); 这句代码，打印出的 8 种地表类型对应的抽样点信息；如下图所示。

在这里，首先可以看到，由于部分地表类型数据对应的抽样点个数比较多，所以对于该类型（如上图中第一种和第四种植被类型）抽样点， print(type_point); 这句代码会报错，是打印不出来结果的。

其次，就是有一点，我这里这个代码暂时 不能直接限定每一种植被类型对应的抽样点个数 ——如上图所示，有的植被类型抽样点可以多到 print(type_point); 这句代码直接报错（也就是超过了 5000 个），而有的植被类型抽样点只有 248 个。这个问题，是 var type_point = type_image.sample() 这句代码中， .sample() 函数不能很好地限制抽样点个数导致的。我的理解是，由于全球范围内，不同植被类型对应的面积不一样，也就是不同植被类型对应的 MODIS 地表分类数据产品的像元总数也就不一样；但是我们在执行 .sample() 函数时，其抽样比例是固定的（例如每 3000 个像素中，选出 1 个作为抽样点）——那么对于在全球 面积占比大 ，也就是 像元个数多 的植被类型，其得到的 抽样点个数自然就多 ；反之自然就少。针对这种情况，如果大家希望对于每一种植被类型，都获取比较一致的抽样点个数，那么就可以单独对这种（或者这几种）植被类型运行上述代码（即修改 landcoverTypes 中的类型数值），然后调整 .sample() 函数中的 numPixels 参数——这个参数越大，对应的抽样点个数就越多；反之，抽样点个数则越少。大家可以结合实际情况，自行调整这个参数，直到能够获取到合适数量的抽样点。