geotrellis使用（四）geotrellis数据处理部分细节

2021-07-03 23:59:41

前面写了几篇博客介绍了Geotrellis的简单使用，具体链接在文后，今天我主要介绍一下Geotrellis在数据处理的过程中需要注意的细节，或者一些简单的经验技巧以供参考。

一、直接操作本地Geotiff

   如果不想将tiff数据切割成瓦片存放到集群中，也可以直接使用Geotrellis操作本地geotiff文件，可以直接使用SinglebandGeoTiff读取单波段的tiff，使用MultibandGeoTiff读取多波段tiff。

val geotiff = SinglebandGeoTiff("data/test.tif")

   然后使用geotiff.tile就可以像处理普通瓦片那样处理整幅tiff图像。

二、Geotiff数据处理需要注意的细节

如果需要将geotiff数据切割并上传到集群首先需要处理的是geotiff的数据类型、无数据值等元数据信息，即前期处理数据的时候需要将tiff文件处理到合适的情况以方便在程序中使用。

与数据类型和无数据值相关的属性是Tile类的CellType，Geotrellis中定义了与各种类型相对应的CellType类型，具体在geotrellis.raster.CellType类中，当然程序中可以使用tile.convert(ShortConstantNoDataCellType)将瓦片的类型转换为其他形式。

三、获取瓦片编号或者瓦片的范围（Extent）

  将数据上传到集群后，一般可以使用LayerReader将整层的瓦片信息全部读出，

val r = reader.read[SpatialKey, Tile, TileLayerMetadata[SpatialKey]](layerId)

   其中reader为LayerRander的实例，可以是AccumuloLayerReader等，具体看用户将瓦片数据存放到什么位置，layerId是存放信息的实例，包含存放的layer名称以及第几层，然后就可以使用r.metadata.mapTransform函数获取

瓦片的范围或者瓦片的编号，如果该函数的参数是一个key（瓦片编号实例），结果就是瓦片的Extent，如果参数是一个point，算出来的就是包含该点的瓦片的key。

四、数据的重投影

程序中如果需要对tile进行点、线、面的相交取值等处理就必须使用与tile相同的投影方式，否则处理过程中会出现错误，可以使用ReProject首先对点、线、面进行重投影。

Reproject(geo, LatLng, WebMercator)

其中geo表示需要进行重投影的点、线、面，LatLng是原始投影方式，WebMercator是需要转换到的投影方式。Geotrellis中定义了一个CRS类用于记录投影信息。LatLng和WebMercator继承了CRS类，是定义好的4326和3857投影方式，其他投影类型可以使用CRS类中提供的fromEpsgCode等方法进行设置。

五、获取某个坐标点对应的值

如果我们想要获取某个坐标点所对应的数据的值，有两种方式，第一种是使用LayerReader先读取整层瓦片信息，然后根据偏移得到改点的值，具体方法如下：

    val r = reader.read[SpatialKey, Tile, TileLayerMetadata[SpatialKey]](layerId)

    val mapTransform = r.metadata.mapTransform

    val key = r.metadata.mapTransform(point)

    val dataValues: Seq[Double] = r.asRasters().lookup(key).map(_.getDoubleValueAtPoint(point))

    val value = if(dataValues == null || dataValues.length <= 0) 0 else dataValues.head

其中reader和layerId的意义与上文相同，同样key就是根据坐标点的偏移算出的瓦片坐标，然后在所有瓦片中查找此key并且获取该坐标点的值，若多个瓦片均包含该坐标点会获取多个值，取出第一个。

第二种方式是使用ValueReader直接找到包含改点的瓦片，然后根据偏移得到此点的数据，具体代码如下：

    val key = attributeStore.readMetadata[TileLayerMetadata[SpatialKey]](layerId).mapTransform(point)

    val (col, row) = attributeStore.readMetadata[TileLayerMetadata[SpatialKey]](layerId).toRasterExtent().mapToGrid(point)

    val tile: Tile =  tileReader.reader[SpatialKey, Tile](layerId).read(key)

    val tileCol = col - key.col * tile.cols

    val tileRow = row - key.row * tile.rows

    println(s"tileCol=${tileCol}   tileRow = ${tileRow}")

    tile.get(tileCol, tileRow)

其中attributeStore是元数据信息，与用户数据上传的位置有关，key是从元数据中根据坐标点偏移算出的瓦片编号，tilReader是ValueReader实例，col、row是根据偏移算出的坐标点在整个数据集中的像素偏移值，然后通过减去编号乘以瓦片像素个数来获取相对当前瓦片的偏移，从而实现获取当前坐标点的数据值。

两种方式均能得到坐标点对应的值，但是其效率却相差几十倍，在我自己的测试中，使用ValueReader取到数据值大概需要20ms，而使用layerReader则大概需要6000ms，我猜测应当是使用LayerReader的方式会在所有瓦片中lookup，而ValueReader则直接获取单个瓦片，所以效率存在差别。

六、结束语

本文简单记录了近期使用Geotrellis过程中遇到的一些问题，及其解决方案，目前项目只用到了栅格数据，所以只是针对Raster模块，后续会探索其他模块功能，并随时将心得发布到博客园中，欢迎大家共同探讨。

七、参考文献

一、geotrellis使用初探
二、geotrellis使用（二）geotrellis-chatta-demo以及geotrellis框架数据读取方式初探
三、geotrellis使用（三）geotrellis数据处理过程分析四、geotrellis使用（四）geotrellis数据处理部分细节

码农公寓