1. 引入

我们会使用Geotrellis自然是因为它能利用spark进行高效的分布式计算,而分布式计算无论在数据读取还是计算上都与常规计算有许多不同.

我们就以计算一张影像的NDVI结果为例,看看Geotrellis是如何在spark上读取计算的:

package demo

import geotrellis.layer._
import geotrellis.util._
import geotrellis.spark._
import geotrellis.spark.store.hadoop.HadoopGeoTiffRDD
import geotrellis.raster._
import geotrellis.raster.io.geotiff._
import geotrellis.raster.ResampleMethods.NearestNeighbor
import org.apache.spark._
import org.apache.hadoop.fs.Path

object MyApp extends App {

    val sparkConf: SparkConf = new SparkConf().setMaster("local[4]").setAppName("ndvi")
    implicit val sc: SparkContext = new SparkContext(sparkConf)
    // 1. 数据读取
    val geoTiffRDD = HadoopGeoTiffRDD.spatialMultiband(new Path("/Users/temp/clip.tif"))
    // 2. 提取数据集的元数据并转换索引
    val (_, metadata) = geoTiffRDD.collectMetadata[SpatialKey](FloatingLayoutScheme())
    val tiledRDD = geoTiffRDD.tileToLayout[SpatialKey](metadata, NearestNeighbor)
    // 3. 计算NDVI
    val raster = tiledRDD.withContext { rdd =>
      rdd.mapValues { tile =>
        tile.convert(DoubleConstantNoDataCellType).combineDouble(0, 1) { (r, nir) =>
          (nir - r) / (nir + r)
        }
      }.reduceByKey(_.localMax(_))
    }.stitch
    // 4. 写入结果
    GeoTiff(raster, metadata.crs).write("/Users/temp/ndvi.tif")
}

代码分为四个步骤:

1. 数据读取
2. 提取元数据
3. 计算NDVI
4. 写入结果

按顺序我们先从数据读取开始.

这部分的主体是HadoopGeoTiffRDD.spatialMultiband方法,整理后代码如下:

[代码位于HadoopGeoTiffRDD.scala]

  def apply[ProjectedExtent, ProjectedExtent, MultibandTile](
    path: Path = new Path("/Users/temp/clip.tif"),
    uriToKey: (URI, I) => K = (_: URI, key: MultibandTile) => key,
    options: Options = Options.DEFAULT,
  )(implicit sc: SparkContext, rr: RasterReader[Options, (ProjectedExtent, MultibandTile)]): RDD[(ProjectedExtent, MultibandTile)] = {

    val conf = new SerializableConfiguration(configuration(path, options))
    val pathString = path.toString 
    
    // 初始化读取器(没有实际动作)
    val infoReader = HadoopGeoTiffInfoReader(pathString, conf, options.tiffExtensions)
    
    // 1. 读取影像数据地址
    val files: RDD[URI] =infoReader.geoTiffInfoRDD.map(new URI(_))
    
    // 2. 实际读取影像数据
    infoReader.readWindows(
      files,
      uriToKey,
      maxTileSize,
      options.partitionBytes.getOrElse(DefaultPartitionBytes),
      options,
      geometry)
  }

代码分为两个步骤:

1. 读取影像数据地址:从HDFS中读取目录中全部影像文件的URI
2. 实际读取影像数据:调用readWindows方法,将URI读取为RDD[(ProjectedExtent, MultibandTile)]

按顺序我们先从读取影像数据地址开始.

2. 读取影像数据地址

读取数据是计算的起点,而分布式计算往往需要分布式的存储介质.在本例中存储介质是默认的HDFS.

这部分的主体是infoReader.geoTiffInfoRDD方法,整理后代码如下:

[代码位于HadoopGeoTiffInfoReader.scala]

def geoTiffInfoRDD(implicit sc: SparkContext): RDD[String] = {
    // 4. 将数组转换为RDD数组
    sc.parallelize(
        // 1. 列出该目录下全部文件
        HdfsUtils.listFiles(new Path(path), config.value)
        // 2. 转换为字符串
        .map({ path => path.toString }) 
        // 3. 筛选影像文件
        .filter({ path => tiffExtensions.exists({ e => path.endsWith(e) }) })
    )
}

