[大数据]How it works(23) Geotrellis是如何在Spark上计算的(D) 地图代数之Local类算子

1. 引入

在上一节,我们使用Tile对象的Map方法实现了计算NDVI的功能.但对于一些更复杂的功能,Map(f)的操作无法胜任或实现起来极为麻烦.好在Geotrellis实现了较为丰富的地图代数功能,能为我们节省不少时间.

官方文档对此部分着墨不多,不过我们可以知道,内置的地图代数算子主要分为三大类:

• 局部运算(Local):作用于单个像元
• 焦点运算(Focal):作用于邻域内的所有像元
• 分区运算(Zonal):作用于区域内的所有像元

还有不属于三大类的水文相关算法(hydrology)

地图代数逻辑的实现位于geotrellis.raster.mapalgebra包中,操作类型为Tile对象,同时,如果操作的是Seq[(key,Tile)])或Rdd[(k,Tile)]对象,则分别位于geotrellis.layer.mapalgebra和geotrellis.spark.mapalgebra,实现了对基本逻辑的封装.

我们以最基础且简单的Abs(取绝对值)算子为例,分析Geotrellis如何实现算子逻辑,如何将其绑定到Tile,Seq[(k,Tile)]和Rdd[(k,Tile)]对象上,并从中推断出Local类算子的一般结构.

2. Abs算子

2.1 算子的实现

代码位于[geotrellis/raster/mapalgebra/local/Abs.scala]

// 将算子封装为单例对象并扩展Serializable特质可以方便的在多节点间传递
object Abs extends Serializable {
  // 对Tile的每个像素值调用abs方法
  def apply(r: Tile): Tile =
    r.dualMap { z: Int => if(isNoData(z)) z else z.abs }
              { z: Double => z.abs }
}

函数调用了Tile.dualMap,一个简单的柯里化函数,实现了按Tile数据类型选择适当的Map方法:

代码位于[geotrellis/raster/Tile.scala]

def dualMap(f: Int => Int)(g: Double => Double): Tile =
    if (cellType.isFloatingPoint) mapDouble(g) else map(f)

dualMap又是对map和mapDouble方法的封装.他们定义于具体的Tile类中,我们以较为简单的ArrayTile为例:

代码位于[geotrellis/raster/ArrayTile.scala]

def map(f: Int=>Int): Tile = {
    val output = ArrayTile.alloc(cellType, cols, rows)
    var i = 0
    val len = size
    while (i < len) {
      output(i) = f(apply(i))
      i += 1
    }
    output
}

def mapDouble(f: Double => Double): Tile = {
    val tile = ArrayTile.alloc(cellType, cols, rows)
    var i = 0
    val len = size
    while (i < len) {
      // 处理Double的时候使用updateDouble,可以使用Double.NaN代替Nodata
      tile.updateDouble(i, f(applyDouble(i)))
      i += 1
    }
    tile
}

综上,Abs算子执行流程如下:

1. 申请一个与原来尺寸相同的Tile
2. 定义两个Abs算子,一个面向Int型值,一个面向Double型值
3. 区分Tile内原始数据类型是否是浮点型数值:
   • 原始数据非浮点型,执行Int型Abs算子
   • 原始数据为浮点型,执行Double型Abs算子
4. 用算子结果填充新Tile并返回

2.2 挂载到Tile对象

LocalMethods特质中定义了一个localAbs函数,是对Abs算子的封装:

代码位于[geotrellis/raster/mapalgebra/local/LocalMethods.scala]

trait LocalMethods extends MethodExtensions[Tile]
                      with AddMethods { 
    def localAbs(): Tile =
        Abs(self)
}

LocalMethods又挂载到隐式类中,并通过Implicits特质暴露,实际通过隐式对象导出使用:

代码位于[geotrellis/raster/mapalgebra/local/Implicits.scala]

// 生成一个隐式转换对象
object Implicits extends Implicits

trait Implicits {
  implicit class withTileLocalMethods(val self: Tile) extends LocalMethods
}

通过引入import geotrellis.raster.mapalgebra.local._或import geotrellis.raster._,就能导入地图代数的全部Local类算子.此时Tile对象就可以执行localAbs算子了.

2.3 挂载到RDD[(K,Tile)]对象

挂载到Seq[(k,Tile)]对象的处理流程基本类似,在此不再展开.

在Spark环境中,我们经常会操作RDD[(K,Tile)]对象,当然我们也可以通过展开RDD为Tile的形式使用绑定在Tile上的函数,但是Geotrellis其实已经为我们做到了这一点.

由于算子核心都是一样的,因此只需引用基于Tile的Abs方法,并且将其应用于Rdd对象即可:

