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

前言

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

1、寫入數據——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 實現了三種，分別爲 S三、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 文件。

2、讀取數據

數據作了 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 函數將其轉爲瓦片。

3、總結

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

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