谷歌地球引擎Google Earth Engine针对不同地表类型分别自动生成随机采样点的方法

news2024/9/25 1:22:08

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

本文是谷歌地球引擎（Google Earth Engine，GEE）系列教学文章的第22篇，更多GEE文章请参考专栏：GEE学习与应用（https://blog.csdn.net/zhebushibiaoshifu/category_11081040.html）。

首先，来具体明确一下本文的需求。我们现在有一个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参数——这个参数越大，对应的抽样点个数就越多；反之，抽样点个数则越少。大家可以结合实际情况，自行调整这个参数，直到能够获取到合适数量的抽样点。