代码位于[geotrellis/spark/mapalgebra/local/LocalTileRDDMethods.scala]

import geotrellis.raster.mapalgebra.local._
// ...省略
trait LocalTileRDDMethods[K] extends TileRDDMethods[K] 
 // ...省略
{
    def localAbs() =
        self.mapValues { r => Abs(r) }
    // ...省略
}

其中,LocalTileRDDMethods特质继承了TileRDDMethods特质.定义如下:

代码位于[geotrellis/spark/mapalgebra/local/TileRDDMethods.scala]

// private[local] 仅在当前包可用,被其他包导入时不会污染环境
private[local] trait TileRDDMethods[K] extends MethodExtensions[RDD[(K, Tile)]] {
    // 定义一个上下文绑定的隐式对象,以便在擦除泛型k后指引编译器识别正确的类型
    implicit val keyClassTag: ClassTag[K]
}

同挂载到Tile对象上一样,挂载到RDD[(k,Tile)]也通过隐式对象暴露出来:

代码位于[geotrellis/spark/mapalgebra/local/Implicits.scala]

object Implicits extends Implicits

trait Implicits {
  implicit class withTileLocalMethods(val self: Tile) extends LocalMethods
  // ...省略
}

使用时引入import geotrellis.spark.mapalgebra.local._或import geotrellis.spark._即可.

至此,我们就大致明白Geotrellis是如何实现Abs算子并将算子挂载到常用对象了.我们可以将对Abs算子的分析推广到其他Local类算子.

3. 其他Local类算子

3.1 整体结构

简单来说,可以将全部Local类算子进行如下分类:

• 一元算子
  • 普通
    • Abs:绝对值
    • Ceil:向上取整
    • Round:向下取整
    • Log:取以10为底的对数(命名为Log10)和以自然对数为底的对数(命名为Log)
    • Sqrt:开平方
    • Negate:取负值
    • Not:按位取反
    • Defined:如果不是Nodata,则置为1,否则置0
    • Undefined:如果是Nodata,则置为1,否则置0
  • 三角函数相关
    • Sin/Cos/Tan
    • Acos/Asin/Atan/Atan2
  • 双曲函数相关
    • Sinh/Cosh/Tanh
• 二元算子
  • 二元计算
    • Add:A加B
    • Substract:A减B
    • Multiply:A乘B
    • Divide:A除以B,如果A或B有一为Nodata或B为0时,返回Nodata
    • Max:返回A和B中较大的那个,如果A或B有一为Nodata,返回Nodata
    • Min:返回A和B中较小的那个,如果A或B有一为Nodata,返回Nodata
    • And:A按位与B(两者都将被转换为整数),如果A或B有一为Nodata,返回Nodata
    • Or:A按位或B(两者都将被转换为整数),如果A或B有一为Nodata,返回Nodata
    • Xor:A按位异或B(两者都将被转换为整数),如果A或B有一为Nodata,返回Nodata
    • Pow:A的B次方
  • 二元比较
    • Equal:如果A与B相等,则返回True,否则为False
    • Unequal:如果A与B不相等,则返回True,否则为False
    • Greater:如果A比B大,则返回True,否则为False
    • GreaterOrEqual:如果A比B大或相等,则返回True,否则为False
    • Less:如果A比B小,则返回True,否则为False
    • LessOrEqual:如果A比B小或相等,则返回True,否则为False
• 复杂算子
  • Mask/InverseMask
  • IfCell
  • Majority/Minority
  • MaxN/MinN
  • Mean
  • Variety
  • Variance

3.2 二元计算算子-以Add算子为例

3.2.1 算子的实现

代码位于[geotrellis/raster/mapalgebra/local/Add.scala]

object Add extends LocalTileBinaryOp {
  def combine(z1: Int, z2: Int) =
    if (isNoData(z1) || isNoData(z2)) NODATA
    else z1 + z2

  def combine(z1: Double, z2: Double) =
    if (isNoData(z1) || isNoData(z2)) Double.NaN
    else z1 + z2
}

trait AddMethods extends MethodExtensions[Tile] {
  def localAdd(i: Int): Tile = Add(self, i)
  def +(i: Int): Tile = localAdd(i)
  def +:(i: Int): Tile = localAdd(i)
  def localAdd(d: Double): Tile = Add(self, d)
  def +(d: Double): Tile = localAdd(d)
  def +:(d: Double): Tile = localAdd(d)
  def localAdd(r: Tile): Tile = Add(self, r)
  def +(r: Tile): Tile = localAdd(r)
  def localAdd(rs: Traversable[Tile]): Tile = Add(self +: rs.toSeq)
  def +(rs: Traversable[Tile]): Tile = localAdd(rs)
}

