三十八. 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”)<span class="hljs-keyword">val</span> sc = <span class="hljs-keyword">new</span> <span class="hljs-type">SparkContext</span>(conf)
<span class="hljs-keyword">try</span> {run(sc)
} <span class="hljs-keyword">finally</span> {sc.stop()
}
}def run(implicit sc: SparkContext) = {
val inputRdd: RDD[(ProjectedExtent, MultibandTile)] =
sc.hadoopMultibandGeoTiffRDD(inputPath)<span class="hljs-keyword">val</span> (_, rasterMetaData) =
<span class="hljs-type">TileLayerMetadata</span>.fromRdd(inputRdd, <span class="hljs-type">FloatingLayoutScheme</span>(<span class="hljs-number">256</span>))<span class="hljs-keyword">val</span> tiled: <span class="hljs-type">RDD</span>[(<span class="hljs-type">SpatialKey</span>, <span class="hljs-type">MultibandTile</span>)] =
inputRdd.tileToLayout(rasterMetaData.cellType, rasterMetaData.layout, <span class="hljs-type">Bilinear</span>).repartition(<span class="hljs-number">100</span>)<span class="hljs-keyword">val</span> layoutScheme = <span class="hljs-type">ZoomedLayoutScheme</span>(<span class="hljs-type">WebMercator</span>, tileSize = <span class="hljs-number">256</span>)<span class="hljs-keyword">val</span> (zoom, reprojected): (<span class="hljs-type">Int</span>, <span class="hljs-type">RDD</span>[(<span class="hljs-type">SpatialKey</span>, <span class="hljs-type">MultibandTile</span>)] <span class="hljs-keyword">with</span> <span class="hljs-type">Metadata</span>[<span class="hljs-type">TileLayerMetadata</span>[<span class="hljs-type">SpatialKey</span>]]) =<span class="hljs-type">MultibandTileLayerRDD</span>(tiled, rasterMetaData).reproject(<span class="hljs-type">WebMercator</span>, layoutScheme, <span class="hljs-type">Bilinear</span>)<span class="hljs-keyword">val</span> attributeStore = <span class="hljs-type">FileAttributeStore</span>(outputPath)<span class="hljs-keyword">val</span> writer = <span class="hljs-type">FileCOGLayerWriter</span>(attributeStore)<span class="hljs-keyword">val</span> layerName = <span class="hljs-string">"layername"</span><span class="hljs-keyword">val</span> cogLayer =<span class="hljs-type">COGLayer</span>.fromLayerRDD(reprojected,zoom,compression = <span class="hljs-type">NoCompression</span>,maxTileSize = <span class="hljs-number">4096</span>)<span class="hljs-keyword">val</span> keyIndexes =cogLayer.metadata.zoomRangeInfos.map { <span class="hljs-keyword">case</span> (zr, bounds) =&gt; zr -&gt; <span class="hljs-type">ZCurveKeyIndexMethod</span>.createIndex(bounds) }.toMapwriter.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] = {
<span class="hljs-keyword">if</span>(minZoom.getOrElse(<span class="hljs-type">Double</span>.<span class="hljs-type">NaN</span>) != baseZoom.toDouble) {<span class="hljs-keyword">if</span>(rdd.metadata.layout.tileCols != rdd.metadata.layout.tileRows) {sys.error(<span class="hljs-string">"Cannot create Pyramided COG layer for non-square tiles."</span>)}<span class="hljs-keyword">if</span>(!isPowerOfTwo(rdd.metadata.layout.tileCols)) {sys.error(<span class="hljs-string">"Cannot create Pyramided COG layer for tile sizes that are not power-of-two."</span>)}
}<span class="hljs-keyword">val</span> layoutScheme =<span class="hljs-type">ZoomedLayoutScheme</span>(rdd.metadata.crs, rdd.metadata.layout.tileCols)<span class="hljs-keyword">if</span>(rdd.metadata.layout != layoutScheme.levelForZoom(baseZoom).layout) {sys.error(<span class="hljs-string">s"Tile Layout of layer and ZoomedLayoutScheme do not match. <span class="hljs-subst">${rdd.metadata.layout}</span> != <span class="hljs-subst">${layoutScheme.levelForZoom(baseZoom).layout}</span>"</span>)
}<span class="hljs-keyword">val</span> keyBounds =rdd.metadata.bounds <span class="hljs-keyword">match</span> {<span class="hljs-keyword">case</span> kb: <span class="hljs-type">KeyBounds</span>[<span class="hljs-type">K</span>] =&gt; kb<span class="hljs-keyword">case</span> _ =&gt; sys.error(<span class="hljs-string">s"Cannot create COGLayer with empty Bounds"</span>)}<span class="hljs-keyword">val</span> cogLayerMetadata: <span class="hljs-type">COGLayerMetadata</span>[<span class="hljs-type">K</span>] =<span class="hljs-type">COGLayerMetadata</span>(rdd.metadata.cellType,rdd.metadata.extent,rdd.metadata.crs,keyBounds,layoutScheme,baseZoom,minZoom.getOrElse(<span class="hljs-number">0</span>),maxTileSize)<span class="hljs-keyword">val</span> layers: <span class="hljs-type">Map</span>[<span class="hljs-type">ZoomRange</span>, <span class="hljs-type">RDD</span>[(<span class="hljs-type">K</span>, <span class="hljs-type">GeoTiff</span>[<span class="hljs-type">V</span>])]] =cogLayerMetadata.zoomRanges.sorted(<span class="hljs-type">Ordering</span>[<span class="hljs-type">ZoomRange</span>].reverse).foldLeft(<span class="hljs-type">List</span>[(<span class="hljs-type">ZoomRange</span>, <span class="hljs-type">RDD</span>[(<span class="hljs-type">K</span>, <span class="hljs-type">GeoTiff</span>[<span class="hljs-type">V</span>])])]()) { <span class="hljs-keyword">case</span> (acc, range) =&gt;<span class="hljs-keyword">if</span>(acc.isEmpty) {<span class="hljs-type">List</span>(range -&gt; generateGeoTiffRDD(rdd, range, layoutScheme, cogLayerMetadata.cellType, compression))} <span class="hljs-keyword">else</span> {<span class="hljs-keyword">val</span> previousLayer: <span class="hljs-type">RDD</span>[(<span class="hljs-type">K</span>, <span class="hljs-type">V</span>)] = acc.head._2.mapValues { tiff =&gt;<span class="hljs-keyword">if</span>(tiff.overviews.nonEmpty) tiff.overviews.last.tile<span class="hljs-keyword">else</span> tiff.tile}<span class="hljs-keyword">val</span> tmd: <span class="hljs-type">TileLayerMetadata</span>[<span class="hljs-type">K</span>] = cogLayerMetadata.tileLayerMetadata(range.maxZoom + <span class="hljs-number">1</span>)<span class="hljs-keyword">val</span> upsampledPreviousLayer =<span class="hljs-type">Pyramid</span>.up(<span class="hljs-type">ContextRDD</span>(previousLayer, tmd), layoutScheme, range.maxZoom + <span class="hljs-number">1</span>)._2<span class="hljs-keyword">val</span> rzz = generateGeoTiffRDD(upsampledPreviousLayer, range, layoutScheme, cogLayerMetadata.cellType, compression)(range -&gt; rzz) :: acc}}.toMap<span class="hljs-type">COGLayer</span>(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使用 COG 写入和读取相关推荐

OpenCV学习笔记（三十六）——Kalman滤波做运动目标跟踪 OpenCV学习笔记（三十七）——实用函数、系统函数、宏core OpenCV学习笔记（三十八）——显示当前FPS OpenC
OpenCV学习笔记(三十六)--Kalman滤波做运动目标跟踪 kalman滤波大家都很熟悉,其基本思想就是先不考虑输入信号和观测噪声的影响,得到状态变量和输出信号的估计值,再用输出信号的估计误差加 ...
Android项目实战（三十八）：2017最新将AndroidLibrary提交到JCenter仓库（图文教程）...
Android项目实战(三十八):2017最新将AndroidLibrary提交到JCenter仓库(图文教程) 原文:Android项目实战(三十八):2017最新将AndroidLibrary ...
JavaScript学习（三十八）—面向过程与面向对象
JavaScript学习(三十八)-面向过程与面向对象一.程序设计语言中的两大编程思想:面向对象.面向过程 (一).面向过程就是指完成某个需求的时候,先分析出完成该需求时所需要经历的步骤有哪些,然 ...
Android版疯狂填字第三关,iOS/安卓版《疯狂填字》答案攻略第三十八关
<疯狂填字>,最创新的填字玩法,挑战你的脑细胞,现在就下载.疯狂填字是最早的在线中文填字游戏,现在你可以在苹果手机上玩填字也可以在安卓手机上面玩,既打发了时间,又增长了知识,你准备好挑战了 ...
第五章第三十八题（十进制转换八进制）(Decimal to octal)
第五章第三十八题(十进制转换八进制)(Decimal to octal) **5.38(十进制转换为八进制)编写程序,提示用户输入一个十进制整数,然后显示对应的八进制值.在这个程序中不要使用Java的 ...
三十八、Fluent融化凝固模型参数设置依据
1. 融化凝固模型概述 1.1 模型原理我们在Chapter37分享了Fluent融化凝固模型案例,前文只是介绍了Fluent中的操作过程. 不知道大家会不会觉得很奇怪,Fluent模拟融化和凝固, ...
CCNA实验三十八 ZFW(区域防火墙)
CCNA实验三十八 ZFW(区域防火墙) 环境:Windows XP .Packet Tracert5.3 目的:了解ZFW的原理与基本配置说明: ZFW(Zone-Based Policy Fir ...
左耳听风第三十八周
左耳听风第三十八周每周完成一个ARTS: 每周至少做一个 leetcode 的算法题.阅读并点评至少一篇英文技术文章.学习至少一个技术技巧.分享一篇有观点和思考的技术文章.(也就是 Algorit ...
视频教程-三十八课时零基础matlab精通优化算法-Matlab
三十八课时零基础matlab精通优化算法图像和算法等领域有多年研究和项目经验:指导发表科技核心期刊经验丰富:多次指导数学建模爱好者参赛. 宋星星 ¥100.00 立即订阅扫码下载「CSDN程序员学 ...

三十八. geotrellis使用 COG 写入和读取

三十八. geotrellis使用 COG 写入和读取相关推荐

最新文章

热门文章