geotrellis使用（十）緩沖區分析以及多種類型要素柵格化

目錄

1. 前言
2. 緩沖區分析
3. 多種類型要素柵格化
4. 總結
5. 參考鏈接

一、前言

       上兩篇文章介紹了如何使用Geotrellis進行矢量數據柵格化以及柵格渲染，本文主要介紹柵格化過程中常用到的緩沖區分析以及同一范圍內的多種類型要素柵格化。

       本文主要記錄今天過程中碰到的兩個問題，第一個問題就是線狀要素在進行柵格化的時候只有單個像素，看不出應有的效果；第二個問題就是同一地區的數據既包含面狀要素，又包含了線狀要素，普通方式只能柵格化成兩套數據。下面我為大家介紹解決這兩個問題的方法（當然若有人有更好的方法，歡迎交流）。

二、緩沖區分析

       緩沖區分析在百度百科中的定義為：

緩沖區分析是指以點、線、面實體為基礎，自動建立其周圍一定寬度范圍內的緩沖區多邊形圖層，然后建立該圖層與目標圖層的疊加，進行分析而得到所需結果。它是用來解決鄰近度問題的空間分析工具之一。鄰近度描述了地理空間中兩個地物距離相近的程度。

       當然本文並不是教大家如何解決鄰近度問題，只是簡單的說明如何能夠在柵格化的過程中將線狀要素能夠多外擴幾個像素。

自己實現外擴像素

       由於本人非地理信息專業出身（甚至非計算機專業出身，沒辦法，置身碼農把青春賣！）所以在遇到問題的時候並不懂什么緩沖區分析的高大上的詞匯。首先想到的是我可以在矢量化的過程中外擴幾個像素，這樣不就實現了增強的效果，但是有個問題就是我如何知道線段的方向，先將就着來，我把線段點上下左右的像素全部賦予與改點相同的值，這樣可以不用考慮方向，並且應該能達到效果。
       說干就干，再一次認真研讀Geotrellis的Reasterizer.scala的源代碼，冥思苦想一陣之后，想到了方法，主要是要重寫賦值的方法，實現代碼如下：

def rasterize(geom: Geometry, rasterExtent: RasterExtent, value: Int) ={
      val cols = rasterExtent.cols
      val array = Array.ofDim[Int](rasterExtent.cols * rasterExtent.rows).fill(NODATA)
      val f = (col: Int, row: Int) => {
        array(row * cols + col) = value
        if (col > 0)
          array(row * cols + col - 1) = value
        if (col < cols - 1)
          array(row * cols + col + 1) = value
        if (row > 0)
          array((row - 1) * cols + col) = value
        if (row < rasterExtent.rows - 1)
          array((row + 1) * cols + col) = value
      }
      Rasterizer.foreachCellByGeometry(geom, rasterExtent)(f)
      ArrayTile(array, rasterExtent.cols, rasterExtent.rows)
    }

簡單說來就是之前f函數中只有array(row * cols + col) = value一條語句，即實現當前點的像素點賦值，那么加上了判斷不是邊界之后，給上下左右的像素點都賦值即可實現，運行起來。

       得到的結果雖然看起來有點丑，但是總算解決了這個問題，然后把結果拿給老板看，老板什么話也沒說，默默的甩給我https://gitter.im/geotrellis/geotrellis/archives/2016/02/22 這么一個網址。好吧，老板果然是老板，這里也要介紹一下https://gitter.im/geotrellis/geotrellis/，這是Github中的Geotrellis項目交流群，在里面咨詢問題，會有懂的人甚至作者解答，有點考驗英語基礎。

使用buffer函數

       在那個網頁中，上來就有這么一段代碼：

val points = Seq(
  Point(re.gridToMap(100,100)).buffer(30),
  Point(re.gridToMap(200,200)).buffer(30),
  Point(re.gridToMap(300,300)).buffer(30),
  Point(re.gridToMap(400,400)).buffer(30),
  Point(re.gridToMap(500,500)).buffer(30)
)

       根據這段代碼尤其是buffer名稱，可以知道其實在Geotrellis中緩沖區分析就是使對象調用buffer函數即可，參數表示緩沖的距離。趕緊拿來試驗，非常成功，但是這里面卻有幾個需要注意的問題。

