How it works(26) Geotrellis是如何在Spark上计算的(G) 导出结果到GeoTiff文件

默而识之者

已于 2022-05-11 08:21:57 修改

阅读量737

点赞数

分类专栏：源码阅读 geotrellis 文章标签： geotrellis

于 2022-05-11 08:18:56 首次发布

本文链接：https://blog.csdn.net/u011609113/article/details/124702024

版权

源码阅读同时被 2 个专栏收录

26 篇文章 15 订阅

订阅专栏

geotrellis

11 篇文章 1 订阅

订阅专栏

1. 引入

在过去几章,我们从NDVI计算入手,深入到Geotrellis中了解了内置的各种算子.如今我们回归最初计算NDVI的DEMO代码,看一下如何将计算后的NDVI结果以GeoTiff的格式输出:

val ndviTiledRDD: TileLayerRDD[SpatialKey] = // ... 省略计算ndvi的步骤
// 1. 将TileLayerRDD[SpatialKey]对象拼接为一个Raster[Tile]对象
val raster: Raster[Tile] = ndviTiledRDD.stitch()
// 2. 最终导出Geotiff文件到本地
GeoTiff(raster, metadata.crs).write("/temp/t.tif")

将计算结果输出为Geotiff文件需要两个步骤:

将Rdd对象拼接,并转换为Raster对象
通过Raster对象创建Geotiff对象,并写入本地

我们从第一步开始.

2. 拼接RDD

NDVI计算后的结果可视为若干小Tile组成的列表.为了输出为Geotiff,我们需要将它们镶嵌到一个巨大的Tile中,需要调用stitch方法:

代码位于[geotrellis/saprk/stitch/StitchRDDMethods.scala]

abstract class SpatialTileLayoutRDDStitchMethods[
  V <: CellGrid[Int]: Stitcher: * => TilePrototypeMethods[V],
  M: GetComponent[*, LayoutDefinition]
] extends MethodExtensions[RDD[(SpatialKey, V)] with Metadata[M]] {

  def stitch(): Raster[V] = {
    // ... 省略 ...
  }
  // ... 省略 ...
}

stitch方法依旧是通过定义隐式转换对RDD对象进行扩展,因为原本Rdd对象并不存在stitch方法:

代码位于[geotrellis/saprk/stitch/Implicits.scala]

trait Implicits {
  implicit class withSpatialTileLayoutRDDMethods[
    V <: CellGrid[Int]: Stitcher: * => TilePrototypeMethods[V],
    M: GetComponent[*, LayoutDefinition]
  ](
    val self: RDD[(SpatialKey, V)] with Metadata[M]
  ) extends SpatialTileLayoutRDDStitchMethods[V, M]

  // ... 省略 ...
}

我们来梳理一下stitch方法通过隐式转换调用的顺序:

ndviTiledRDD的实际类型为RDD[(SpatialKey, Tile)] with Metadata[TileLayerMetadata[SpatialKey]]
RDD类型中没有名为stitch的方法,而在隐式类withSpatialTileLayoutRDDMethods中具有
Scala开始在上下文中寻找RDD向withSpatialTileLayoutRDDMethods类型转换是否可行
在这里,withSpatialTileLayoutRDDMethods的调用使用了Value calsses定义方法,因此变为需要检查Tile类型是否符合定义中对于V的限制和TileLayerMetadata[SpatialKey]是否符合定义中对于M的限制.
经验证,存在Tile向Stitcher[Tile]的转换,其他的类型限制在前文中也验证过,即类型符合定义中对类型的限制.因此可以调用withSpatialTileLayoutRDDMethods类中的stitch方法.

我们再来看stitch方法的实体:

代码位于[geotrellis/saprk/stitch/StitchRDDMethods.scala]

