How it works(24) Geotrellis是如何在Spark上计算的(E) 地图代数之游标类算子

1. 引入

上一章我们研究了Geotrellis中Local算子的实现.Local类算子种类较多,实现也较简单,其核心为:Tile中某个位置的值由另一个Tile中相同位置的值通过计算得到.

Geotrellis中也存在更加复杂的Focal类算子,其核心为:Tile中某个像元位置的值由该像元周围若干像元(邻域)的值经过特定计算得到.可见,Focal类算子有两个关键点:

• 定义邻域
• 定义如何根据邻域内容计算结果

我们首先看邻域是如何定义的.

2. 邻域定义

在Geotrellis中,邻域描述了某个像元周围哪些位置的像元参与到计算中.这些周围区域的形状可能简单也可能复杂,为了不针对各种形状制定复杂的结构从而影响计算效率,Geotrellis通过这样一种方式描述邻域:

代码位于[geotrellis/raster/mapalgebra/focal/Neighborhood.scala]

trait Neighborhood extends Serializable {
  // 聚焦范围,范围为1则对应3x3的聚焦邻域,即向中心点周围8方向各扩展1个单位
  val extent:Int 

  // 是否具有遮罩(只有最基础的正方形邻域没有遮罩)
  val hasMask:Boolean

  /** 
   * 定义具体的邻域:
   * 通过定义不同位置的返回值确定邻域形状,返回true代表被遮罩覆盖,不参与计算
   * (0,0) 代表左上角, (extent*2+1,extent*2+1) 代表右下角
   */
  def mask(col:Int,row:Int):Boolean = { false }
}

