geotrellis使用（三十八）COG 寫入和讀取

前言

上一篇中簡單介紹了 COG 的概念和 Geotrellis 中引入 COG 的原因及簡單的原理，本文為大家介紹如何在 Geotrellis 中使用 COG 來寫入和讀取 GeoTIFF數據。

一、寫入數據——ETL

1.1 實現方案

其實這與之前的普通 ETL 操作在概念上是相似的，都是將原始數據轉換成系統能用的數據的過程，這是寬泛的 ETL 的定義。在 Geotrellis 中實現很簡單，與之前代碼基本一致，只要切換一下 writer 類型以及最后建立金字塔額時候略有不同。實現方案如下：

val inputPath = "file://" + new File("in.tif").getAbsolutePath
val outputPath = "/your/catalog/path"

def main(args: Array[String]): Unit = {
    // Setup Spark to use Kryo serializer.
    val conf =
      new SparkConf()
        .setMaster("local[*]")
        .setAppName("Spark Tiler")
        .set("spark.serializer", "org.apache.spark.serializer.KryoSerializer")
        .set("spark.kryo.registrator", "geotrellis.spark.io.kryo.KryoRegistrator")
    
    val sc = new SparkContext(conf)
    try {
      run(sc)
    } finally {
      sc.stop()
    }
}

def run(implicit sc: SparkContext) = {
    val inputRdd: RDD[(ProjectedExtent, MultibandTile)] =
    sc.hadoopMultibandGeoTiffRDD(inputPath)
    
    val (_, rasterMetaData) =
    TileLayerMetadata.fromRdd(inputRdd, FloatingLayoutScheme(256))
    
    val tiled: RDD[(SpatialKey, MultibandTile)] =
    inputRdd
      .tileToLayout(rasterMetaData.cellType, rasterMetaData.layout, Bilinear)
      .repartition(100)

    val layoutScheme = ZoomedLayoutScheme(WebMercator, tileSize = 256)
    
    val (zoom, reprojected): (Int, RDD[(SpatialKey, MultibandTile)] with Metadata[TileLayerMetadata[SpatialKey]]) =
      MultibandTileLayerRDD(tiled, rasterMetaData)
        .reproject(WebMercator, layoutScheme, Bilinear)
    
    val attributeStore = FileAttributeStore(outputPath)
    
    val writer = FileCOGLayerWriter(attributeStore)
    
    val layerName = "layername"

    val cogLayer =
      COGLayer.fromLayerRDD(
        reprojected,
        zoom,
        compression = NoCompression,
        maxTileSize = 4096
      )
    
    val keyIndexes =
      cogLayer.metadata.zoomRangeInfos.
        map { case (zr, bounds) => zr -> ZCurveKeyIndexMethod.createIndex(bounds) }.
        toMap
    
    writer.writeCOGLayer(layerName, cogLayer, keyIndexes)
}

執行 main 函數即可實現 COG 方式的 ETL 操作，其他部分與之前介紹過的 ingest 相同，主要區別在於 writer，此處為 FileCOGLayerWriter 實例，從名字可以看出這是一個文件系統的 COG writer，目前 Geotrellis 實現了三種，分別為 S3、Hadoop、File，這三種后端理論上都是對大量小文件支持不好的。

1.2 背后邏輯

下面來詳細分析一下 Geotrellis 中 COG 實現原理。

首先看上面的 cogLayer 對象：

val cogLayer =
  COGLayer.fromLayerRDD(
    reprojected,
    zoom,
    compression = NoCompression,
    maxTileSize = 4096
  )

cogLayer 對象是 COGLayer 實例，fromLayerRDD 源碼如下：