def stitch(): Raster[V] = {
    // 1. 获取拼接后的大Tile,第一个瓦片的行列号及其偏移尺寸
    val (tile, (kx, ky), (offsx, offsy)) = TileLayoutStitcher.stitch(self)
    // 2. 从TileLayerMetadata中获取数据的布局及坐标转换
    val layout = self.metadata.getComponent[LayoutDefinition]
    val mapTransform = layout.mapTransform
    // 3. 从行列号反算回坐标
    val nwTileEx = mapTransform(kx, ky)
    val base = nwTileEx.southEast
    // 4. 得到实际的左上角坐标
    val (ulx, uly) = (base.x - offsx.toDouble * layout.cellwidth, base.y + offsy * layout.cellheight)
    // 5. 得到最终的Raster[V]对象
    Raster(tile, Extent(ulx, uly - tile.rows * layout.cellheight, ulx + tile.cols * layout.cellwidth, uly))
}

拼接过程包含了两部分:

把若干小Tile拼接为一个大Tile
获取拼接后的大Tile的尺寸与其他信息

第一部分的逻辑位于TileLayoutStitcher.stitch方法中.

2.1 拼接Tile

代码位于[geotrellis/layer/stitch/TileLayoutStitcher.scala]

object TileLayoutStitcher {
  def stitch[
    V <: CellGrid[Int]: Stitcher
  ](tiles: Traversable[(Product2[Int, Int], V)]): (V, (Int, Int), (Int, Int)) = {
    val colWidths = mutable.Map.empty[Int, Int]
    val rowHeights = mutable.Map.empty[Int, Int]
    // 1. 验证传入的tiles
    tiles.foreach{ case (key, tile) =>
      // 验证当前tile的宽度与其他同列tile的宽度是否一致
      val curWidth = colWidths.getOrElseUpdate(key._1, tile.cols)
      assert(curWidth == tile.cols, "Tiles in a layout column must have the same width")
      val curHeight = rowHeights.getOrElseUpdate(key._2, tile.rows)
      assert(curHeight == tile.rows, "Tiles in a layout row must have the same height")
    }
    
    // 2. 整理tile的布局,计算实际的行列偏移量
    // 如果colWidths存储的是[(0->256),(2->256),(1->256)]
    // colPos的结果就是[(0->0),(1->256),(2->512)]
    // width的结果就是768
    val (colPos, width) = colWidths.toSeq.sorted.foldLeft( (Map.empty[Int, Int], 0) ){
      case ((positions, acc), (col, w)) => (positions + (col -> acc), acc + w)
    }
    val (rowPos, height) = rowHeights.toSeq.sorted.foldLeft( (Map.empty[Int, Int], 0) ){
      case ((positions, acc), (row, h)) => (positions + (row -> acc), acc + h)
    }
    // 3. 拼接全部tile
    // 显式定义隐式对象
    val stitcher = implicitly[Stitcher[V]]
    val result = stitcher.stitch(tiles.map{ case (key, tile) => tile -> (colPos(key._1), rowPos(key._2)) }.toIterable, width, height)
    // 最终返回相关信息
    val (minx, miny) = (colWidths.keys.min, rowHeights.keys.min)
    (result, (minx, miny), (colWidths(minx), rowHeights(miny)))
  }
}

拼接Tile分为3步:

验证传入的全部tile是否是整齐的,这个要求其实比较宽松:
- 不要求tile具有统一的大小,但同一行的tile需要具有一样的高度,同一列的tile需要具有一样的宽度
- 不要求tile必须紧密排列,但因为拼接的时候会转换为紧凑的形式,因此如果tile为非连续,可能最终得不到预期的效果
整理tile的布局信息,为后面拼接做准备
正式开始拼接.主要是调用了stitcher.stitch方法:

代码位于[geotrellis/raster/stitch/Stitcher.scala]

implicit object TileStitcher extends Stitcher[Tile] {
    def stitch(pieces: Iterable[(Tile, (Int, Int))], cols: Int, rows: Int): Tile = {
      // 按拼接后的尺寸申请一个大Tile
      val result = ArrayTile.empty(pieces.head._1.cellType, cols, rows)
      // updateColMin/updateRowMin为我们计算后的实际偏移量
      // 因此每一个局部Tile将会被添加到大Tile中的正确位置上
      for((tile, (updateColMin, updateRowMin)) <- pieces) {
        result.update(updateColMin, updateRowMin, tile)
      }
      result
    }
}

