【开GEE倒车】如何按着Google Earth Engine的脑袋搞逐像元计算

AUSTIN

研究方向真空球形鸡

需求

Google Earth Engine（GEE）是一款新兴的地理信息云计算平台。它的推出，说是为地学数据的处理计算带来了革命性的改观也不为过。从源数据到生产产品一切过程集成在云端，云的便捷性吸引了众多研究者跃跃欲试。GEE所设计的计算方式是完全的并行计算，这使得计算效率得到极大提升。

为了适应GEE并行计算的工作方式，我们需要在编程上进行一系列调整，如图像处理时，创建和调用的对象一般都应为整幅图片（ee.Image）或图像集（ee.ImageCollection），而不是每个像元。不同于图像处理的一般功能，GEE的图像上不提供按像元编号的索引。尽管后台处理时像元仍是基本处理单元，但像元并不显式地出现在代码中。按行列号索引处理像元，这种在MATLAB、R等工具中非常常见的编程思路，在GEE中是不可行的。

这种规定方便了后台的并行计算，但有时会对用户自定义部署算法带来困难。

在GEE上处理栅格数据，有丰富的简单函数可供调用和组合，实现覆盖大多数任务目标的算法。在简单的图像处理函数以外，GEE为一些特定的需求，如时间序列分析，提供了有限的内置算法（ee.Algorithm）。

然而，当用户希望求解某些GEE不提供内置函数的复杂问题时（例如正演模拟问题），不能进行像元索引的限制迫使用户必须硬着头皮，将算法改写为建立在GEE的对象之上的代码格式。其实这种情况下GEE提供的回旋余地还算勉强够用。但当我们进一步地，需要利用正演模拟模型来构建反演模型时（例如求解参数最优化问题），GEE的对象不够灵活的问题彻底突显，靠基于ee.Image和ee.ImageCollection转写算法的路子几乎行不通。

因为我的研究涉及较为复杂的生态模拟问题，随便一个生态模型都需要考虑众多的生物物理和生物化学过程，并且因此我迫切地需要在像元级别实现正向的模拟算法，以及反向的算法参数优化。GEE的免费计算资源实在是太香了，叫我不能随随便便就弃用，于是我开始探索一种在像元级别进行处理、且同时满足GEE并行计算编程格式的技术路线。

我承认这种做法可以说是相当的另类。放弃掉GEE推荐的图像处理编程方式，意味着会损失不少计算效率。但是我在改写算法的路上已经是走投无路，外加打小就蔫淘，另类不另类的那是后话，能跑起来就行了。

实现

此路不通我就绕路走。我们考虑将图像转成一个要素集（ee.FeatureCollection），其中容纳由每个像元转换而成的矢量点对象（ee.Geometry.Point）。既然图像对象（ee.Image和ee.ImageCollection）不支持索引像元，那么我们就利用 栅格转要素 的GEE内置函数，让像元所包含的数据转移到对应位置的矢量点上，由矢量点来代替像元作为计算的基本单元。在借助矢量点完成处理后，再对更新了运算结果的矢量点阵逆向进行一遍 要素转栅格 的过程，还原为图像进行输出。

下图是这个思路的概念图。需要注意的是，数据处理的算法是部署在单个矢量点之上的，所以这套方法与其说是“逐个像元作处理“，不如严谨一点地说是“逐个像元位置作处理”。

步骤1: 矢量点阵会继承图像每个像元的中心位置和数据，步骤2: 数据处理算法在每个矢量点上对数据做运算，步骤3: 用处理后的点阵重建图像

以下GEE JavaScript代码提供了这种思路下的一个示例。本例调用SRTM数字高程数据，在一个小块研究区上实现一个简单的数值筛选：高程若小于等于120米则置为零，若大于120米则保留原值。虽然对于这一滤波运算使用图像掩膜的标准方法会更有效率，但我们仅以此例示范逐像元处理的技术路线。该滤波操作没有实际应用意义，仅用于示范。

// 作者：知乎用户Austin (https://www.zhihu.com/people/zpzuo)
// 单位：波士顿大学
// 代码最后修改时间：2021年5月29日
// 开放一切学术交流使用，转载请注明出处
// 本代码是使用内置函数 .sample() 和 .reduceToImage() 在一个小块
// 示范区域内进行图像和点要素集之间无损、无偏转换的样例。数据存储
// 为点要素集的阶段进行了一次数值筛选操作，用以示范GEE逐像元处理
// 的技术路线。
// 加载感兴趣区域roi和数字高程数据
// 一块100余个像点的小区域
var roi = ee.Geometry.Polygon(
        [[[-71.65781521358379, 42.881740327133244],
          [-71.65591084518321, 42.87907123590606],
          [-71.65169977702983, 42.88032128122396],
          [-71.65436589279064, 42.88199189973334]]]);
// SRTM数字高程（DEM）数据
var image = ee.Image("USGS/SRTMGL1_003"); // 波段名为 elevation
// 利用 .sample() 和 .map() 实现栅格转点阵并且传递像元值
var points = image.addBands(ee.Image.pixelLonLat()) // 为图像加入像元经纬度两个层
  // 用 .sample() 将 ee.Image 转为 ee.FeatureCollection
  .sample({ 
    region: roi, // 以roi为采样区，不设置其它参数则采样全部像元
    geometries: true // 使生成的点阵中的每个点都置于像元中心
  // 用 .map() 对要素集中所有的要素逐个处理
  .map(function(sample){
    var lat = sample.get('latitude'); 
    var lon = sample.get('longitude');
    var value = sample.get('elevation'); 
    // 硬定义返回的要素为矢量点，点的位置由经纬度创建，数据为上一行的value
    return ee.Feature(ee.Geometry.Point([lon, lat]), {value: value}); 
// 对点阵中的每个点逐个进行数值筛选操作
// 用 .map() 对要素集中所有的要素逐个处理
var ptsComputed = points.map(function(point){
  // 获取点要素存储的原数据
  var attribute = ee.Number(point.get('value')); 
  // 数据小于等于120则重置为0
  var computedAttribute = ee.Algorithms.If(attribute.lte(120), 0, attribute); 
  // 修改点要素的数据并返回
  return point.set('value', computedAttribute); 
// 将滤波后的点阵重新转为图像，使用 .reduceToImage() 实现要素转栅格
var image2 = ptsComputed.reduceToImage({
  properties: ['value'],
  reducer: ee.Reducer.mean() // 这里随意选择一个reducer，因为不存在要素重叠所以没有区别
// 图像重采样，若无这一步图像将缺失投影而无法显示
.reproject(image.projection().crs(), image.projection().getInfo().transform )
// 重命名波段
.rename('computedAttribute');
// 可视化参数
var band_viz = {
  min: 90,
  max: 125,
  palette: ['0D0887', '5B02A3', '9A179B', 'CB4678',
            'EB7852', 'FBB32F', 'F0F921']};