1. 緩沖距離

       此處的緩沖距離經過實際測試發現與當前數據的坐標系相一致，即如果是WGS84地理坐標系，那么此處緩沖距離就是以經緯度為單位，大地坐標系此處就是以米為單位。

2. 緩沖類型

       一般情況下只需要給點、線要素使用緩沖即可，這里就可以使用模式匹配，如下：

val geom = WKT.read(pro.getValue.toString) match {
                case geom: Point        =>    geom.buffer(bufferDistance)
                case geom: Line         =>    geom.buffer(bufferDistance)
                case geom: MultiLine    =>    geom.buffer(bufferDistance).toGeometry().get
                case geom               =>    geom
              }

       這里就僅為Point、Line以及MultiLine類型進行了緩沖區設置，其他需要轉換的可以用同樣的方式進行匹配，展示一下最終的效果。

       其實查看buffer函數的定義，不難發現該函數實現的就是將要點線要素轉換成了面要素。

       以上就實現了緩沖區分析，下面進行下一個主題多種類型要素柵格化。

三、多種類型要素柵格化

       同一個區域數據即包含面狀要素又包含線狀要素，顯然在shape文件中以及數據庫中我們都沒有辦法將其進行合並，而如果我們又不想得到兩套柵格化的數據該如何是好呢？

       其實方法也很簡單，只需要將要素拼接到同一個GeometryCollection中然后統一獲取其RasterExtent即可，實現代碼如下：

val features = mutable.ListBuffer[Geometry]()
    for (path <- paths) {
      val file = new File(path)
      if(file.exists()) {
        val shpDataStore = new ShapefileDataStore(new File(path).toURI().toURL())
        shpDataStore.setCharset(Charset.forName(charset))
        val featureSource = shpDataStore.getFeatureSource()
        val itertor = featureSource.getFeatures().features()
        while (itertor.hasNext()) {
          val feature = itertor.next()

          val p = feature.getProperties()
          val it = p.iterator()

          while (it.hasNext()) {
            val pro = it.next()
            if (pro.getName.getLocalPart.equals("the_geom")) {//get all geom from shp
              val geom = WKT.read(pro.getValue.toString) match {
                case geom: Point        =>    geom.buffer(resolution * bufferDistance)
                case geom: Line         =>    geom.buffer(resolution * bufferDistance)
                case geom: MultiLine    =>    geom.buffer(resolution * bufferDistance).toGeometry().get //0.0054932 * 7
                case geom               =>    geom
              }
              features += geom
            }
          }
        }
        itertor.close()
        shpDataStore.dispose()
      } else
        println(s"the file ${path} isn't exist")
    }

       以上代碼實現的是逐個循環需要柵格化的文件，然后將每個geometry對象添加到features中，剩下的在前面的文章中已經介紹過，不再贅述。

四、總結

       以上講述了如何進行緩沖區分析以及多種類型要素柵格化。雖然實現方法比較較難，但是在剛碰到這些問題的時候確實會讓人摸不着頭腦，本文簡單記錄之，僅為整理思路以及方便以后使用，如果能夠幫助到一些苦苦探索的人當然是更好的。最后感謝在工作過程中給予了重大幫助和指導的吳老板！

五、參考鏈接

一、geotrellis使用初探
二、geotrellis使用（二）geotrellis-chatta-demo以及geotrellis框架數據讀取方式初探
三、geotrellis使用（三）geotrellis數據處理過程分析
四、geotrellis使用（四）geotrellis數據處理部分細節
五、geotrellis使用（五）使用scala操作Accumulo
六、geotrellis使用（六）Scala並發（並行）編程
七、geotrellis使用（七）記錄一次慘痛的bug調試經歷以及求DEM坡度實踐
八、geotrellis使用（八）矢量數據柵格化
九、geotrellis使用（九）使用geotrellis進行柵格渲染
十、geotrellis使用（十）緩沖區分析以及多種類型要素柵格化