2.2 计算拼接后Tile的信息

拼接完Tile后,我们还需要得出地理范围.前文我们说过,为了进行快速的计算,Geotrellis将实际的坐标值转换为行列号(即SpatialKey),而行列号与坐标值互转的信息则存储在元数据中(即Metadata[TileLayerMetadata[SpatialKey]]).

在索引转换篇中,我们已经讨论了如何从坐标信息转换到行列号信息.我们通过调用mapTransform的ExtentToKey方法,将实际的坐标范围转换为行列号.此时我们要反过来,调用KeyToExtent,将行列号转换回坐标范围.

整理上文的SpatialTileLayoutRDDStitchMethods.stitch方法,获得拼接后Tile坐标范围的逻辑如下:

调用TileLayoutStitcher.stitch方法,得到了左上角第一个瓦片的行列号和宽与高
从元数据中获取mapTransform对象,为了进行行列号到坐标的转换
反算第一个瓦片的行列号代表的坐标值,得到第一个瓦片右下角的坐标值
根据第一个瓦片的宽高和像元尺寸,计算得到左上角坐标
根据左上角坐标,像元尺寸和拼接后Tile的宽高,得到拼接后Tile的坐标范围

至此,有了拼接后的Tile及其范围,我们就能得到预期中的Raster[Tile]对象了.

3. 输出Geotiff

输出Geotiff的逻辑较为简单,即先创建一个GeoTiff对象,再将其写入本地文件:

GeoTiff(raster, metadata.crs).write("/temp/t.tif")

我们先看一下,创建Geotiff对象的过程中发生了什么.

3.1 创建Geotiff对象

创建Geotiff对象中最关键的一步是将便于计算的Tile形式转换回Geotiff的内部结构Segment:

代码位于[geotrellis/raster/io/geotiff/SinglebandGeoTiff.scala]

// 在本例中,ndvi计算结果为单波段数据,因此以SinglebandGeoTiff对象为例
case class SinglebandGeoTiff(
  tile: Tile,
  extent: Extent,
  crs: CRS,
  tags: Tags,
  options: GeoTiffOptions,
  overviews: List[GeoTiff[Tile]] = Nil
) extends GeoTiff[Tile] {
    val cellType = tile.cellType

    // ... 省略

    def imageData: GeoTiffImageData =
        tile match {
          case gtt: GeoTiffTile => gtt
          // tile对象转换扩展成了GeoTiffImageData的GeoTiffTile对象
          case _ => tile.toGeoTiffTile(options)
        }
 
    // ... 省略
 
}

在读取元数据篇中提到,我们使用诸如LazySegmentBytes或ArraySegmentBytes的方式,将Geotiff的二进制流读取为Segement对象,如今转出Geotiff,则需要从Tile对象中得到Segment.

tile.toGeoTiffTile的核心代码如下:

代码位于[geotrellis/raster/io/geotiff/GeoTiffTile.scala]

object GeoTiffTile
{
  // ... 省略
  