代码分为四个步骤:

1. 利用HDFS提供的接口,获取某个目录下下全部文件.
2. 因为HDFS接口返回的是Path类型对象,因而需要提取里面的字符串,方便后面使用.
3. 根据后缀名筛选影像文件.
4. 将数组并行化,转换为RDD数组.

我们可以提出两个问题:

1. Geotrellis规划数据的思路为何?
   • 从代码中不难看出,Geotrellis从一开始就关注于数据集,在HDFS中,数据集就是特定的文件夹(尽管单个文件路径也支持),这样处理大量数据时,我们就无需将他们的文件名一一指明.
2. 为何要先读取影像的地址,而不是直接读取数据?
   • 在这一阶段Geotrellis并没有进行分布式读取,获取需要处理的全部影像的实际路径是在单节点完成的低成本操作,无需也不应该并行化.
   • 先整理好要处理的影像,将其转换为RDD数组,有了RDD对象,Geotrellis就能使用Spark进行操作了.在后面我们也能看到这种思路.即先规划任务,再实际执行计算.

准备好待处理的影像地址,就可以实际读取影像数据了.

3. 实际读取影像数据

实际读取影像就是我们要了解的核心功能了.我们要尝试回答一个最根本的问题:Geotrellis如何实现分布式的数据读取的?

这部分的主体是infoReader.readWindows方法,整理后代码如下:

[代码位于GeoTiffInfoReader.scala]

def readWindows[O, I, K, V](
    files: RDD[URI],
    uriToKey: (URI, I) => K,
    maxSize: Int = 256,
    partitionBytes: Long = 128 * 1024 * 1024,
    rasterReaderOptions: O,
    geometry: Option[Geometry]
  )(
    implicit sc: SparkContext, rr: RasterReader[O, (I, V)]
  ): RDD[(K, V)] = {
    
    // 2. 规划读取窗口
    val windows: RDD[(URI, Array[GridBounds[Int]])] =
      files.flatMap({ uri =>
        windowsByPartition(
          info = getGeoTiffInfo(uri.toString),// 1. 读取影像信息
          maxSize = maxSize,
          partitionBytes = partitionBytes,
          geometry = geometry
        ).map { windows => (uri, windows) }
      })
    
    // 为了减少重复读取操作,进行持久化
    windows.persist()
    
    // 如有必要则重分区
    val repartition = {
      // 代码中的第一个action算子,实际计算从此开始
      val windowCount = windows.count.toInt
      // 如果规划的分区数大于实际的分区数,则按照我们规划的分区数来重分区
      if (windowCount > windows.partitions.length) {
        windows.repartition(windowCount)
      }
      else windows
    }
    
    // 3. 开始读取数据
    val result = repartition.flatMap { case (uri, windows) =>
      val info = getGeoTiffInfo(uri.toString)
      rr.readWindows(windows, info, rasterReaderOptions).map { case (k, v) =>
        uriToKey(uri, k) -> v
      }
    }
    
    // 取消持久化
    windows.unpersist()
    result
}

代码分为三个主要步骤:

1. 读取影像信息.
2. 规划读取窗口.
3. 开始读取数据.

按照顺序,我们从读取影像信息开始.

3.1 读取影像信息

这部分的主体是getGeoTiffInfo方法,整理后代码如下:

[代码位于HadoopGeoTiffInfoReader.scala]

def getGeoTiffInfo(uri: String): GeoTiffInfo = {
    // 1. 创建HDFS块读取器,读取指定数据
    val rr = HdfsRangeReader(new Path(uri), config.value)
    
    // 准备ovr参数
    val ovrPath = new Path(s"${uri}.ovr")
    val ovrReader: Option[ByteReader] =
      if (HdfsUtils.pathExists(ovrPath, config.value)) Some(HdfsRangeReader(ovrPath, config.value)) else None
    
    // 2. 读取影像信息
    // HDFSRangeReader通过隐式转换,转换为所需的StreamingByteReader类型
    GeoTiffInfo.read(rr, streaming, true, ovrReader)
}