Geotrellis定义了5种邻域(#符号表示在邻域中最终参与计算,.符号表示不参与计算,extent设为4):

正方形(Square)邻域      十字形(Nesw)邻域        圆形(Circle)邻域

# # # # # # # # #       . . . . # . . . .       . . . . # . . . .
# # # # # # # # #       . . . . # . . . .       . . # # # # # . . 
# # # # # # # # #       . . . . # . . . .       . # # # # # # # . 
# # # # # # # # #       . . . . # . . . .       . # # # # # # # .
# # # # # # # # #       # # # # # # # # #       # # # # # # # # #
# # # # # # # # #       . . . . # . . . .       . # # # # # # # .
# # # # # # # # #       . . . . # . . . .       . . # # # # # . .
# # # # # # # # #       . . . . # . . . .       . . # # # # # . .
# # # # # # # # #       . . . . # . . . .       . . . . # . . . . 

环形(Annulus)邻域       90°-180°楔形(Wedge)邻域   180°-90°的楔形邻域

. . . . # . . . .       . . . . # . . . .       . . . . # . . . . 
. . # # # # # . .       . . # # # . . . .       . . . . # # # . . 
. # # # # # # # .       . # # # # . . . .       . . . . # # # # . 
. # # . . . # # .       . # # # # . . . .       . . . . # # # # . 
# # # . . . # # #       # # # # # . . . .       # # # # # # # # # 
. # # . . . # # .       . . . . . . . . .       . # # # # # # # . 
. # # # # # # # .       . . . . . . . . .       . # # # # # # # . 
. . # # # # # . .       . . . . . . . . .       . . # # # # # . . 
. . . . # . . . .       . . . . . . . . .       . . . . # . . . . 

了解完邻域的定义,我们就需要看看Focal类算子如何根据邻域中的内容计算结果.

3. Focal类算子的整体架构

Focal算子可以划分为4大类:

• CursorCalculation(游标类算子):
  • Max(最大值)/Min(最小值)
  • StandardDeviation(标准差)
  • Mean(均值)/Median(中值)/Mode(众数)/sum(总和)/moran(莫兰指数)
• KernelCalculation(核算子)
  • Convolve(卷积)
• CellwiseCalculation(像元尺度算子)
  • Conway(生命游戏)
  • Mean(均值)/Median(中值)/Mode(众数)/sum(总和)/moran(莫兰指数)
• SurfacePointCalculation(表面点算子)
  • aspect(坡向)
  • slope(坡度)
  • hillshade(山体阴影)

其中部分统计算子同时具有游标类算子和像元尺度算子两种实现.

这4大类算子都继承自抽象的Focal算子类(FocalCalculation):

代码位于[geotrellis/raster/mapalgebra/focal/FocalCalculation.scala]

trait Resulting[T] {
  val copyOriginalValue: (Int, Int, Int, Int) => Unit
  def result: T
}

abstract class FocalCalculation[T](
    val r: Tile, n: Neighborhood, analysisArea: Option[GridBounds[Int]], val target: TargetCell)
  extends Resulting[T]
{
  // 计算范围.默认为Tile自身的范围 
  val bounds: GridBounds[Int] = analysisArea.getOrElse(GridBounds(r))

  def execute(): T
}

其中包含了全部Focal算子需要具备的三个核心:

• 输入:定义了支持的输入参数
• 运算:约定实现execute方法
• 输出:约定实现对结果集操作的方法
  • 不同的算子按需混入不同类型的Resulting特质,实现操作具体数据类型的结果集
  • 我们以最具有一般性的ArrayTileResult为例(ArrayTileResult适用于任何不需要改变原始数据类型的算子):

// 通过自类型注释(self-type annotation),限定自己必须是FocalCalculation类型
// 保证其中具有r,bounds等对象
trait ArrayTileResult extends Resulting[Tile] { self: FocalCalculation[Tile] =>
  def resultCellType: DataType with NoDataHandling = r.cellType
  val cols: Int = bounds.width
  val rows: Int = bounds.height
  // 定义返回的对象,类型为传入tile对象的数据类型
  val resultTile: MutableArrayTile = ArrayTile.empty(resultCellType, cols, rows)
  val copyOriginalValue: (Int, Int, Int, Int) => Unit =
    if(!r.cellType.isFloatingPoint) {
      { (focusCol: Int, focusRow: Int, col: Int, row: Int) =>
        resultTile.set(col, row, r.get(focusCol, focusRow))
      }
    } else {
      { (focusCol: Int, focusRow: Int, col: Int, row: Int) =>
        resultTile.setDouble(col, row, r.getDouble(focusCol, focusRow))
      }
    }
  // 将结果集指向约定值,在算子中我们操作的也是这个值
  // 使用def预定义再通过后期指向的方式是为了通过编译
  def result = resultTile
}

4大类中最关键的算子是游标类算子(CursorCalculation),核算子是其特化,而像元尺度算子和表面点算子可以看做其简化.因此我们从游标类算子开始分析.

4. 游标类算子-以Max算子为例

4.1 抽象游标类算子的定义

首先来分析抽象的游标类算子是如何定义的:

代码位于[geotrellis/raster/mapalgebra/focal/FocalMethods.scala]

abstract class CursorCalculation[T](tile: Tile, n: Neighborhood, val analysisArea: Option[GridBounds[Int]], target: TargetCell)
  extends FocalCalculation[T](tile, n, analysisArea, target)
{
  def traversalStrategy = TraversalStrategy.DEFAULT

  def execute(): T = {
    // 1.
    val cursor = Cursor(tile, n, bounds)
    // 2.
    val calcFunc: () => Unit =
      target match {
        case TargetCell.All =>
          { () => calc(tile, cursor) }
        case TargetCell.Data =>
          { () =>
            calc(tile, cursor)
            if(!isData(r.get(cursor.focusCol, cursor.focusRow))){
              copyOriginalValue(cursor.focusCol, cursor.focusRow, cursor.col, cursor.row)
            }
          }
        case TargetCell.NoData =>
          { () =>
            calc(tile, cursor)
            if(!isNoData(r.get(cursor.focusCol, cursor.focusRow))) {
              copyOriginalValue(cursor.focusCol, cursor.focusRow, cursor.col, cursor.row)
            }
          }
      }
    // 3.
    CursorStrategy.execute(cursor, calcFunc, bounds, traversalStrategy)
    result
  }

  def calc(tile: Tile, cursor: Cursor): Unit
}

游标算子的核心是execute方法,在该方法中主要执行了如下3步:

1. 定义一个游标对象(Cursor)
2. 定义calcFunc方法
3. 将一系列参数传入CursorStrategy.execute方法中执行

我们按顺序先从游标对象的定义开始分析.

4.2 游标对象的定义

Focal类算子是一种滑动窗口型算子,而滑动的那个"窗口"就是游标对象,它定义了这个窗口的大小与形状,是我们进行滑动和读取数据的基本粒度.

4.2.1 游标对象的初始化与滑动

游标对象进行滑动前首先需要将其初始化:

代码位于[geotrellis/raster/mapalgebra/focal/Cursor.scala]

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {
 
 // 记录总的行列数
 private val rows = r.rows
 private val cols = r.cols

 // 偏移量,如果不指定限制范围,则偏移量为0
 val analysisOffsetCols = analysisArea.colMin
 val analysisOffsetRows = analysisArea.rowMin
 
 // 游标的边长
 private val d = 2 * extent + 1

 // 记录当前游标所在区域的行列范围
 private var _colmin = 0
 private var _colmax = 0
 private var _rowmin = 0
 private var _rowmax = 0

 // 相当于对上面的min/max定义了一个公开的"get"方法
 protected def colmin = _colmin
 protected def colmax = _colmax
 protected def rowmin = _rowmin
 protected def rowmax = _rowmax

 // 记录在当前移动过程中,加入的/移除的行列号
 private var addedCol = 0
 private var removedCol = 0
 private var addedRow = 0
 private var removedRow = 0
   
 // 当前移动方向
 var movement = NoMovement

 // 当前聚焦位置及其"get"方法
 private var _col = 0
 private var _row = 0
 def focusCol = _col
 def focusRow = _row

 def isReset = movement == NoMovement

 // 当前聚焦的行列号的真实值
 def col = _col - analysisOffsetCols
 def row = _row - analysisOffsetRows

 // 设置起点
 def centerOn(col: Int, row: Int) = {
   movement = NoMovement
   // 设置当前要处理的行列号为起点行列号
   _col = col
   _row = row
   // 计算在该起点下,最大最小行列号都是多少
   _colmin = max(0, _col - extent)
   _colmax = min(cols - 1, _col + extent)
   _rowmin = max(0, _row - extent)
   _rowmax = min(rows - 1, _row + extent)
 }
 
 // ... 省略其他
}

初始化过程实现了三个功能:

• 限定行列范围
• 设定滑动起点
• 定义辅助变量

我们以指定参数[Tile=5x4,extent=1,不指定分析范围,初始位置为(0,0)]为例:

初始化后各个变量为:

变量	初始值	变量	初始值
cols	5	rows	4
analysisOffsetCols	0	analysisOffsetRows	0
_colmin	0	_rowmin	0
_colmax	1	_rowmax	1
_col	0	_row	0
addedCol	0	addedRow	0
removedCol	0	removedRow	0

初始化完成后就可以进行滑动(move)操作了:

代码位于[geotrellis/raster/mapalgebra/focal/Cursor.scala]

sealed trait Movement { val isVertical: Boolean }
// 定义滑动方向标志位
object Movement {
 val Up = new Movement { val isVertical = true }
 val Down = new Movement { val isVertical = true }
 val Left = new Movement { val isVertical = false }
 val Right = new Movement { val isVertical = false }
 val NoMovement = new Movement { val isVertical = false }
}

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {
   // ... 省略初始化阶段
   
   // 滑动游标
   def move(m: Movement) = {
       movement = m
       m match {
         case Up =>
           addedRow = _rowmin - 1
           removedRow = _row + extent
           _row -= 1
         // ... 省略类似的向下移动逻辑
         case Left =>
           addedCol = _colmin - 1
           removedCol = _col + extent
           _col -= 1
         // ... 省略类似的向右移动逻辑
         case _ =>
       }
       
       // 重新计算滑动后的新值
       _colmin = max(0, _col - extent)
       _colmax = min(cols - 1, _col + extent)
       _rowmin = max(0, _row - extent)
       _rowmax = min(rows - 1, _row + extent)
 }
   
   // ...省略其他
 
}

我们以向右滑动为例:

此时对比初始值:

变量	初始值	滑动后
_col	0	1
_colmax	1	2
addedCol	0	2
removedCol	0	-1(无意义)

这些变量都是为了方便进行遍历操作而准备的.

4.2.2 简单邻域下的遍历操作

遍历操作使我们可以获取该游标中的每一个像元的值.在上文中我们介绍的几种邻域中,最简单的是正方形邻域,因为没有需要跳过的位置,因此处理起来也最简单.我们以正方形邻域为例,看看遍历操作的一般逻辑.

针对游标对象,有两类遍历方式:

• 完全遍历:遍历游标中的全部位置.
• 变化区域遍历:只遍历上一次滑动操作后新加入/退出的行列,这在某些算子中可以避免重复计算

核心代码如下:

代码位于[geotrellis/raster/mapalgebra/focal/Cursor.scala]

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {
    
    // ... 省略初始化阶段
    
    // 遍历全部位置
    protected def foreach(f: (Int, Int)=>Unit): Unit = {
          var y = _rowmin
          var x = 0
          // 遍历全部的行列
          while(y <= _rowmax) {
            x = _colmin
            while(x <= _colmax) {
              f(x, y)
              x += 1
            }
            y += 1
          }
        }
    
    // 只遍历移动后新加入的位置
    protected def foreachAdded(f: (Int, Int)=>Unit): Unit = {
        // 如果不移动,则与foreach无异
        if(movement == NoMovement) {
          foreach(f)
        } else if (movement.isVertical) {
          // 垂直方向移动后,只遍历新加入的那行
          // 排除无意义的addedRow值
          if(0 <= addedRow && addedRow < rows) {
              var x = _colmin
              while(x <= _colmax) {
                f(x, addedRow)
                x += 1
              }
          }
        } else { // 水平方向移动后,只遍历新加入的那列
          if(0 <= addedCol && addedCol < cols) {
              var y = _rowmin
              while(y <= _rowmax) {
                f(addedCol, y)
                y += 1
              }
          }
        }
      }
  
    // ... 省略与foreachAdded相似的foreachRemoved
    
    // ... 省略其他

4.2.3 复杂邻域下的遍历操作与游标遮罩的定义

游标遮罩(CursorMask)是一个辅助遍历器,辅助游标遍历非正方形的复杂邻域.

相比正方形邻域,复杂邻域需要面对的问题就是其中有若干区域在遍历中应被遮罩遮蔽,不参与计算.比较笨的方法是每个位置都进行判断,但必然会产生大量重复计算.为避免浪费计算资源,游标遮罩采取空间换时间的策略,设计了"快查表"机制,尽可能在不占用很多额外空间的前提下只计算需要计算的位置.

在当前语境下,“遮罩"等同于"不参与计算”,“非遮罩"等同于"参与计算”

4.2.3.1 "快查表"的数据结构

游标遮罩为"快查表"设计了一种便于计算的用一维数组表示二维数据的数据结构:遮罩集合(MaskSet):

代码位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

private class MaskSet {
    // d=2*extent+1,与Cursor中意义一致
    private var data = Array.ofDim[Int](d*(d+1))
    // 插入一个值的时候需要更新两个值
    def add(i:Int,v:Int) = {
      val len = data(i*(d+1)) + 1
      data(i*(d+1) + len) = v
      data(i*(d+1)) = len
    }
    // 遍历时需要指定所在行
    def foreachAt(i:Int)(f:Int=>Unit) = {
      val len = data(i*(d+1))
      if(len > 0) {
        var x = 0
        while(x < len) {
          f(data(i*(d+1) + x + 1))
          x += 1
        }
      }
    }
}

遮罩集合虽然使用1维数组存储数据,但它与游标的尺寸是对应的.一个尺寸为3x3的游标,对应的遮罩集合中data数组的长度就是3x(3+1)=12,即我们可以将其看做一个3x4的二维矩阵叠放到一个一维数组中,额外申请的一列用来存储每一行中元素的数量.这种数据结构的好处是:在一个数据稀疏的矩阵中,可以避免大量不必要的判断.

下面我们来分析在游标遮罩对象中,各种"快查表"都是如何被初始化的.

4.2.3.2 "快查表"的初始化

代码位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

class CursorMask(d:Int,f:(Int,Int)=>Boolean) {
    // ...省略MaskSet的定义
    
    // 下文中的最左/最右/位置等都是指在游标范围内的相对位置
    // 需要通过反算得到在游标移动后在tile中的真实位置
    
    // 记录所有参与计算的非遮罩位置
    private val unmasked = new MaskSet
    
    // 记录最左列/最右列哪一行参与计算,因为只有1列,所以无需定义为遮罩集合的形式
    private var westColumnUnmasked = Array.ofDim[Int](d+1)
    private var eastColumnUnmasked = Array.ofDim[Int](d+1)
    
    // 记录向某个方向移动后改变遮罩性质的全部位置
    private val maskedAfterMoveLeft = new MaskSet
    private val unmaskedAfterMoveLeft = new MaskSet
    private val maskedAfterMoveUp = new MaskSet
    private val unmaskedAfterMoveUp = new MaskSet
    
    var x = 0
    var y = 0
    var len = 0
    var isMasked = false
    var leftIsMasked = false
    // 遍历游标中的全部位置
    while(y < d) {
        x = 0
        len = 0
        while(x < d) {
          isMasked = f(x,y)
          if(!isMasked) {
            // 记录所有参与计算的非遮罩位置
            unmasked.add(y,x)
            // 记录最左/最右两列的信息,因为他们比较特殊,所以特别记录
            if(x == 0) {
              val len = westColumnUnmasked(0) + 1
              westColumnUnmasked(len) = y
              westColumnUnmasked(0) = len
            }
            if(x == d - 1) {
              val len = eastColumnUnmasked(0) + 1
              eastColumnUnmasked(len) = y
              eastColumnUnmasked(0) = len
            }
          }
        
          // 记录随着向左移动,是否发生遮罩状态变化
          if(x > 0) {
            val leftIsMasked = f(x-1,y)
            // 向左移动将变为非遮罩状态
            if(leftIsMasked && !isMasked) {
              unmaskedAfterMoveLeft.add(y,x)
            } else if(!leftIsMasked && isMasked) {
              maskedAfterMoveLeft.add(y,x)
            }
          }
          // ... 省略类似的向上移动逻辑
          x += 1
        }
        y += 1
    }
    
    // ... 省略其他
}

我们以十字型邻域(3x3)为例:

 . # .
 # # #
 . # .

初始化后,各个快查表的内容为:

unmasked
[1, 1, null, null, 3, 0, 1, 2, 1, 1, null, null]
westColumnUnmasked
[1, 1, 0, 0]
eastColumnUnmasked
[1, 1, 0, 0]
maskedAfterMoveUp
[null, null, null, null, null, null, null, null, 2, 0, 2, null]
unmaskedAfterMoveUp
[null, null, null, null, 2, 0, 2, null, null, null, null, null]
maskedAfterMoveLeft
[1, 2, null, null, null, null, null, null, 1, 2, null, null]
unmaskedAfterMoveLeft
[1, 1, null, null, null, null, null, null, 1, 1, null, null]

这些快查表中内容现在看比较抽象,我们将在遍历操作中看看它们如何起作用.

4.2.3.3 使用快查表遍历数据

与简单邻域的遍历一样,复杂邻域也有完全遍历和变化区域遍历两种.我们先从完全遍历开始:

代码位于[geotrellis/raster/mapalgebra/focal/Cursor.scala]

object Cursor {
  def apply(r: Tile, n: Neighborhood, analysisArea: GridBounds[Int]): Cursor = {
    val result = new Cursor(r, analysisArea, n.extent)
    // 复杂邻域时设置游标遮罩对象
    if(n.hasMask) { result.setMask(n.mask) }
    result
  }
}

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {

    // ... 省略
    
    def setMask(f: (Int, Int) => Boolean) = {
        hasMask = true
        mask = new CursorMask(d, f)
    }
  
    protected def foreach(f: (Int, Int)=>Unit): Unit = {
        if(!hasMask) {
            // ... 省略简单邻域的情况
        } else {
            var y = 0
            while(y < d) {
                // 遍历某一行中的全部列中参与计算的位置
                mask.foreachX(y) { x =>
                  // 反算回真实位置
                  val xTile = x + (_col - extent)
                  val yTile = y + (_row - extent)
                  // 避免边缘位置越界
                  if(_colmin <= xTile && xTile <= _colmax && _rowmin <= yTile && yTile <= _rowmax) {
                    f(xTile, yTile)
                  }
                }
                y += 1
            }
        }
    }
    
    // ... 省略其他
}

核心为对游标遮罩对象foreachX方法的调用:

代码位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

class CursorMask(d:Int,f:(Int,Int)=>Boolean) {
    // ... 省略
    
    // 固定行,遍历列
    def foreachX(row:Int)(f:Int=>Unit) = {
        // 遍历全部非遮罩位置
        unmasked.foreachAt(row)(f)
    }
    
    // 省略
}

变化区域遍历则复杂许多,因为不仅要处理边缘部分的增减,还要处理因为遮罩区域变化引起的增减.我们以foreachAdded为例:

代码位于[geotrellis/raster/mapalgebra/focal/Cursor.scala]

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {

    // ... 省略

    protected def foreachAdded(f: (Int, Int)=>Unit): Unit = {
        // 先处理边缘部分
        // 垂直移动的情况
        if (movement.isVertical) {
          if(0 <= addedRow && addedRow < rows) {
            if(!hasMask) {
                // ... 省略简单邻域的情况
            } else {
              // addedRow - (_row - extent)=在Cursor中row的相对变化值
              mask.foreachX(addedRow - (_row - extent)) { x =>
                // 再从相对值反算回来
                val xTile = x + (_col - extent)
                if(0 <= xTile && xTile <= cols) {
                  f(xTile, addedRow)
                }
              }
              
              // 等同于下面的方法
              /*
              if(movement == Up) {
                mask.foreachX(0) { x =>
                  val xTile = x + (_col - extent)
                  if(0 <= xTile && xTile < cols) {
                    f(xTile, addedRow)
                  }
                }
              }
              else {
                mask.foreachX(d-1) { x =>
                  val xTile = x + (_col - extent)
                  if(0 <= xTile && xTile < cols) {
                    f(xTile, addedRow)
                  }
                }
              }
              */
            }
          }
        } else { // 水平移动的情况
          if(0 <= addedCol && addedCol < cols) {
            if(!hasMask) {
                // ... 省略简单邻域的情况
            } else {
              if(movement == Left) {
               // 向左移动,左侧的非遮罩位置会变为新加入位置
                mask.foreachWestColumn { y =>
                  // 反算回实际行
                  val yTile = y + (_row - extent)
                  if(0 <= yTile && yTile < rows) {
                    f(addedCol, yTile)
                  }
                }
              } else { // 向右移动
                // ...省略相似逻辑
              }
            }
          }
        }
        
        // 处理游标内部的其他变化区域
        if(hasMask) {
          mask.foreachUnmasked(movement) { (x, y) =>
            val xTile = x + (_col - extent)
            val yTile = y + (_row - extent)
            if(0 <= xTile && xTile < cols && 0 <= yTile && yTile < rows) {
              f(xTile, yTile)
            }
          }
        }
    }
    
    // ... 省略其他

}

核心为对游标遮罩foreachWestColumn和foreachUnmasked方法的调用:

代码位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

class CursorMask(d:Int,f:(Int,Int)=>Boolean) {
    // ... 省略
    
    def foreachWestColumn(f:Int=>Unit) = {
        val len = westColumnUnmasked(0)
        if(len > 0) {
          var y = 1
          // 遍历游标最左列中每一个非遮罩位置
          while(y <= len) {
            f(westColumnUnmasked(y))
            y += 1
          }
        }
    }
    
    private def foreach(xOffset:Int,yOffset:Int,startY:Int,set:MaskSet)(f:(Int,Int)=>Unit) = {
        var y = startY
        while(y < d) {
          set.foreachAt(y) { x => f(x - xOffset, y - yOffset) }
          y += 1
        }
    }
    
    def foreachUnmasked(mv:Movement)(f:(Int,Int)=>Unit): Unit = {
        mv match {
          case Left => foreach(0,0,0,unmaskedAfterMoveLeft)(f)
          case Right => foreach(1,0,0,maskedAfterMoveLeft)(f)
          case Up => foreach(0,0,1,unmaskedAfterMoveUp)(f)
          case Down => foreach(0,1,1,maskedAfterMoveUp)(f)
          case _ =>
        }
    }
    
    // ... 省略
}

变化区域遍历的流程是这样的:

1. 通过遍历"边缘快查表",得到最左边/最右边新增的像元位置
2. 通过遍历"变化快查表",得到其他位置新增的像元位置

为什么不能直接通过遍历"变化快查表"得到包括左右边缘的变化信息呢?因为当初定义"变化快查表"的时候,没法处理边缘区域这些移动后会越界的位置,因此无法记录边缘区域的变化.

在上面定义"快查表"的时候,我们可能有些疑惑:为什么只需要定义诸如unmaskedAfterMoveLeft之类向左/向上的方法,却不需要定义对应的unmaskedAfterMoveRight等向右/向下的方法呢?我们从foreachUnmasked方法中就可以看到答案:unmaskedAfterMoveRight与偏移一列后的maskedAfterMoveLeft表达了同样的意思.

至此,我们就大致了解了游标对象提供了哪些功能,以及实现原理,基于这些原理与机制,我们需要看看如何在这些基础之上,构建各种实际的算子.

4.3 calcFunc方法与Max算子

我们再回到抽象Focal类的execute方法中.在定义完游标对象后,又定义了calcFunc方法.

代码位于[geotrellis/raster/mapalgebra/focal/FocalCalculation.scala]

sealed trait TargetCell extends Serializable

object TargetCell {
  object NoData extends TargetCell
  object Data extends TargetCell
  object All extends TargetCell
}

// ... 省略 

val calcFunc: () => Unit =
  target match {
    case TargetCell.All =>
      { () => calc(tile, cursor) }
    case TargetCell.Data =>
      { () =>
        calc(tile, cursor)
        // 把原始的Nodata值复制回来,避免算子使Nodata值发生改变
        if(!isData(r.get(cursor.focusCol, cursor.focusRow))){
          copyOriginalValue(cursor.focusCol, cursor.focusRow, cursor.col, cursor.row)
        }
      }
    case TargetCell.NoData =>
      // ... 省略与TargetCell.Data类似的逻辑
}

// ... 省略

其中的核心是对calc方法的调用,它们对应着具体算子的实现.我们以Max算子为例:

代码位于[geotrellis/raster/mapalgebra/focal/Max.scala]

def calc(r: Tile, cursor: Cursor) = {
  var m = Int.MinValue
  cursor.allCells.foreach { (x, y) =>
    val v = r.get(x, y)
    if(v > m) { m = v }
  }
  // 注意这里的col/row (col = focusCol - analysisOffsetCols)
  // 即结果的写入永远从0,0开始,没有偏移量
  resultTile.set(cursor.col, cursor.row, m)
}

Max算子的实现很简单:遍历当前游标中的每一个像元,得到最大值,然后将最大值赋值给游标当前的聚焦位置.

calcFunc还受目标策略(TargetCell)影响,如果用户指定当前算子只针对非Nodata值或Nodata值,在计算完后还需要用旧值将原先的非Nodata值或Nodata值的位置恢复,避免算子使Nodata值发生改变,造成与预期不一样的结果.这部分的逻辑主要调用了copyOriginalValue方法,我们再来分析方法的实现:

代码位于[geotrellis/raster/mapalgebra/focal/FocalCalculation.scala]

val copyOriginalValue: (Int, Int, Int, Int) => Unit =
if(!r.cellType.isFloatingPoint) {
    { (focusCol: Int, focusRow: Int, col: Int, row: Int) =>
        resultTile.set(col, row, r.get(focusCol, focusRow))
    }
} else {
    // ... 省略类似的逻辑
}

传入的值分两种:focusCol/focusRow,col/row (col = focusCol - analysisOffsetCols),我们将其还原到calcFunc的场景(假定为非浮点类数值):

case TargetCell.Data =>
  { () =>
    calc(tile, cursor)
    if(!isData(r.get(cursor.focusCol, cursor.focusRow))){
      resultTile.set(cursor.col, cursor.row, r.get(cursor.focusCol, cursor.focusRow))
    }
  }

恢复Nodata值的逻辑是这样:

1. 检测当前游标聚焦点位置是否为Nodata值
2. 如果是Nodata值,则将其赋予给当前游标聚焦点减去分析区偏移的位置.这也与结果的写入逻辑相一致.

4.4 游标策略的定义

如果说游标对象和calc方法定义了游标在滑动到Tile中某一个位置时的静态行为,游标策略(CursorStrategy)则描述了游标如何在整个Tile上滑动的动态过程.因此在抽象游标算子的最后一步,会将一系列参数传入CursorStrategy.execute方法中执行.我们来分析CursorStrategy.execute方法的定义

代码位于[geotrellis/raster/mapalgebra/focal/FocalStrategy.scala]

// 定义了3种游标滑动策略
sealed trait TraversalStrategy
// 之字形
case object ZigZagTraversalStrategy extends TraversalStrategy
// 线型
case object ScanLineTraversalStrategy extends TraversalStrategy
// 螺旋型
case object SpiralZagTraversalStrategy extends TraversalStrategy
// 默认策略为之字形
object TraversalStrategy {
  def DEFAULT: TraversalStrategy = ZigZagTraversalStrategy
}

object CursorStrategy {
  
  // ... 省略
  
  def execute(
    cursor: Cursor,
    calc: () => Unit,
    analysisArea: GridBounds[Int],
    traversalStrategy: TraversalStrategy
  ): Unit = {
    traversalStrategy match {
      case ZigZagTraversalStrategy => handleZigZag(analysisArea, cursor, calc)
      case ScanLineTraversalStrategy => handleScanLine(analysisArea, cursor, calc)
      case SpiralZagTraversalStrategy => handleSpiralZag(analysisArea, cursor, calc)
    }
  }
  
  // ...省略其他
}

可见,execute方法的核心就是使用不同的策略操游标进行滑动的过程.我们以默认的之字形策略为例:

代码位于[geotrellis/raster/mapalgebra/focal/FocalStrategy.scala]

  private def handleZigZag(analysisArea: GridBounds[Int], cursor: Cursor, calc: () => Unit) = {
    val colMax = analysisArea.colMax
    val rowMax = analysisArea.rowMax
    val colMin = analysisArea.colMin
    val rowMin = analysisArea.rowMin

    var col = colMin
    var row = rowMin

    var direction = 1
    // 初始化游标的位置
    cursor.centerOn(col, row)

    while(row <= rowMax) {
      // 每滑动一次都进行一次计算,并记录结果
      calc()
      // 先向右滑动
      col += direction
      // 移动越界了,就开始换行并转向
      if(col < colMin || colMax < col) {
        // 滑动到正下方位置
        direction *= -1
        row += 1
        col += direction
        cursor.move(Movement.Down)
      } else {
        if(direction == 1) { cursor.move(Movement.Right) }
        else { cursor.move(Movement.Left) }
      }
    }
  }

不同策略的滑动轨迹如下:

之字形策略              线型策略                螺旋型策略
→   →   →   ↓           →   →   →   →           →   →   →   ↓
↓   ←   ←   ←           →   →   →   →           →   →   →   ↓
→   →   →   →           →   →   →   →           ↑   ←   ←   ←

5. 总结

我们再整体梳理一下Max算子:

代码位于[geotrellis/raster/mapalgebra/focal/FocalMethods.scala]

object Max {
  def calculation(tile: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): FocalCalculation[Tile] = {
    if (tile.cellType.isFloatingPoint) {
      new CursorCalculation[Tile](tile, n, bounds, target)
        with ArrayTileResult
      {
        def calc(r: Tile, cursor: Cursor) = {
            // ... 省略
        }
      }

    } else {
      // ... 省略处理Int类型数据时类似的逻辑
    }
  }

  def apply(tile: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): Tile =
    calculation(tile, n, bounds, target).execute()
}

1. Max算子接收原始Tile和邻域信息作为参数,输出结果Tile.本质上是对一个Focal类算子执行execute方法.
2. 根据最大值算法的特点,Max算子实现为Focal类算子的一个子类:游标类算子,并且根据其返回值的特性,混入了ArrayTileResult特质.
3. 为成为一个实际的游标类算子,Max实现了预定的calc方法,Max算法的核心就在这里.
4. 用户实际调用Max算子时,游标类算子会创建一个滑动游标,在用户传入的Tile上按照指定的策略滑动,每次滑动后都会计算当前滑动窗口内通过定义的calc方法得到的结果值,并存入结果Tile中.
5. 当滑动完成,用户会得到最终的结果Tile.

Max算子通过类似Local类算子的方法挂载:

代码位于[geotrellis/raster/mapalgebra/focal/FocalMethods.scala]

trait FocalMethods extends MethodExtensions[Tile] {
  // ... 省略
  def focalMax(n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): Tile = {
    Max(self, n, bounds, target)
  }
  // ... 省略
}

至此, 我们也大致了解了游标类算子的结构与运行流程.