相比简单的Abs算子,Add算子具有以下复杂性:

• 定义的Add单例对象扩展了LocalTileBinaryOp特质,其中定义了面向Double和Int类型的两个combine函数.
• 单独封装了一个名为AddMethods的特质,其中定义了若干操作符.

我们首先分析一下LocalTileBinaryOp特质.

3.2.2 LocalTileBinaryOp特质

Geotrellis中的二元操作需要满足什么样的设计?

• 需要区分Int和Double类型的数据(在Geotrellis中经常会遇到)
• 支持交换操作值,即:数值A [操作符] 数值B和数值B [操作符] 数值A都应该支持(当然可以拥有不同的意义,比如除运算)

我们来看一下Geotrellis是如何实现的:

代码位于[geotrellis/raster/mapalgebra/local/LocalTileBinaryOp.scala]

trait LocalTileBinaryOp extends Serializable {
  val name = {
    // 获取实际操作符名称
    val n = getClass.getSimpleName
    if(n.endsWith("$")) n.substring(0, n.length - 1)
    else n
  }

  // 兼容Int和Double型输入值
  def apply(r: Tile, c: Int): Tile =
    r.dualMap(combine(_, c))({
      val d = i2d(c)
      (z: Double) => combine(z, d)
    })
  def apply(r: Tile, c: Double): Tile =
    r.dualMap({z => d2i(combine(i2d(z), c))})(combine(_, c))

  // 兼容交换数值顺序
  def apply(c: Int, r: Tile): Tile =
    r.dualMap(combine(c, _))({
      val d = i2d(c)
      (z: Double) => combine(d, z)
    })
  def apply(c: Double, r: Tile): Tile =
    r.dualMap({z => d2i(combine(c, i2d(z)))})(combine(c, _))

  // 支持操作符都是Tile的情况
  def apply(r1: Tile, r2: Tile): Tile = 
    r1.dualCombine(r2)(combine)(combine)
  def apply(rs: Traversable[Tile]): Tile = 
    new TileReducer(combine)(combine)(rs)
    
  // 预定义Int型/Double型算子
  def combine(z1: Int, z2: Int): Int
  def combine(z1: Double, z2: Double): Double
}

综上,一个二元计算算子执行流程如下:

1. 申请一个与原来尺寸相同的Tile
2. 区分Tile内原始数据类型是否是浮点型数值
3. 区分传入值是否为Double类型
4. 最终可能出现4种情况:
   1. 原始数据非浮点型,传入值为Int型,则直接执行Int型combine算子
   2. 原始数据非浮点型,传入值为Double型,则将传入值转换为Int后再执行Double型combine算子(在算子中原始的非浮点数据将转换为Double型以配合算子数据类型定义)
   3. 原始数据为浮点型,传入值为Int型,则将传入值转换为Double后再执行Int型combine算子(在算子中原始的浮点数据将转换为Int型以配合算子数据类型定义)
   4. 原始数据为浮点型,传入值为Double型,则直接执行Double型combine算子
5. 用算子结果填充新Tile并返回

因此我们只需如同Add单例那样定义Int/Double类型的combine函数并扩展LocalTileBinaryOp特质,即可得到对应的二元计算算子.

3.2.2.1 传入值为Tile类型的情况

除了常见的传入值为Int/Double的情况,LocalTileBinaryOp还实现了传入值是Tile的情况,即对于每一个位置,其传入值不是固定的,而是来自于另一个Tile的相同位置:

def apply(r1: Tile, r2: Tile): Tile = 
    r1.dualCombine(r2)(combine)(combine)

核心为Tile.dualCombine方法:

代码位于[geotrellis/raster/Tile.scala]

  def dualCombine(r2: Tile)(f: (Int, Int) => Int)(g: (Double, Double) => Double): Tile =
    if (cellType.union(r2.cellType).isFloatingPoint) combineDouble(r2)(g) else combine(r2)(f)

dualCombine又是对combine/combineDouble的封装.他们的定义位于具体的Tile类中,我们依然以较为简单的ArrayTile为例:

代码位于[geotrellis/raster/ArrayTile.scala]

  def combine(other: ArrayTile)(f: (Int, Int) => Int): ArrayTile = {
    // 确保当前Tile与传入值Tile具有同样的尺寸
    (this, other).assertEqualDimensions()

    val output = ArrayTile.alloc(cellType.union(other.cellType), cols, rows)
    var i = 0
    val len = size
    while (i < len) {
      // 本Tile在i位置的值和传入Tile的i位置的值计算结果作为输出值
      output(i) = f(apply(i), other(i))
      i += 1
    }
    output
  }

  // ...省略处理Double类型的类似逻辑