  def apply(tile: Tile, options: GeoTiffOptions): GeoTiffTile = {
    // 获取数据类型
    val bandType = BandType.forCellType(tile.cellType)
    // 获取Segment布局.默认为256x256的瓦片状布局(可选条带状)
    val segmentLayout = GeoTiffSegmentLayout(tile.cols, tile.rows, options.storageMethod, BandInterleave, bandType)
    // 计算原始大Tile以256x256的瓦片分割后的Segment数量
    val segmentCount = segmentLayout.tileLayout.layoutCols * segmentLayout.tileLayout.layoutRows
    // 默认采用无压缩模式
    val compressor = options.compression.createCompressor(segmentCount)
    // 先将原始Tile切分为Segment布局
    val segmentBytes = Array.ofDim[Array[Byte]](segmentCount)
    val segmentTiles: Seq[Tile] =
      options.storageMethod match {
        case _: Tiled => tile.split(segmentLayout.tileLayout)
        case _: Striped => tile.split(segmentLayout.tileLayout, Split.Options(extend = false))
      }
    // 将切割后的每一块Tile转换为字节数据
    cfor(0)(_ < segmentCount, _ + 1) { i =>
      val bytes = segmentTiles(i).toBytes
      segmentBytes(i) = compressor.compress(bytes, i)
    }

    apply(new ArraySegmentBytes(segmentBytes), compressor.createDecompressor, segmentLayout, options.compression, tile.cellType)
  }
}

将原始Tile切分为Segment布局的核心逻辑如下:

代码位于[geotrellis/raster/split/SinglebandTileSplitMethods.scala]

def split(tileLayout: TileLayout, options: Options): Seq[Tile] = {
    // 默认尺寸为256x256
    val tileCols = tileLayout.tileCols
    val tileRows = tileLayout.tileRows
    // 申请预期的存储空间
    val tiles = Array.ofDim[Tile](tileLayout.layoutCols * tileLayout.layoutRows)
    cfor(0)(_ < tileLayout.layoutRows, _ + 1) { layoutRow =>
      cfor(0)(_ < tileLayout.layoutCols, _ + 1) { layoutCol =>
        // 计算当前处理的行列号
        val firstCol = layoutCol * tileCols
        val lastCol = {
          val x = firstCol + tileCols - 1
          if(!options.extend && x > self.cols - 1) self.cols - 1
          else x
        }
        val firstRow = layoutRow * tileRows
        val lastRow = {
          val x = firstRow + tileRows - 1
          if(!options.extend && x > self.rows - 1) self.rows - 1
          else x
        }
        // 切出指定行列号范围的Tile
        val gb = GridBounds(firstCol, firstRow, lastCol, lastRow)
        tiles(layoutRow * tileLayout.layoutCols + layoutCol) =
          if(options.cropped) CroppedTile(self, gb)
          else CroppedTile(self, gb).toArrayTile
      }
    }
    
    tiles
}

3.2 导出Geotiff文件

导出文件实际调用了GeotiffWriter.write方法:

代码位于[geotrellis/raster/io/geotiff/writer/GeoTiffWriter.scala]

object GeoTiffWriter {
  def write(geoTiff: GeoTiffData, path: String, optimizedOrder: Boolean): Unit = {
    val fos = new FileOutputStream(new File(path))
    try {
      val dos = new DataOutputStream(fos)
      try {
        new GeoTiffWriter(geoTiff, dos).write(optimizedOrder)
      } finally {
        dos.close
      }
    } finally {
      fos.close
    }
  }
}

class GeoTiffWriter(geoTiff: GeoTiffData, dos: DataOutputStream) {
  // 初始化转换为Byte的方法,区分大小端,因为这涉及最初的字节该填充什么值
  implicit val toBytes: ToBytes =
    if(geoTiff.imageData.decompressor.byteOrder == ByteOrder.BIG_ENDIAN) {
      BigEndianToBytes
    } else {
      LittleEndianToBytes
    }
  
  // 拼接全部IFD列表,即Geotiff本身的IFD及其全部金字塔的IFD
  // IFD即图像文件目录,实际存储了TIff文件的各种元数据和图像实体等信息,每层金字塔都具有独立的IFD
  lazy val IFDs: List[GeoTiffData] = geoTiff :: geoTiff.overviews

  // 封装向数据流中写入各种类型的方法,即写入数据同时推动指针向前
  var index: Int = 0
  def writeByte(value: Byte): Unit = { dos.writeByte(value.toInt); index += 1 }
  def writeBytes(value: Array[Byte]): Unit =  { dos.write(value, 0, value.length); index += value.length }