代码分为两个步骤:

1. 创建HDFS块读取器(RangeReader).
   • 在分布式的存储机制中,一般都会采取块读取的方式,因为相比本地读取,有固定的时间成本花在了建立网络连接上下文上,因此应尽量减少读取频次,提高单词读取量,即以"块"为读取粒度.采用合适大小的块可以取得资源占用与效率的平衡.
2. 调用GeoTiffInfo.read读取元数据.
   • 代码中存在将RangeReader(块读取)转换为ByteReader(Byte读取)的逻辑.其实现比较复杂,在此不展开,但核心思想较为简单:把块读取包装为具有缓冲区的流读取,底层读取一个块,上层读取该块中的字节,读完了底层就再读一个块,上层接着读取接下来的字节,这样对于上层的使用者来说其行为与标准Byte读取无异,只是读取速度时快时慢而已.
   • 值得一提的是,最终返回的GeoTiffInfo对象,其默认为streaming模式,即按需读取模式,这个配置使得后面读取影像实际值时效率更高.

准备好影像的元数据信息,就可以调用规划读取窗口的方法了.

3.2 规划读取窗口

实际读取的过程依然符合前述的先规划任务,再实际执行计算思路,因此要先规划读取任务.

这部分的主体是windowsByPartition方法,整理后代码如下:

[代码位于GeoTiffInfoReader.scala]

  private def windowsByPartition(
    info: GeoTiffInfo,
    maxSize: Int,
    partitionBytes: Long,
    geometry: Option[Geometry]
  ): Array[Array[GridBounds[Int]]] = {
    
    // 1. 列出全部的window
    val: Array[GridBounds[Int]] windows = info.segmentLayout.listWindows(maxSize)
    
    // 2. 对window分组
    info.segmentLayout.partitionWindowsBySegments(windows, partitionBytes / math.max(info.cellType.bytes, 1))
  }

代码分为两个步骤:

1. 根据影像的元数据,列出全部的Window.
2. 将这些Window对象进行分组,形成若干Window集合

按照顺序,我们先从列出全部Window开始

3.2.1 列出全部Window

这部分的主体是info.segmentLayout.listWindows方法,整理后代码如下:

[代码位于GeoTiffInfoReader.scala]

def listWindows(maxSize: Int): Array[GridBounds[Int]] = {
    // 回顾读取元数据的部分:
    // tileCols/Rows代表一个Segment的一列/行有多少个像素
    val segCols = tileLayout.tileCols
    val segRows = tileLayout.tileRows
    
    //1. 计算一个window的大小
    val colSize: Int =
      if (maxSize >= segCols * 2) { 
        // 如果实际的Segment的tilescol小于等于128,colsize取不大于256的tilecols的最大倍数.
        math.floor(maxSize.toDouble / segCols).toInt * segCols
      } else if (maxSize >= segCols) {
        // tilecols在128-256之间则取tilecols实际值
        segCols
      } else //tilecols大于256,则取tilecols不大于256的最大约数
      bestWindowSize(maxSize, segCols)

    val rowSize: Int =
      if (maxSize >= segRows * 2) {
        math.floor(maxSize.toDouble / segRows).toInt * segRows
      } else if (maxSize >= segRows) {
        segRows
      } else bestWindowSize(maxSize, segRows)
    
    // 2. 根据计算的大小生成window集合
    val windows = listWindows(colSize, rowSize)

    windows
}