3.2.2.2 传入值为Traversable[Tile]类型的情况

传入值为Traversable[Tile]类型传入值为Tile类型的复杂化:

def apply(rs: Traversable[Tile]): Tile = 
    // 第一个combine对应Int型,第二个对应Double型
    new TileReducer(combine)(combine)(rs)

核心为TileReducer对象:

代码位于[geotrellis/raster/mapalgebra/local/TileReducer.scala]

class TileReducer(handle: (Int, Int)=>Int)(handleDouble: (Double, Double)=>Double) {
  def apply(seq: Traversable[Tile]): Tile =
    apply(seq.toList)

  def apply(list: List[Tile]): Tile =
    handleTiles(list)
  
  // 尾递归优化注释
  @tailrec final def reduce(d: Tile, rasters: List[Tile]): Tile = {
    // 检测未处理的Tile的状态
    rasters match {
      // 没有未处理Tile,不再计算,结束递归
      case Nil => d
      // 还有起码一个未处理的Tile
      case r :: rs => if (r.cellType.isFloatingPoint) {
        reduceDouble(d.combineDouble(r)(handleDouble), rs)
      } else {
        // 最终转换为传入值是单个Tile的情况
        reduce(d.combine(r)(handle), rs)
      }
    }
  }

  // ...省略处理Double类型的类似逻辑

  def handleTiles(rasters: List[Tile]) = {
    // 取列表中的第一个元素
    val (r :: rs) = rasters

    if (r.cellType.isFloatingPoint) {
      reduceDouble(r, rs)
    } else {
      reduce(r, rs)
    }
  }
}

其原理为用递归的方式将针对Traversable[Tile]对象的操作化简为多次针对Tile的操作.

3.2.3 AddMethods特质

AddMethods特质封装了许多操作符,方便我们的实际使用,因为相比一元算子,二元算子兼容的传入值更为复杂:

trait AddMethods extends MethodExtensions[Tile] {
  def localAdd(i: Int): Tile = Add(self, i)
  def +(i: Int): Tile = localAdd(i)
  def +:(i: Int): Tile = localAdd(i)
  def localAdd(d: Double): Tile = Add(self, d)
  def +(d: Double): Tile = localAdd(d)
  def +:(d: Double): Tile = localAdd(d)
  def localAdd(r: Tile): Tile = Add(self, r)
  def +(r: Tile): Tile = localAdd(r)
  def localAdd(rs: Traversable[Tile]): Tile = Add(self +: rs.toSeq)
  def +(rs: Traversable[Tile]): Tile = localAdd(rs)
}

同时也扩展了MethodExtensions特质:

代码位于[geotrellis/util/MethodExtensions.scala]

// 协变类型T
trait MethodExtensions[+T] extends Serializable {
  def self: T
}

定义了名为self的变量,且self的实际类型可以为Tile的各种子类,如ArrayTile/ArrayMultiBandTile等.

最终,AddMethods通过with的方式挂载到了上文提到的LocalMethods特质中,可以通过形如tile_1 + tile_2的形式使用localAdd算子.

同时,因为Add算子支持传入值为Traversable[Tile]的类型,因此在LocalSeqMethods特质中也注册了一个操作符:

代码位于[geotrellis/raster/mapalgebra/local/LocalSeqMethods.scala]

trait LocalSeqMethods extends MethodExtensions[Traversable[Tile]] {
  def localAdd(): Tile =
    Add(self)
  def +(): Tile = localAdd()
}

这样也可以通过形如+List(tile_1, tile_2, tile_3)的形式在Traversable[Tile]对象上使用localAdd算子.

3.3 二元比较算子-以Equal算子为例

3.3.1 算子的实现

代码位于[geotrellis/raster/mapalgebra/local/Equal.scala]

object Equal extends LocalTileComparatorOp {
  def compare(z1:Int,z2:Int):Boolean =
    if(z1 == z2) true else false

  def compare(z1:Double,z2:Double):Boolean =
    if(isNoData(z1)) { if(isNoData(z2)) true else false }
    else {
      if(isNoData(z2)) { false }
      else {
        if(z1 == z2) true
        else false
      }
    }
}

trait EqualMethods extends MethodExtensions[Tile] {
  def localEqual(i: Int): Tile = Equal(self, i)
  def localEqual(d: Double): Tile = Equal(self, d)
  def localEqual(r: Tile): Tile = Equal(self, r)
}