  def writeShort(value: Int): Unit = { writeBytes(toBytes(value.toShort)) }
  def writeInt(value: Int): Unit = { writeBytes(toBytes(value)) }
  def writeLong(value: Long): Unit = { writeBytes(toBytes(value)) }
  def writeFloat(value: Float): Unit = { writeBytes(toBytes(value)) }
  def writeDouble(value: Double): Unit = { writeBytes(toBytes(value)) }
  
  // 向原始数据及其金字塔追加TIFF_TAG信息
  // 对于TIFF文件,其内嵌的金字塔也需要有对应的整套TIFF_TAG
  private def append(list: List[GeoTiffData]): Unit = {
    val overviewsIter = (geoTiff +: geoTiff.overviews).toIterator
    overviewsIter.foreach(append(_, !overviewsIter.hasNext))
  }
  
  // 核心方法
  // 向数据流中追加某个GeoTiff数据的TIFF_TAG
  private def append(geoTiff: GeoTiffData, last: Boolean = true): Unit = {
    // 获取该geotiff的全部TIFF_TAG
    val (fieldValues, offsetFieldValueBuilder) = TiffTagFieldValue.collect(geoTiff)
    // 获取Geotiff分片信息
    val segments = geoTiff.imageData.segmentBytes
    val segmentCount = segments.size
    val segmentBytesCount = (0 until segmentCount).map(segments.getSegmentByteCount).sum
    
    // 每一个TIFF_TAG字段描述的长度为2+2+4+4=12.详细的描述可见[读取元数据]篇的3.4节
    // 额外+1是为了补足后面动态生成的与图像有关的TIFF_TAG
    val tagFieldByteCount = (fieldValues.length + 1) * 12 
    
    // 统计全部TIFF_TAG内容的长度
    val tagDataByteCount = {
      var s = 0
      cfor(0)(_ < fieldValues.length, _ + 1) { i =>
        val len = fieldValues(i).value.length
        // 只有长度大于4的内容会占用独立的存储位置,否则将被直接存储在字段描述中
        if(len > 4) {
          s += fieldValues(i).value.length
        }
      }
      // 为 offsetFieldValue 内容预留位置
      if(segmentCount > 1) {
        s += (segmentCount * 4)
      }
      s
    }

    // 计算与偏移相关(offsetFieldValue)的TIFF_TAG
    val offsetFieldValue = {
      val imageDataStartOffset =
        index +
          2 + 
          4 + 
          tagFieldByteCount + tagDataByteCount

      val offsets = Array.ofDim[Int](segmentCount)
      var offset = imageDataStartOffset
      segments.getSegments(0 until segmentCount).foreach { case (i, segment) =>
        offsets(i) = offset
        offset += segment.length
      }
      offsetFieldValueBuilder(offsets)
    }

    // 将TIFF_TAG按标准TAG_CODE排序
    val sortedTagFieldValues = (offsetFieldValue :: fieldValues.toList).sortBy(_.tag).toArray

    // 写入TIFF_TAG数量
    writeShort(sortedTagFieldValues.length)

    // 计算TIFF_TAG内容偏移
    val tagDataStartOffset =
      index +
        4 + // 存储下一个IFD地址所预留的位置
        tagFieldByteCount

    // 写入全部TIFF_TAG描述信息
    var tagDataOffset = tagDataStartOffset
    cfor(0)(_ < sortedTagFieldValues.length, _ + 1) { i =>
      val TiffTagFieldValue(tag, fieldType, length, value) = sortedTagFieldValues(i)
      // 写入TIFF_TAG描述
      writeShort(tag)
      writeShort(fieldType)
      writeInt(length)
      // 长度大于4存储在TIFF_TAG内容存储区
      if(value.length > 4) {
        // 指向偏移位置
        writeInt(tagDataOffset)
        tagDataOffset += value.length
      } else {
        // 否则直接存储在TIFF_TAG描述中
        var i = 0
        while(i < value.length) {
          writeByte(value(i))
          i += 1
        }
        // 不足4位的要补齐
        while(i < 4) {
          writeByte(0.toByte)
          i += 1
        }
      }
    }

    // 写入结束标志位
    if(last) writeInt(0)
    else writeInt(tagDataOffset + segmentBytesCount)

    // 在内容存储区写入长度大于4的TIFF_TAG内容
    cfor(0)(_ < sortedTagFieldValues.length, _ + 1) { i =>
      val TiffTagFieldValue(tag, fieldType, length, value) = sortedTagFieldValues(i)
      if(value.length > 4) {
        writeBytes(value)
      }
    }

    // 写入图像
    segments.getSegments(0 until segmentCount).foreach { case (_, segment) =>
      writeBytes(segment)
    }
  }

