How it works(26) Geotrellis是如何在S

作者: 默而识之者 | 来源:发表于2022-05-11 08:12 被阅读0次

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中还有许多领域与功能没有涉及,但对数据的处理思路,其实都是一脉相承的.

How it works(26) Geotrellis是如何在S

1. 引入

2. 拼接RDD

2.1 拼接Tile

2.2 计算拼接后Tile的信息

3. 输出Geotiff

3.1 创建Geotiff对象

3.2 导出Geotiff文件

4. 总结

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

GIS后端

Geotrellis