Equal算子与Add算子的显著不同在于扩展了LocalTileComparatorOp特质,与之相对应,定义了两个返回布尔值的compare函数.

3.3.2 LocalTileComparatorOp特质

与LocalTileBinaryOp特质相同,LocalTileComparatorOp特质也是二元操作,不同的是,LocalTileComparatorOp特质更加简单,因为二元比较算子最终永远返回一个Bit类型的二值Tile.

核心代码如下:

代码位于[geotrellis/raster/mapalgebra/local/LocalTileComparatorOp.scala]

trait LocalTileComparatorOp extends Serializable {

  def apply(r: Tile, c: Int): Tile = {
    // 申请一个bit类型的Tile承载结果
    val tile = BitArrayTile.ofDim(r.cols, r.rows)
    if (r.cellType.isFloatingPoint) {
      val cons = c.toDouble
      // 直接遍历全部像元
      cfor(0)(_ < r.rows, _ + 1) { row =>
        cfor(0)(_ < r.cols, _ + 1) { col =>
          // 根据compare结果直接赋值
          tile.set(col, row, if (compare(r.getDouble(col, row), cons)) 1 else 0)
        }
      }
    } else {
      // ...省略处理int类型的近似逻辑
    }
    tile
  }

  // ...省略

  def apply(r1: Tile, r2: Tile): Tile = {
    // 验证两个Tile是否具有相同尺寸
    Traversable(r1, r2).assertEqualDimensions()
    val Dimensions(cols, rows) = r1.dimensions
    val tile = BitArrayTile.ofDim(cols, rows)

    cfor(0)(_ < r1.rows, _ + 1) { row =>
      cfor(0)(_ < r1.cols, _ + 1) { col =>
        if (r1.cellType.isFloatingPoint || r2.cellType.isFloatingPoint) {
          // 用两个Tile相同位置的数值做对比,直接赋值
          tile.set(col, row, if (compare(r1.getDouble(col, row), r2.getDouble(col, row))) 1 else 0)
        } else {
          tile.set(col, row, if (compare(r1.get(col, row), r2.get(col, row))) 1 else 0)
        }
      }
    }
    tile
  }

  def compare(z1: Int, z2: Int): Boolean
  def compare(z1: Double, z2: Double): Boolean
}

从代码中可以看出,二元比较算子是简化版的二元计算算子,其实可以用二元计算算子模拟二元比较算子,但直接使用二元比较算子在内存使用和计算速度上都更有优势.

3.4 Mask(掩膜)/InverseMask(反掩膜)

Mask/InverseMask互为相反算子,我们以Mask为例.

MasK算子实现了这样一种功能:准备原始Tile和掩膜Tile,掩膜值和结果值4个输入参数,对于原始Tile中的某个位置,如果其对应的掩膜Tile的值为掩膜值,则设为结果值,否则保持原值:

代码位于[geotrellis/raster/mapalgebra/local/Mask.scala]

object Mask extends Serializable {
  // 掩膜Tile应为Int类型,比如二值图
  def apply(r1: Tile, r2: Tile, readMask: Int, writeMask: Int): Tile =
    r1.dualCombine(r2) { (z1,z2) => if (z2 == readMask) writeMask else z1 }
                       { (z1,z2) => if (d2i(z2) == readMask) i2d(writeMask) else z1 }
}

代码本身比较简单,但值得注意的是,它并没有直接挂载到Tile对象上,而是只挂载到Seq[(K,Tile)]和RDD[(K,Tile)]对象上:

3.4.1 挂载到Seq[(K,Tile)]对象

Mask算子通过注册到LocalTileCollectionMethods特质上,最终实现挂载到Seq[(K,Tile)]对象上:

代码位于[geotrellis/raster/mapalgebra/local/Mask.scala]

trait LocalTileCollectionMethods[K] extends MethodExtensions[Seq[(K, Tile)]]
{
    def localMask(other: Seq[(K, Tile)], readMask: Int, writeMask: Int): Seq[(K, Tile)] =
        self.combineValues(other) {
          case (r1, r2) => Mask(r1, r2, readMask, writeMask)
    }
    
    // ...省略
}

代码的核心为combineValues方法:

代码位于[geotrellis/layer/mapalgebra/CollectionMethods.scala]

abstract class CollectionCombineMethods[K, V] extends MethodExtensions[Seq[(K, V)]] {
  def combineValues[R](other: Seq[(K, V)])(f: (V, V) => R): Seq[(K, R)] =
    (self ++ other).groupBy(_._1).mapValues { case Seq((_, v1), (_, v2)) => f(v1, v2) }.toSeq

