geotrellis使用（三）geotrellis数据处理过程分析

之前简单介绍了geotrellis的工作过程以及一个简单的demo，最近在此demo的基础上实现了SRTM DEM数据的实时分析以及高程实时处理，下面我就以我实现的上述功能为例，简单介绍一下geotrellis的数据处理过程。

一、原始数据处理

geotrellis支持geotiff的栅格数据（矢量数据还未研究），可以将geotiff直接缓存至hadoop框架下的Accumulo NOSQL数据库，并建立金字塔等，具体处理过程在geotrellis.spark.etl.Etl类中。具体代码如下：

 1 def ingest[
 2     I: Component[?, ProjectedExtent]: TypeTag: ? => TilerKeyMethods[I, K],
 3     K: SpatialComponent: Boundable: TypeTag,
 4     V <: CellGrid: TypeTag: Stitcher: (? => TileReprojectMethods[V]): (? => CropMethods[V]): (? => TileMergeMethods[V]): (? => TilePrototypeMethods[V])
 5   ](
 6     args: Seq[String], keyIndexMethod: KeyIndexMethod[K], modules: Seq[TypedModule] = Etl.defaultModules
 7   )(implicit sc: SparkContext) = {
 8     implicit def classTagK = ClassTag(typeTag[K].mirror.runtimeClass(typeTag[K].tpe)).asInstanceOf[ClassTag[K]]
 9     implicit def classTagV = ClassTag(typeTag[V].mirror.runtimeClass(typeTag[V].tpe)).asInstanceOf[ClassTag[V]]
10 
11     /* parse command line arguments */
12     val etl = Etl(args)
13     /* load source tiles using input module specified */
14     val sourceTiles = etl.load[I, V]
15     /* perform the reprojection and mosaicing step to fit tiles to LayoutScheme specified */
16     val (zoom, tiled) = etl.tile(sourceTiles)
17     /* save and optionally pyramid the mosaiced layer */
18     etl.save[K, V](LayerId(etl.conf.layerName(), zoom), tiled, keyIndexMethod)
19   
重要的就是参数args，geotrellis根据不同的参数将数据进行不同的处理。具体的参数信息在https://github.com/geotrellis/geotrellis/blob/master/docs/spark-etl/spark-etl-intro.md
中均有介绍，这里介绍一些重要的配置。
1、--layoutScheme      layoutScheme有tms和floating两种选项，如果用floating切瓦片的时候只有0层，切记这一点，因为调用瓦片的时候跟层有很大关系；用tms会建立金字塔。相当于用floating处理的就是原始数据只将数据切割成256*256的块，层为0（具体x、y编号不需要操心，geotrellis会自动计算），用tms会将数据从最大层（此最大层根据数据的分辨率计算得出）切到第一层，调用的时候直接根据层进行调用。
2、--pyramid              加上此参数在layoutScheme=tms的时候系统会建立金字塔
3、-I path=file:/..        如果此处的路径为文件，则单独导入此文件，如果为文件夹，则一次将整个路径导入，并且会自动拼接，瓦片不会有缝隙，这一点非常漂亮，此处只能用漂亮来形容，geotrellis不但能够分布式瓦片切割，还能自动拼接，实在是漂亮。
4、--layer                   此参数用于区分不同的数据，取数据的时候根据此项区分不同的数据。
通过简单的调用ingest方法就能进行分布式瓦片切割，不得不说geotrllis提供了很多强大的功能。
 

二、发起服务
要对外提供数据，系统首先要能够发起服务，geotrellis建立一个服务也很容易，只需要使用以下语句系统遍自动的在host和相应的port上发起服务。

1 IO(Http) ! Http.Bind(service, host, port)

具体路由信息需要在service类中定义。service类需要继承Actor方法，并覆盖父类的receive方法。

 1 override def receive = runRoute(serviceRoute)
 2 
 3 def serviceRoute = get {
 4   pathPrefix("gt") {
 5       pathPrefix("tms")(tms) ~
 6       path("geoTiff")(geoTiff)
 7   } ~
 8     pathEndOrSingleSlash {
 9       getFromFile(staticPath + "/index.html")
10     } ~
11     pathPrefix("") {
12       getFromDirectory(staticPath)
13     }
14 }
以上就是建立了service的路由匹配表以及具体的控制器。当只请求IP及相应端口时会请求index.html，请求gt/tms时交给tms控制器，gt/geotiff交给geotiff控制器，其他会去匹配静态地址，如图片、
js、css等。

三、瓦片调用

调取数据最简单的方式就是显示瓦片。前端使用openlayer、leaflet均可。以leaftlet为例，在js中添加以下代码：
1 WOLayer = new L.tileLayer(server + 
2                     'gt/tms/{z}/{x}/{y}', {
3                     format: 'image/png',
4                 });
5 WOLayer.addTo(map);

前台便会请求后台的tms控制器，tms控制器定义如下：

tms获取到请求的x、y、z、值，并从Accumulo中取出相应的瓦片交给leaftlet，leaflet将瓦片数据放到合适的位置，便完成了瓦片的加载，从Accumulo中取出瓦片的的大致代码如下：

1 val tile: Tile = tileReader.reader[SpatialKey, Tile](LayerId(LayerName, zoom)).read(key)

其中tileReader是一个AccumuloValueReader对象，很明显看出此对象是一个有关Accumulo的对象，其中包含Accumulo的用户密码等。LayerName就是上文中导入数据时候设置的layer参数对应的值。key是个SpatialKey对象，val key = SpatialKey(x, y)，记录了瓦片x、y编号值。读到瓦片之后将数据发送到前台的代码如下：

1 respondWithMediaType(MediaTypes.`image/png`) {
2         val result = tile.renderPng().bytes
3         complete(result)
4 }

其实就是调用Tile类的renderPng方法，然后将Png数据转换成bytes发送到前端。

 

四、高级瓦片调用

当然如果只是简单的调用瓦片，那就没有必要非要使用geotrellis了，很多工具包括arcgis、tilemill等都包含此功能，使用geotrellis不仅是其基于Spark框架能分布式运行，而是geotrellis提供了强大的分布式计算能力，比如我们想要划定区域内的瓦片，而此区域不是标准的矩形，即不是请求完整的瓦片，这时候采用普通的框架很难完成，而采用geotrellis却得心应手，只需要使用以下代码即可：
1 val maskedTile = {
2      val poly = maskz.parseGeoJson[Polygon]
3      val extent: Extent = attributeStore.read[TileLayerMetadata[SpatialKey]](LayerId(LayerName, zoom), Fields.metadata).mapTransform(key)
4      tile.mask(extent, poly.geom)
5 }

其中maskz是前端想要显示内容的区域（Polygon），attributeStore是AccumuloAttributeStore对象，同样可以看出是一个操作Accumulo的对象，attributeStore主要完成的功能就是读取当前瓦片的extent即外接矩形范围。通过调用Tile类的mask方法将请求的polygon与extent做交集，只取相交的部分的数据，再将此数据发到前端，在前端便能看到只显示设定区域内瓦片的效果。

 

五、统计分析

如果只是进行区域内瓦片显示，明显意义也不大（哈哈，王婆卖瓜），geotrellis还能完成各种复杂的基于数据的统计分析（只有你想不到的，没有你做不到的）。比如我现在做的一个demo就是统计分析给定区域内（Polygon）的高程信息（包含最大值、最小值、平均值）。
首先将DEM数据使用Etl.ingest方法导入Accumulo，注意此时就可以将--layoutScheme设置为floating，这样就不需要建立金字塔，只取第0层数据即可，即节省存储空间、切割时间又保证数据的一致性。
 1 val layerId = LayerId(layer, 0)
 2 val raster = reader.read[SpatialKey, Tile, TileLayerMetadata[SpatialKey]](layerId)
 3 val masked = raster.mask(polygon)
 4 val mapTransform = masked.metadata.mapTransform
 5 val maps = masked map { case (k: SpatialKey, tile: Tile) =>
 6     val extent: Extent = mapTransform(k)
 7     val hist: Histogram[Int] = tile.polygonalHistogram(extent, extent.toPolygon())
 8 
 9     var max, min = hist.maxValue().getOrElse(0)
10     var count:Long = 0
11     var sum : Double = 0
12     hist.foreach((s1:Int, s2:Long) => {
13
14         if (max < s1) max = s1
15         if (min > s1) min = s1
16         sum += s1 * s2
17         count += s2
18     })
19     (max, min, sum, count)
20 }
21 val (max, min, sum, count) = maps reduce { case ((z1, a1, s1, c1), (z2, a2, s2, c2)) => (Math.max(z1, z2), Math.min(a1, a2), s1 + s2, c1 + c2) }
22 val avg = sum / count

val layerId = LayerId(layer, 0)表示取的是导入数据的第0层，由于使用floating方式此处必须是0。reader是一个AccumuloLayerReader对象，此处与上面的AccumuloVlaueReader不同之处在于上文中取固定key值得瓦片，此处需要根据范围进行选择，masked就是根据polygon筛选出的结果，是一个RDD[(SpatialKey, Tile)]对象，即存储着范围内的所有瓦片以及其编号信息。对masked进行map操作，获取其单个瓦片的extent，以及polygon内的统计信息，算出最大值，最小值以及高程加权和。最后对结果进行reduce操作，获取整体的最大值、最小值、平均值。（此处平均值算法可能不妥，希望有更好建议的能够留言，感激！）。将计算到的结果发到前端，前端就能实时显示统计分析结果。

 

六、结尾

geotrellis的功能非常强大，此处只是冰山一脚，后续还会进行相关研究，经验心得会及时总结到这里，以使自己理解的更加透彻，如果能帮助到其他人也是极好的！

七、参考链接

一、geotrellis使用初探
二、geotrellis使用（二）geotrellis-chatta-demo以及geotrellis框架数据读取方式初探
三、geotrellis使用（三）geotrellis数据处理过程分析
posted @ 2016-04-22 21:38  shoufengwei  阅读(5513)  评论(5)  编辑  收藏  举报