def fromLayerRDD[
    K: SpatialComponent: Ordering: JsonFormat: ClassTag,
    V <: CellGrid: ClassTag: ? => TileMergeMethods[V]: ? => TilePrototypeMethods[V]: ? => TileCropMethods[V]: GeoTiffBuilder
](
    rdd: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
    baseZoom: Int,
    compression: Compression = Deflate,
    maxTileSize: Int = 4096,
    minZoom: Option[Int] = None
): COGLayer[K, V] = {

    if(minZoom.getOrElse(Double.NaN) != baseZoom.toDouble) {
      if(rdd.metadata.layout.tileCols != rdd.metadata.layout.tileRows) {
        sys.error("Cannot create Pyramided COG layer for non-square tiles.")
      }
    
      if(!isPowerOfTwo(rdd.metadata.layout.tileCols)) {
        sys.error("Cannot create Pyramided COG layer for tile sizes that are not power-of-two.")
      }
    }
    
    val layoutScheme =
      ZoomedLayoutScheme(rdd.metadata.crs, rdd.metadata.layout.tileCols)
    
    if(rdd.metadata.layout != layoutScheme.levelForZoom(baseZoom).layout) {
      sys.error(s"Tile Layout of layer and ZoomedLayoutScheme do not match. ${rdd.metadata.layout} != ${layoutScheme.levelForZoom(baseZoom).layout}")
    }
    
    val keyBounds =
      rdd.metadata.bounds match {
        case kb: KeyBounds[K] => kb
        case _ => sys.error(s"Cannot create COGLayer with empty Bounds")
      }
    
    val cogLayerMetadata: COGLayerMetadata[K] =
      COGLayerMetadata(
        rdd.metadata.cellType,
        rdd.metadata.extent,
        rdd.metadata.crs,
        keyBounds,
        layoutScheme,
        baseZoom,
        minZoom.getOrElse(0),
        maxTileSize
      )
    
    val layers: Map[ZoomRange, RDD[(K, GeoTiff[V])]] =
      cogLayerMetadata.zoomRanges.
        sorted(Ordering[ZoomRange].reverse).
        foldLeft(List[(ZoomRange, RDD[(K, GeoTiff[V])])]()) { case (acc, range) =>
          if(acc.isEmpty) {
            List(range -> generateGeoTiffRDD(rdd, range, layoutScheme, cogLayerMetadata.cellType, compression))
          } else {
            val previousLayer: RDD[(K, V)] = acc.head._2.mapValues { tiff =>
              if(tiff.overviews.nonEmpty) tiff.overviews.last.tile
              else tiff.tile
            }
    
            val tmd: TileLayerMetadata[K] = cogLayerMetadata.tileLayerMetadata(range.maxZoom + 1)
            val upsampledPreviousLayer =
              Pyramid.up(ContextRDD(previousLayer, tmd), layoutScheme, range.maxZoom + 1)._2
    
            val rzz = generateGeoTiffRDD(upsampledPreviousLayer, range, layoutScheme, cogLayerMetadata.cellType, compression)
    
            (range -> rzz) :: acc
          }
        }.
        toMap
    
    COGLayer(layers, cogLayerMetadata)
}

此函數返回類型正是 COGLayer，其兩個屬性分別為 layers 和 cogLayerMetadata。

cogLayerMetadata 是 COGLayerMetadata 對象，表示 COG 層的元數據信息，包含每層對應的瓦片范圍等，這個與傳統的元數據很接近，唯一不同的在於此處使用了 ZommRange 的概念，即“ 1 層”可能對應多個 Zoom，而不再是 1 對 1 的關系，這樣能夠更進一步的節省存儲空間，在處理速度和存儲空間上做了綜合考慮。

layers 是 Map[ZoomRange, RDD[(K, GeoTiff[V])]] 對象，ZoomRange 即為上述元數據中的每層的 zoom 最大和最小值，RDD[(K, GeoTiff[V])] 是 spark rdd 對象，即每一個層級范圍對應一個 Tiff 對象，從此可以看出，COG 方式 ETL 后每層存儲的不再是 Tile，而是 Tiff 文件，這個 Tiff 文件是 COG 類型的，當用戶請求某個瓦片的時候直接從對應的 Tiff 文件中讀出需要的部分即可。

最后調用 writer.writeCOGLayer(layerName, cogLayer, keyIndexes) 即可將元數據信息和 Tiff 數據寫入相應的位置，完成 ETL 過程。

此處還需要注意的是為了防止單個 Tiff 文件過大， Geotrellis 對每一層進行了分割處理，這樣每一層可能會得到多個 Tiff 文件，而為了達到 COG 的真實效果，又引入了 GDAL 中 VRT 的概念（參見http://www.gdal.org/gdal_vrttut.html），其中很詳細的講述了 VRT 的格式和意義，簡單來說 VRT 就是將多個 Tiff 文件合並成一個虛擬的 Tiff 文件。

二、讀取數據

數據做了 ETL 后，就可以讀取出來並進行相應的處理。

2.1 實現方案

其實現方式也基本與傳統方式相同，代碼如下：

val catalogPath = new java.io.File("/your/catalog/path").getAbsolutePath
val fileValueReader = FileCOGValueReader(catalogPath)
val key = SpatialKey(x, y)
val tile = fileValueReader.reader(LayerId("layername", z)).read(key)

這樣就能根據瓦片的 x、y 編號和 z（zoom）取到對應的瓦片。

2.2 原理

主要代碼在 COGValueReader 類中的 baseReader 方法中 read 方法，如下：

GeoTiffReader[V].read(uri, decompress = false, streaming = true)
  .getOverview(overviewIndex)
  .crop(gridBounds)
  .tile

傳統方式存儲的是切割好的瓦片，可以直接定位到確定的瓦片，這里是完全符合 COG 方式的讀取方式。getOverview 獲取到對應層（z）的 Tiff 文件，crop 對 Tiff 根據需要的范圍（x、y）進行切割，tile 函數將其轉為瓦片。

三、總結

本文介紹了如何在 Geotrellis 中如何進行 COG 方式的 ETL 操作，實現了全新的數據寫入和讀取方式。

Geotrellis系列文章鏈接地址http://www.cnblogs.com/shoufengwei/p/5619419.html