  // ...省略
}

根据combineValues的代码,可以归纳出Mask算子的流程:

1. 将原始Tile与掩膜Tile融合到一个序列中
2. 根据其索引(可以是空间索引,也可以是时空索引),将它们归为两两一组
3. 每个匹配上索引的原始Tile/掩膜Tile对进行Mask算子操作,最终返回结果

值得注意的是:掩膜Tile必须与原始Tile具有一致的索引,即要求全部掩膜和全部原始Tile都要完全对应,不完全对应会报错.
我估计这也是为何只在Seq[(k,Tile)]对象和RDD[(k,Tile)]对象中挂载了Mask算子,因为具有索引可以从逻辑上保证匹配,也在逻辑上更合理.

3.4.2 挂载到RDD[(K,Tile)]对象

逻辑基本与挂载到Seq[(k,Tile)]一致:

trait LocalTileRDDMethods[K] extends TileRDDMethods[K] {
    def localMask(other: RDD[(K, Tile)], readMask: Int, writeMask: Int, partitioner: Option[Partitioner] = None): RDD[(K, Tile)] =
        self.combineValues(other, partitioner) {
            case (r1, r2) => Mask(r1, r2, readMask, writeMask)
    }
    
    // ... 省略
}

代码的核心为combineValues方法:

代码位于[geotrellis/spark/mapalgebra/CollectionMethods.scala]

abstract class CombineMethods[K: ClassTag, V: ClassTag] extends MethodExtensions[RDD[(K, V)]] {
  def combineValues[R: ClassTag](other: RDD[(K, V)])(f: (V, V) => R): RDD[(K, R)] = combineValues(other, None)(f)
  def combineValues[R: ClassTag](other: RDD[(K, V)], partitioner: Option[Partitioner])(f: (V, V) => R): RDD[(K, R)] =
    partitioner
      .fold(self.join(other))(self.join(other, _))
      .mapValues { case (tile1, tile2) => f(tile1, tile2) }

  // ... 省略
}

不同的地方为使用了RDD.join实现了按索引进行归类.

3.5 IfCell(条件)

IfCell算子实现了这样一种功能:将Tile中全部符合给定条件(一个输入像素值返回布尔值的函数)的像素设为给定值.

具体又分为以下几种:

• 将Tile中全部符合给定条件的像素值设为给定值,否则保留原始值
• 将Tile中全部符合给定条件的像素值设为给定值A,不符合的设为给定值B
• 可以提供一个额外Tile,与原始Tile一起作为作为条件判断参数

IfCell的实现基本可以看做是将稍微复杂一些的二元比较算子,在此不再展开.

3.6 Majority(多数值)/Minority(少数值)

Majority/Minority互为相反算子.我们以Majority为例.

Majority算子实现了这样一种功能:给定一系列尺寸相同的Tile,针对Tile中的每一个位置,取这一系列Tile中该位置出现最多(或第N多)的值,最终输出一个Tile:

代码位于[geotrellis/spark/mapalgebra/Majority.scala]

object Majority extends Serializable {
  
  // 支持一个或多个Tile作为参数
  def apply(r: Tile*): Tile = apply(0, r)

  def apply(level: Int, r: Tile*): Tile = apply(level, r)
  
  // 默认为计算出现最多次
  def apply(rs: Traversable[Tile])(implicit d: DI): Tile = apply(0, rs)
  
  // DI为DummyImplicit
  def apply(level: Int, rs: Traversable[Tile])(implicit d: DI): Tile = {
    rs.assertEqualDimensions

    val layerCount = rs.toSeq.length
    if(layerCount == 0) {
      sys.error("Can't compute majority of empty sequence")
    } else {
      // 计算合适的结果类型,避免溢出
      val newCellType = rs.map(_.cellType).reduce(_.union(_))
      // 以第一个Tile的尺寸为准
      val Dimensions(cols, rows) = rs.head.dimensions
      val tile = ArrayTile.alloc(newCellType, cols, rows)

      if(newCellType.isFloatingPoint) {
        val counts = mutable.Map[Double, Int]()

        cfor(0)(_ < rows, _ + 1) { row =>
          cfor(0)(_ < cols, _ + 1) { col =>
            counts.clear
            // 收集该位置全部Tile的像素值
            for(r <- rs) {
              val v = r.getDouble(col, row)
              // 排除Nodata值的情况
              if(isData(v)) {
                if(!counts.contains(v)) {
                  counts(v) = 1
                } else {
                  counts(v) += 1
                }
              }
            }
            // 按出现次数排序
            val sorted =
              counts.keys
                .toSeq
                .sortBy { k => counts(k) }
                .toList
            val len = sorted.length - 1
            // 只有在数值足够多的时候赋值,否则设为NaN
            val m =
              if(len >= level) { sorted(len-level) }
              else { Double.NaN }
            tile.setDouble(col, row, m)
          }
        }
      } else {
       // ... 省略处理Int类型的类似逻辑
      }
      tile
    }
  }
}