def listWindows(cols: Int, rows: Int): Array[GridBounds[Int]] = {
    // 给全部window预留位置
    val result = scala.collection.mutable.ArrayBuilder.make[GridBounds[Int]]
    result.sizeHint((totalCols / cols) * (totalRows / rows))
    // 双层遍历生成全部的GridBounds
    // 这种先从上到下,再从左到右的window排布方式并不符合一般数据在Tiff文件中的实际排布
    // 但因为Tiff文件并没有限定通过何种排列方式,因此不能假定其实际顺序
    cfor(0)(_ < totalCols, _ + cols) { col =>
      cfor(0)(_ < totalRows, _ + rows) { row =>
        result +=
          GridBounds(
            col,
            row,
            math.min(col + cols - 1, totalCols - 1),
            math.min(row + rows - 1, totalRows - 1)
          )
      }
    }
    result.result
}

代码分为了两个主要部分:

1. 根据元数据中的Segment信息计算一个Window的合适大小
   • 所谓Window,即一个GridBounds对象.GridBounds可以通过起始位置和结束位置描述Geotiff数据中一片特定的区域,也包含了一些常用的空间判断和空间交互方法,方便使用.

   • Window的尺寸是根据Segment尺寸计算而来的.理论上Window与Segment的尺寸关系会存在以下9种:

     
     
     • 红框为window,黑框为Segment
     • maxSize取默认值为256:
       • 112 代表maxSize >= seg * 2的情况
       • 192 代表 maxSize/2 < seg < maxSize的情况
       • 272 代表 maxSize < seg 的情况
     • 按对应关系,又可以分为以下4种:
       • 一个Segment只对应一个window(最理想,最简单的情况,表中192x192)
       • 一个Segment对应多个window(表中272x272/192x272/272x192)
       • 一个window对应多个Segment(表中112x112/112x192/192x112)
       • 多个window对应多个Segment(最复杂,表中112x272/272x112)
2. 根据计算的大小生成全部Window,这些Window集合所描述的全部区域,就是整个影像的区域.

我们可以提出两个问题:

1. 为什么Geotrellis要定义Window这样一个概念?
2. Segment本身也可以描述Geotiff中一片特定区域.为何不直接使用Segment的尺寸,而是要按一定方法重新规划尺寸呢?

这两个问题的解答来自下面的小节.

3.2.2 按分区限制对window集合分组

为什么要分组?我们不能忘记我们的目的:我们要为实际分布式读取数据规划读取窗口.

在spark中如何高效的计算?我们需要保证spark中每个分区的计算量相似,即保证每个分区处理的数据量是相似的.这就是对Window集合分组的意义:给每个分区分配合适的读取任务.

我们可以提出一个问题:怎么样分配才算合适呢?

这部分的主体是info.segmentLayout.partitionWindowsBySegments方法,整理后代码如下:

[代码位于GeoTiffInfoReader.scala]

// 按spark分区尺寸限制划分window
// 本例中partitionBytes使用spark默认的单个分区最大字节数
info.segmentLayout.partitionWindowsBySegments(windows, partitionBytes / math.max(info.cellType.bytes, 1))

def partitionWindowsBySegments(windows: Seq[GridBounds[Int]], maxPartitionSize: Long): Array[Array[GridBounds[Int]]] = {
    val partition = mutable.ArrayBuilder.make[GridBounds[Int]]
    partition.sizeHintBounded(128, windows)
    var partitionSize: Long = 0l
    var partitionCount: Long = 0l
    val partitions = mutable.ArrayBuilder.make[Array[GridBounds[Int]]]

    def finalizePartition() {
      val res = partition.result
      if (res.nonEmpty) partitions += res
      partition.clear()
      partitionSize = 0l
      partitionCount = 0l
    }

    def addToPartition(window: GridBounds[Int]) {
      partition += window
      partitionSize += window.size
      partitionCount += 1
    }
    
    // 描述整个影像区域的GridBounds
    val sourceBounds = GridBounds(0, 0, totalCols - 1, totalRows - 1)

    // 1. 对window进行重新排序
    val sorted = windows
      .filter(sourceBounds.intersects) // 只筛选与影像相交的window
      .map { window =>
        // 反算与window对应的Segment的序号(SegmentIndex)
        window -> getSegmentIndex(col = window.colMin, row = window.rowMin) 
      } 
      .sortBy(_._2) 
    
    // 2. 根据分区限制对Window集合进行分组
    for ((window, _) <- sorted) {
      if ((partitionCount == 0) || (partitionSize + window.size) < maxPartitionSize) {
        addToPartition(window)
      } else {
        finalizePartition()
        addToPartition(window)
      }
    }

    finalizePartition()
    partitions.result
  }

