geotrellis使用（三）geotrellis數據處理過程分析

之前簡單介紹了geotrellis的工作過程以及一個簡單的demo，最近在此demo的基礎上實現了SRTM DEM數據的實時分析以及高程實時處理，下面我就以我實現的上述功能為例，簡單介紹一下geotrellis的數據處理過程。

一、原始數據處理

geotrellis支持geotiff的柵格數據（矢量數據還未研究），可以將geotiff直接緩存至hadoop框架下的Accumulo NOSQL數據庫，並建立金字塔等，具體處理過程在geotrellis.spark.etl.Etl類中。具體代碼如下：

def ingest[
    I: Component[?, ProjectedExtent]: TypeTag: ? => TilerKeyMethods[I, K],
    K: SpatialComponent: Boundable: TypeTag,
    V <: CellGrid: TypeTag: Stitcher: (? => TileReprojectMethods[V]): (? => CropMethods[V]): (? => TileMergeMethods[V]): (? => TilePrototypeMethods[V])
  ](
    args: Seq[String], keyIndexMethod: KeyIndexMethod[K], modules: Seq[TypedModule] = Etl.defaultModules
  )(implicit sc: SparkContext) = {
    implicit def classTagK = ClassTag(typeTag[K].mirror.runtimeClass(typeTag[K].tpe)).asInstanceOf[ClassTag[K]]
    implicit def classTagV = ClassTag(typeTag[V].mirror.runtimeClass(typeTag[V].tpe)).asInstanceOf[ClassTag[V]]

    /* parse command line arguments */
    val etl = Etl(args)
    /* load source tiles using input module specified */
    val sourceTiles = etl.load[I, V]
    /* perform the reprojection and mosaicing step to fit tiles to LayoutScheme specified */
    val (zoom, tiled) = etl.tile(sourceTiles)
    /* save and optionally pyramid the mosaiced layer */
    etl.save[K, V](LayerId(etl.conf.layerName(), zoom), tiled, keyIndexMethod)
  

重要的就是參數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上發起服務。

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

具體路由信息需要在service類中定義。service類需要繼承Actor方法，並覆蓋父類的receive方法。

override def receive = runRoute(serviceRoute)

def serviceRoute = get {
  pathPrefix("gt") {
      pathPrefix("tms")(tms) ~
      path("geoTiff")(geoTiff)
  } ~
    pathEndOrSingleSlash {
      getFromFile(staticPath + "/index.html")
    } ~
    pathPrefix("") {
      getFromDirectory(staticPath)
    }
}

以上就是建立了service的路由匹配表以及具體的控制器。當只請求IP及相應端口時會請求index.html，請求gt/tms時交給tms控制器，gt/geotiff交給geotiff控制器，其他會去匹配靜態地址，如圖片、
js、css等。

三、瓦片調用

調取數據最簡單的方式就是顯示瓦片。前端使用openlayer、leaflet均可。以leaftlet為例，在js中添加以下代碼：

WOLayer = new L.tileLayer(server + 
                    'gt/tms/{z}/{x}/{y}', {
                    format: 'image/png',
                });
WOLayer.addTo(map);

前台便會請求后台的tms控制器，tms控制器定義如下：

tms獲取到請求的x、y、z、值，並從Accumulo中取出相應的瓦片交給leaftlet，leaflet將瓦片數據放到合適的位置，便完成了瓦片的加載，從Accumulo中取出瓦片的的大致代碼如下：

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編號值。讀到瓦片之后將數據發送到前台的代碼如下：

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

其實就是調用Tile類的renderPng方法，然后將Png數據轉換成bytes發送到前端。

 

四、高級瓦片調用

當然如果只是簡單的調用瓦片，那就沒有必要非要使用geotrellis了，很多工具包括arcgis、tilemill等都包含此功能，使用geotrellis不僅是其基於Spark框架能分布式運行，而是geotrellis提供了強大的分布式計算能力，比如我們想要划定區域內的瓦片，而此區域不是標准的矩形，即不是請求完整的瓦片，這時候采用普通的框架很難完成，而采用geotrellis卻得心應手，只需要使用以下代碼即可：

val maskedTile = {
     val poly = maskz.parseGeoJson[Polygon]
     val extent: Extent = attributeStore.read[TileLayerMetadata[SpatialKey]](LayerId(LayerName, zoom), Fields.metadata).mapTransform(key)
     tile.mask(extent, poly.geom)
}

其中maskz是前端想要顯示內容的區域（Polygon），attributeStore是AccumuloAttributeStore對象，同樣可以看出是一個操作Accumulo的對象，attributeStore主要完成的功能就是讀取當前瓦片的extent即外接矩形范圍。通過調用Tile類的mask方法將請求的polygon與extent做交集，只取相交的部分的數據，再將此數據發到前端，在前端便能看到只顯示設定區域內瓦片的效果。

 

五、統計分析

如果只是進行區域內瓦片顯示，明顯意義也不大（哈哈，王婆賣瓜），geotrellis還能完成各種復雜的基於數據的統計分析（只有你想不到的，沒有你做不到的）。比如我現在做的一個demo就是統計分析給定區域內（Polygon）的高程信息（包含最大值、最小值、平均值）。

首先將DEM數據使用Etl.ingest方法導入Accumulo，注意此時就可以將--layoutScheme設置為floating，這樣就不需要建立金字塔，只取第0層數據即可，即節省存儲空間、切割時間又保證數據的一致性。

val layerId = LayerId(layer, 0)
val raster = reader.read[SpatialKey, Tile, TileLayerMetadata[SpatialKey]](layerId)
val masked = raster.mask(polygon)
val mapTransform = masked.metadata.mapTransform
val maps = masked map { case (k: SpatialKey, tile: Tile) =>
    val extent: Extent = mapTransform(k)
    val hist: Histogram[Int] = tile.polygonalHistogram(extent, extent.toPolygon())

    var max, min = hist.maxValue().getOrElse(0)
    var count:Long = 0
    var sum : Double = 0
    hist.foreach((s1:Int, s2:Long) => {

        if (max < s1) max = s1
        if (min > s1) min = s1
        sum += s1 * s2
        count += s2
    })
    (max, min, sum, count)
}
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) }
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數據處理過程分析