从代码中可以看出,Majority算子同时支持Tile和Traversable[Tile]作为参数,因此与前文的Add算子类似,同时挂载到LocalMethods和LocalSeqMethods特质中.

3.7 MaxN(第N大)/MinN(第N小)

MaxN/MinN互为相反算子.我们以MaxN为例.

Majority算子实现了这样一种功能:给定一系列尺寸相同的Tile,针对Tile中的每一个位置,取这一系列Tile中该位置第N大的值,最终输出一个Tile.

算法使用基于快速排序的随机就地快速选择算法:

代码位于[geotrellis/spark/mapalgebra/local/MaxN.scala]

object MaxN extends Serializable {
  // 注意ArrayView被定义为case class
  case class ArrayView[Num](arr: Array[Num], from: Int, until: Int) {
    // 从ArrayView中取值
    def apply(n: Int) =
      if (from + n < until) arr(from + n)
      else throw new ArrayIndexOutOfBoundsException(n)

    // 根据基准值就地分区
    def partitionInPlace(p: Num => Boolean): (ArrayView[Num], ArrayView[Num]) = {
      // 初始化大小值索引
      var upper = until - 1
      var lower = from
      // 最终将序列排序为元素全大于基准值和元素全不大于基准值的左右两个部分
      while (lower < upper) {
        while (lower < until && p(arr(lower))) lower += 1
        while (upper >= from && !p(arr(upper))) upper -= 1
        if (lower < upper) { val tmp = arr(lower); arr(lower) = arr(upper); arr(upper) = tmp }
      }
      // 利用case class的copy方法将排序好的数组拆为两个独立的ArrayView
      (copy(until = lower), copy(from = lower))
    }

    def size = until - from
    def isEmpty = size <= 0
  }
  object ArrayView {
    def apply(arr: Array[Double]) = new ArrayView(arr, 0, arr.size)
    def apply(arr: Array[Int]) = new ArrayView(arr, 0, arr.size)
  }

  def findNthDoubleInPlace(arr: ArrayView[Double], n: Int): Double = {
    if(n >= arr.size) {
      Double.NaN
    } else {
      // 从输入值中随机选择一个作为中点
      val pivot = arr(scala.util.Random.nextInt(arr.size))
      // 分为左右部分,左列表元素全大于基准值,右列表元素全不大于基准值
      val (left, right) = arr partitionInPlace (_ > pivot)
      // 
      if (left.size == n) pivot
      // 左列表为空,说明全部元素都不大于基准值
      else if (left.isEmpty) {
        // 重新计算,看看等于基准值的有多少 
        val (left, right) = arr partitionInPlace (_ == pivot)
        // 等于基准值的数值够多,则说明基准值就是要找的值
        if (left.size > n) pivot
        // 基准值不够,则从剩下的部分中寻找
        else findNthDoubleInPlace(right, n - left.size)
      } else if (left.size < n) findNthDoubleInPlace(right, n - left.size) // 从右列表中寻找
      else findNthDoubleInPlace(left, n) // 左列表不为空,且不满足上面的情况,说明左列表太大.需要进一步缩小范围
    }
  }
  
  // ...省略处理int类型的近似逻辑

  def apply(n: Int, rs: Tile*): Tile =
    apply(n, rs)

  def apply(n: Int, rs: Traversable[Tile])(implicit d: DI): Tile = {
    rs.assertEqualDimensions

    val layerCount = rs.toSeq.length
    if(layerCount < n) {
      sys.error(s"Not enough values to compute Nth")
    } else {
      val newCellType = rs.map(_.cellType).reduce(_.union(_))
      val Dimensions(cols, rows) = rs.head.dimensions
      val tile = ArrayTile.alloc(newCellType, cols, rows)

      cfor(0)(_ < rows, _ + 1) { row =>
        cfor(0)(_ < cols, _ + 1) { col =>
          if(newCellType.isFloatingPoint) {
            // 计算并赋值
            val maxN = findNthDoubleInPlace(ArrayView(rs.map(r => r.getDouble(col, row)).filter(num => !isNoData(num)).toArray), n)
            tile.setDouble(col, row, maxN)
          }else { // integer values
            // ...省略处理int类型的近似逻辑
          }
        }
      }
      tile
    }
  }
}