代码分为两个步骤:

1. 对筛选后的Window进行重新排序.
   • 根据window所代表的区域反算SegmentIndex,按SegmentIndex排序.
   • Window集合的默认顺序来自于listWindows,可能并非按SegmentIndex排序,SegmentIndex是数据存储的顺序,按这个顺序读取相当于顺序读取,理论上是最高效的.
2. 根据分区限制为Window集合分组
   • 将全部window划分为连续的若干组,每一组内全部Window描述区域数据所占的体积都最接近spark分区尺寸极限.

现在我们可以回答上节提出的问题了:

1. 为什么Geotrellis要定义Window这样一个概念?
   • Window是我们读取数据的粒度描述,它通过一定的策略保障了我们在每一个节点都处理接近分区上限的数据.没有这一中间概念,我们无法规划如何高效的分布式读取数据.
2. 定义Window的时候为何不直接使用Segment的尺寸,而是要按一定方法重新规划尺寸呢?
   • 直接原因就是:如果使用Segment原本的尺寸,在保证安全的情况下可能会造成巨大的空间浪费.
   • 在listwindow方法中我们确定了window的尺寸,无论是长还是宽,都在128-256之间.这与我们使用的默认分区大小128MB(128 * 1024 * 1024)和Tiff6.0限定尺寸必须为16的倍数是相关的,由于Window是我们规划读取数据的粒度,只能完整读取整个Window的数据,合适的Window尺寸能保证我们能在不超过分区尺寸限制的情况下尽可能多的读取数据.
   • 如果我们按照Segment的实际尺寸设定window的尺寸,极端情况下会出现这种情况:一个巨大的Window所包含数据的体积恰好比分区限制的一半大一点,因为要保证不超限制,又不能读取不完整的window,只能让该分区浪费接近一半的空间,只读取一个window.
   • 采用特定的策略计算windwos的长宽尺寸,将其控制在128-256之间,能保证浪费的内存空间最小.

最终,通过flatMap和一些其他操作,Geotrellis得到了一个RDD[(URI, Array[GridBounds[Int]])]对象,其含义很明确:

• 一个(URI,Array[GridBounds[Int]])元组描述某个影像的某一特定部分,这部分影像的体积能撑满一个分区空间.
• 若干个具有相同URI的(URI, Array[GridBounds[Int]])元组描述了某个影像的全部数据
  • 最终,它们将通过spark在不同的分区并行读取特定的部分.
• 全部(URI, Array[GridBounds[Int]])描述了全部影像及他们各自全部的分区.

至此,规划读取窗口的工作就做完了,可以进行实际读取了.

3.3 开始读取数据

我们再来回顾一下开始读取数据的代码:

val result = repartition.flatMap { case (uri, windows) =>
  // 因为发生了重分区,因此需要在每个分区重新读取影像信息
  val info = getGeoTiffInfo(uri.toString)
  // 读取数据的核心,rr是一个隐式的rasterReader对象
  rr.readWindows(windows, info, rasterReaderOptions).map { case (k, v) =>
    uriToKey(uri, k) -> v
  }
}

这部分的主体是rasterReader.readWindows方法,整理后代码如下:

[代码位于RasterReader.scala]

implicit def multibandGeoTiffReader = new RasterReader[Options, (ProjectedExtent, MultibandTile)] {
    def readWindows(gbs: Array[GridBounds[Int]], info: GeoTiffInfo, options: Options) = {
      // 创建一个多波段Geotiff对象,由于使用了惰性加载,因此此时没有实际读取任何数据
      val geoTiff = GeoTiffReader.geoTiffMultibandTile(info)
      val re = info.rasterExtent
      // 读取的核心:读取传入的一系列window对象所描述的区域的数据,存入内存
      geoTiff.crop(gbs.filter(_.intersects(geoTiff.dimensions))).map { case (gb, tile) =>
        (ProjectedExtent(re.extentFor(gb, clamp = false), options.crs.getOrElse(info.crs)), tile)
      }
    }
}