  def write(optimizedOrder: Boolean = false): Unit = {
    // 在最开始写入区分大小端的标志位
    if (geoTiff.imageData.decompressor.byteOrder == ByteOrder.BIG_ENDIAN) {
      val m = 'M'.toByte
      writeByte(m)
      writeByte(m)
    } else {
      val i = 'I'.toByte
      writeByte(i)
      writeByte(i)
    }
    // 写入TIFF头标志位,标明是普通Tiff还是BigTiff
    // 这里默认为普通Tiff
    writeShort(Tiff.code.toShort)
    // 写入开始记录TIFF_TAG的标志位,即在第4个位置写入整数4,此时index也变为4
    // 该位置仅在普通Tiff类型时为4
    writeInt(index + 4)
    // 写入全部TIFF_TAG和数据本体
    append(IFDs)
    dos.flush()
  }
}

在TiffTagFieldValue.collect(geoTiff)方法中导出的TIFF_TAG整理如下:

名称	说明	条件
ImageWidth	图像宽度	-
ImageLength	图像长度	-
BitsPerSample	每个像元的bit长度,可以区分不同数据类型	-
Compression	描述压缩形式	-
Predictor	预测器,与压缩相关	-
PhotometricInterpretation	与色彩空间描述相关	-
SamplesPerPixel	波段数	-
SampleFormat	采样格式	-
PlanarConfiguration	存储方式,波段间隔还是像素间隔	-
NewSubfileType	新子文件类型,默认为空	-
ColorMap	色彩空间描述	指定了自定义调色板
GDAL_NODATA	GDAL专用的Nodata表达式	指定了Nodata值
ModelPixelScaleTag	定义仿射变换	-
ModelTiepointTag	仿射变换原点,与ModelPixelScaleTag配合使用	-
GeoKeyDirectoryTag	地理相关TAG扩展,记录CRS等信息	-
GeoDoubleParamsTag	地理相关TAG扩展,存储GeoKeyDirectoryTag中定义的Double类型的TAG	-
DateTime	日期	指定了了timeTag值
GDAL_METADATA	GDAL使用的Metadata	-
TileWidth,TileLength,TileByteCounts,TileOffsets	瓦片式布局参数	Tiff文件以瓦片式布局存储
RowsPerStrip,StripByteCounts,StripOffsets	条带式布局参数	Tiff文件以条带式布局存储

4. 总结

至此,我们就完善了在Geotrellis中,Geotiff文件从读取到处理,最终导出的全部过程.

虽然在Geotrellis中还有许多领域与功能没有涉及,但对数据的处理思路,其实都是一脉相承的.

默而识之者

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
How it works(26) Geotrellis是如何在Spark上计算的(G) 导出结果到GeoTiff文件

1. 引入在过去几章,我们从NDVI计算入手,深入到Geotrellis中了解了内置的各种算子.如今我们回归最初计算NDVI的DEMO代码,看一下如何将计算后的NDVI结果以GeoTiff的格式输出:val ndviTiledRDD: TileLayerRDD[SpatialKey] = // ... 省略计算ndvi的步骤// 1. 将TileLayerRDD[SpatialKey]对象拼接为一个Raster[Tile]对象val raster: Raster[Tile] = ndviTiledR
复制链接

扫一扫