值得注意的是,MaxN算子的输入参数仅支持Traversable[Tile],因此,也只挂载到LocalSeqMethods特质中.

3.8 Mean(均值)/Variance(方差)/Variety(出现不同数字的数量)

Mean/Variance/Variety算子实现的功能顾名思义(具体实现也很简单,在此不再展开),他们也都具有相同的执行流程:

1. 仅支持Traversable[Tile]作为参数
2. 以第一个Tile作为标准尺寸,创建一个空Tile
3. 遍历空Tile的每一个像元
4. 针对每一个像元位置,聚合所有Tile中的像元值,执行Mean/Variance/Variety算子,将结果赋予空Tile的对应位置
5. 最终输出赋值后的Tile

4. Temporal类算子

因为可以包含带有时间属性的索引Key,因此Seq[(K,Tile)]和RDD[(k,Tile)]的Local类算子中额外添加了Temporal类算子,针对时间值实现了Mean/Max/Min/Variance 4种算子.

在RDD[(k,Tile)]对象上的实现逻辑类似,在此不再展开.

我们以Seq[(K,Tile)]对象上的Temporal类算子为例:

代码位于[geotrellis/layer/mapalgebra/local/temporal/LocalTemporalStatistics.scala]

object LocalTemporalStatistics {
  import TemporalWindowHelper._

  // 核心代码
  // 确保key为时空索引
  private[geotrellis] def aggregateWithTemporalWindow[K: SpatialComponent: TemporalComponent](
    sourceSeq: Seq[(K, Tile)],
    windowSize: Int,
    unit: Int,
    start: ZonedDateTime,
    end: ZonedDateTime)(reduceOp: Traversable[Tile] => Tile): Seq[(K, Tile)] = {
    // 1.筛选时间在给定范围内的Tile
    val seq =
      sourceSeq
        .map { case (key, tile) =>
          val SpatialKey(col, row) = key.getComponent[SpatialKey]
          val time = key.getComponent[TemporalKey].time
          val startDiff = getDifferenceByUnit(unit, start, time)
          val endDiff = getDifferenceByUnit(unit, time, end)

          val newKey =
            if (startDiff < 0 || endDiff < 0) {
              (-1, col, row)
            }
            else {
              // 将时间归一到同一个时间窗口中去,避免数值太分散,加大计算量
              val timeDelimiter = startDiff / windowSize
              (timeDelimiter, col, row)
            }

          (newKey, (key, tile))
        }
        .filter { case ((i, col, row), _) => i >= 0 }
    
    seq.groupBy(_._1) // 2.按时间窗口分组
      .map { case (k, v) => (k, v.map(_._2)) } 
      // 3. 拆分为Seq[(newKey,Seq[(k,Tile)])] 
      .map { case (_, iter) =>
        val (keys, tiles) = iter.unzip
        // 不再用窗口值,而是以实际时间值中的最小值作为key
        val key = keys.min(Ordering.by { key: K => key.getComponent[TemporalKey].time })
        // 4. 执行指定的算子操作
        val tile = reduceOp(tiles)
        // 5. 返回结果
        (key, tile)
      }.toSeq
  }

  // 计算时间差值,如果time早于base,则返回负值
  private[geotrellis] def getDifferenceByUnit(unit: Int, base: ZonedDateTime, time: ZonedDateTime): Int =
    Math.toIntExact(unit match {
      case UnitSeconds => SECONDS.between(base, time)
      case UnitMinutes => MINUTES.between(base, time)
      // ... 省略逻辑类似的其他事件单位
      case _ => throw new IllegalStateException(s"Bad unit $unit.")
    })

  private[geotrellis] def minReduceOp(tiles: Traversable[Tile]): Tile = tiles.localMin

  // ... 省略逻辑类似的Max/Mean/Variance
  
  // 最终封装为temporalMin方法
  def temporalMin[K: SpatialComponent: TemporalComponent](
    seq: Seq[(K, Tile)],
    windowSize: Int,
    unit: Int,
    start: ZonedDateTime,
    end: ZonedDateTime): Seq[(K, Tile)] =
    aggregateWithTemporalWindow(seq, windowSize, unit, start, end)(minReduceOp)

  // ... 省略逻辑类似的Max/Mean/Variance

}

5. 总结

在本篇中我们探讨了Geotrellis地图代数中Local类算子,研究了Local类算子的主要结构以及部分算子的实现方法和挂载到Tile/Seq[(k,Tile)]/RDD[(k,Tile)]对象的逻辑,在下一篇我们将继续研究Focal类算子的结构与实现.