这部分的主体是geoTiff.crop方法,整理后代码如下:

[代码位于GeoTiffMultibandTile.scala]

 def crop(windows: Seq[GridBounds[Int]], bandIndices: Array[Int]): Iterator[(GridBounds[Int], ArrayMultibandTile)] = {
    
    val bandSubsetLength = bandIndices.length
    
    // 1. 准备阶段
    case class Chip(
      window: GridBounds[Int],
      bands: Array[MutableArrayTile],
      intersectingSegments: Int,
      var segmentsBurned: Int = 0
    )
    
    val chipsBySegment = mutable.Map.empty[Int, List[Chip]]
    val intersectingSegments = mutable.SortedSet.empty[(Int, Int)]
    
    for (window <- windows) {
      val segments: Array[(Int, Int)] =
        if (hasPixelInterleave)
          getIntersectingSegments(window).map(i => 0 -> i)
        else
          getIntersectingSegments(window, bandIndices)

      val bands = Array.fill(bandSubsetLength)(ArrayTile.empty(cellType, window.width, window.height))
      val chip = Chip(window, bands, segments.length)
      for (segment <- segments.map(_._2)) {
        val tail = chipsBySegment.getOrElse(segment, Nil)
        chipsBySegment.update(segment, chip :: tail)
      }
      intersectingSegments ++= segments
    }
    
    // 2. 实际读取数据
    val chips =
      if (hasPixelInterleave) {
        getSegments(intersectingSegments.map(_._2)).flatMap { case (segmentId, segment) =>
          burnPixelInterleave(segmentId, segment)
        }
      } else {
        val bandSegmentCount = segmentCount / bandCount
        val segmentBandMap = intersectingSegments.map { case (f, s) => (s, f) }.toMap
        val bandIndexToSubsetIndex = bandIndices.zipWithIndex.toMap

        getSegments(intersectingSegments.map(_._2)).flatMap { case (segmentId, segment) =>
          val bandIndex = segmentBandMap(segmentId)
          val subsetBandIndex = bandIndexToSubsetIndex(bandIndex)

          burnBandInterleave(segmentId, subsetBandIndex, segment)
        }
      }

    chips.map { chip =>
      chip.window -> ArrayMultibandTile(chip.bands)
    }
  }

代码分为两个步骤:

1. 准备阶段
2. 实际读取数据

按顺序我们先从准备阶段开始.

3.3.1 准备阶段

这里依旧延续了先规划任务,再实际执行计算的思路.我们传入的是Window对象,而根据上一篇的经验,实际要读取数据需要的是Segment对象.因此,我们需要先规划Window对象与Segment对象的映射关系.

代码分为两部分,我们首先看一下定义部分:

// 定义辅助对象/结构
case class Chip(
    window: GridBounds[Int],
    bands: Array[MutableArrayTile],
    intersectingSegments: Int,
    var segmentsBurned: Int = 0
)

val chipsBySegment = mutable.Map.empty[Int, List[Chip]]
val intersectingSegments = mutable.SortedSet.empty[(Int, Int)]

定义辅助对象/结构如下

• Chip:数据结构,是一个关联辅助类型.用于最终的数据读取
• chipsBySegment:映射表,记录与某一个SegmentIndex关联的全部Chip
• intersectingSegments:非重复数据集,记录与整片区域相关联的全部SegmentIndex.

实际执行部分:

// 遍历Window集合,填充上述的对象
for (window <- windows) {
  // 1. 计算与该window所描述范围关联的每个波段的SegmentIndex
  // segments代表[(波段Index,SegmentIndex)]
  val segments: Array[(Int, Int)] =
    if (hasPixelInterleave) // 要处理间隔方式造成的波段差异
      getIntersectingSegments(window).map(i => 0 -> i)
    else
      getIntersectingSegments(window, bandIndices)
  
  // 2. 创建chip对象
  val bands = Array.fill(bandSubsetLength)(ArrayTile.empty(cellType, window.width, window.height))
  val chip = Chip(window, bands, segments.length)
  
  // 3. 将Chip与Segment关联起来    
  for (segment <- segments.map(_._2)) {
    val tail = chipsBySegment.getOrElse(segment, Nil)
    chipsBySegment.update(segment, chip :: tail)
  }
  
  // 4. 记录全部涉及的Segment
  intersectingSegments ++= segments
}

代码分为四个步骤:

1. 计算与该Window所描述范围关联的每个波段的SegmentIndex.
   • 因为涉及到多波段数据读取,我们需要关注波段间隔带来的读取差异,可以在这里了解不同波段间隔的差异.
   • 该步骤有两个操作:
     • 计算与该Window相对应的全部SegmentIndex(根据上表,存在一个Window对应多个Segment的情况)
     • 因为计算出来的是第一个波段的SegmentIndex,还需要根据规律计算全部波段的SegmentIndex
   • 针对多波段,需要按对应的波段间隔规律计算每一个波段对应的SegmentIndex
2. 创建一个Chip对象.
   • bands属性为全部波段都预留了写入数据的位置
   • intersectingSegments属性记录了该Window对应多少个Segment
3. 将Chip与Segment关联起来
   • chipsBySegment以Segment为键,就需要解决多个Chip对应同一个Segment的情况.因此在这里要使用update方法,不断融合多个Chip,保证数据完整.
4. 记录全部涉及到的Segment
   • intersectingSegments是一个Set,因此不会包含重复数据
   • 最终我们可以根据intersectingSegments,遍历出所有对应的Chip,也就能读取到全部的数据了.

一切准备就绪,可以开始真正的读取阶段了.

3.3.2 开始读取

val chips =
  if (hasPixelInterleave) {
  
    // 从不同数据类型的数据源中根据SegmentIndex读取Segment
    getSegments(intersectingSegments.map(_._2)).flatMap { case (segmentId, segment) =>
      burnPixelInterleave(segmentId, segment)
    }
  } else {
    val bandSegmentCount = segmentCount / bandCount
    val segmentBandMap = intersectingSegments.map { case (f, s) => (s, f) }.toMap
    val bandIndexToSubsetIndex = bandIndices.zipWithIndex.toMap

    getSegments(intersectingSegments.map(_._2)).flatMap { case (segmentId, segment) =>
      val bandIndex = segmentBandMap(segmentId)
      val subsetBandIndex = bandIndexToSubsetIndex(bandIndex)

      burnBandInterleave(segmentId, subsetBandIndex, segment)
    }
  }

// 结果为[(GridBounds[Int], ArrayMultibandTile)]类型的对象
// 在后面的过程中,GridBounds[Int]被转换为ProjectedExtent类型,即我们最终需要的类型
chips.map { chip =>
  chip.window -> ArrayMultibandTile(chip.bands)
}

可以看出,这部分的主体是针对不同波段间隔的burnPixelInterleave方法.整理后代码如下:

// BIP类型
def burnPixelInterleave(segmentId: Int, segment: GeoTiffSegment): List[Chip] = {
  var finished: List[Chip] = Nil
  val segmentBounds = getGridBounds(segmentId)
  val segmentTransform = getSegmentTransform(segmentId)

  // 遍历chip
  for (chip <- chipsBySegment(segmentId)) {
    val gridBounds = chip.window
    // 获取影像中该window对应的位置
    val overlap = gridBounds.intersection(segmentBounds).get
    if (cellType.isFloatingPoint) {
      cfor(overlap.colMin)(_ <= overlap.colMax, _ + 1) { col =>
        cfor(overlap.rowMin)(_ <= overlap.rowMax, _ + 1) { row =>
          cfor(0)(_ < bandSubsetLength, _ + 1) { bandIndex =>
            // 针对BIP类型读取不同波段的数据
            val i = segmentTransform.gridToIndex(col, row, bandIndices(bandIndex))
            // 读取数据,因为数据采用懒加载,只有在这里才真正读取了数据,即每个分区都中读取了自己应该读取的那部分
            val v = segment.getDouble(i)
            chip.bands(bandIndex).setDouble(col - gridBounds.colMin, row - gridBounds.rowMin, v)
          }
        }
      }
    } else {
      cfor(overlap.colMin)(_ <= overlap.colMax, _ + 1) { col =>
        cfor(overlap.rowMin)(_ <= overlap.rowMax, _ + 1) { row =>
          cfor(0)(_ < bandSubsetLength, _ + 1) { bandIndex =>
            val i = segmentTransform.gridToIndex(col, row, bandIndices(bandIndex))
            val v = segment.getInt(i)
            chip.bands(bandIndex).set(col - gridBounds.colMin, row - gridBounds.rowMin, v)
          }
        }
      }
    }
    // 因为存在Window比Segment大的情况,因此需要保证每一个Chip读取的都是完整的数据
    chip.segmentsBurned += 1
    // 因此一个Chip可能会被多个Segment引用,只有当该Chip对应的所有Segment都被读取后,才算整个Chip读取完整
    if (chip.segmentsBurned == chip.intersectingSegments)
      finished = chip :: finished
  }
  // 因为intersectingSegments是非重复的,因此一个Segment读取完毕就可以保证不再有其他的Chip需要该Sement的数据,所以可以删掉
  chipsBySegment.remove(segmentId)
  // 返回的一定是完整读取数据的Chip
  finished
}

// BIL/BSQ类型
def burnBandInterleave(segmentId: Int, subsetBandIndex: Int, segment: GeoTiffSegment): List[Chip] = {
  var finished: List[Chip] = Nil
  val segmentBounds = getGridBounds(segmentId)
  val segmentTransform = getSegmentTransform(segmentId)

  for (chip <- chipsBySegment(segmentId)) {
    val gridBounds = chip.window
    val overlap = gridBounds.intersection(segmentBounds).get

    if (cellType.isFloatingPoint) {
      cfor(overlap.colMin)(_ <= overlap.colMax, _ + 1) { col =>
        cfor(overlap.rowMin)(_ <= overlap.rowMax, _ + 1) { row =>
          val i = segmentTransform.gridToIndex(col, row)
          val v = segment.getDouble(i)
          chip.bands(subsetBandIndex).setDouble(col - gridBounds.colMin, row - gridBounds.rowMin, v)
        }
      }
    } else {
      cfor(overlap.colMin)(_ <= overlap.colMax, _ + 1) { col =>
        cfor(overlap.rowMin)(_ <= overlap.rowMax, _ + 1) { row =>
          val i = segmentTransform.gridToIndex(col, row)
          val v = segment.getInt(i)
          chip.bands(subsetBandIndex).set(col - gridBounds.colMin, row - gridBounds.rowMin, v)
        }
      }
    }

    chip.segmentsBurned += 1
    if (chip.segmentsBurned == chip.intersectingSegments)
      finished = chip :: finished
  }

  chipsBySegment.remove(segmentId)
  finished
}

至此,Geotiff中的数据,就在每个分区中按一定的规则读取到内存中去了

4. 回顾与总结

如果从目的反推行动,将会是这样:

1. 为了要分布式的读取数据,就要规划每个节点应该读取哪些数据.
2. 为了规划每个节点应该读取哪些数据,我们需要对需要处理的数据进行分组.
3. 为了进行分组,我们要知道如何更好的分组.
4. 为了知道最合理的分组策略,我们需要从影像的元数据获取相关信息
5. 为了读取影像的元数据,我们需要知道有哪些影像.
6. 因此我们首先获取了全部待处理影像的路径.

回顾一下整个读取流程,可以说都是在践行先规划任务,再实际执行计算的思路,合理的任务规划与分配对分布式计算最终效率影像极大.