geotrellis使用（八）矢量數據柵格化

目錄

1. 前言
2. 柵格化處理
3. 總結
4. 參考鏈接

一、前言

       首先前幾天學習了一下Markdown，今天將博客園的編輯器改為Markdown，從編寫博客到界面美觀明顯都清爽多了，也能寫出各種樣式的東西了，有關Markdown，網上內容很多，暫且不表，開始進入今天的主題。

       前幾天碰到一個任務，需要將矢量數據導入到Accumulo中，然后通過geotrellis進行調用。這一下又犯難了，之前處理的全是raster數據，通過ETL類可以直接進行導入生成金字塔等，如何將矢量數據導入平台之前未曾碰到，但是大致分析首先需要進行柵格化，因為柵格化之后就可以直接使用Geotrellis進行處理，矢量數據柵格化之前也未遇到過，解決問題就要一步步來，一步步分析，下面就為大家講解我本次實現的過程。

二、柵格化處理

       要想柵格化第一步肯定需要讀取矢量數據。

讀取矢量數據

本文中主要講解shapefile，數據庫部分后面講解。

       首先瀏覽Geotrellis的源代碼，發現一個ShapeFileReader類，貌似直接能解決問題啊，趕緊寫代碼如下:

geotrellis.shapefile.ShapeFileReader.readSimpleFeatures(path)

       滿心歡喜的以為一句話就解決問題了，誰知道一直報如下錯誤：

The following locker still has a lock: read on file:..shp by org.geotools.data.shapefile.shp.ShapefileReader 
The following locker still has a lock: read on file:..shx by org.geotools.data.shapefile.shp.IndexFile 
The following locker still has a lock: read on file:...dbf by org.geotools.data.shapefile.dbf.DbaseFileReader 
Exception in thread "main" java.lang.IllegalArgumentException: Expected requestor org.geotools.data.shapefile.dbf.DbaseFileReader@4ea5b703 to have locked the url but it does not hold the lock for the URL

       實驗了各種方法無果，那么看一下他的源代碼，然后直接拿過來用，發現可以，代碼如下：

/**
    * get the features from shape file by the attrName,default "the_geom"
    * @param path
    * @return mutable.ListBuffer[Geometry]
    */
  def getFeatures(path: String, attrName: String = "the_geom", charset: String = "UTF-8"): mutable.ListBuffer[Geometry] ={
    val features = mutable.ListBuffer[Geometry]()
    var polygon: Option[MultiPolygon] = null
    val shpDataStore = new ShapefileDataStore(new File(path).toURI().toURL())
    shpDataStore.setCharset(Charset.forName(charset))
    val typeName = shpDataStore.getTypeNames()(0)
    val featureSource = shpDataStore.getFeatureSource(typeName)
    val result = featureSource.getFeatures()
    val itertor = result.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(attrName)) {
          features += WKT.read(pro.getValue.toString) //get all geom from shp
        }
      }
    }
    itertor.close()
    shpDataStore.dispose()
    features
  }

       實驗中的shape文件包含一個字段the_geom，里面存儲了空間信息的WKT語句，所以程序中讀出該屬性的值然后使用WKT.read(pro.getValue.toString)將其轉換成Geometry對象。

注意最后需要添加shpDataStore.dispose()否則會同樣報上述文件鎖定的錯誤，所以我猜測此處應該是Geotrellis的一個bug。

通過上述可以得出其實通過數據庫讀取矢量數據也只是個驅動的問題，只要將需要的記錄逐行讀出然后轉化為Geometry對象即可，后面會通過一篇博客詳細說明。

       讀出了矢量數據后，緊接着就是將數據映射到柵格圖像上。

將Geometry數組對象進行柵格化

獲取Geometry數組對象的空間范圍RasterExtent

       柵格化后的數據仍然包含了投影、空間范圍等空間信息以及分辨率、圖像尺寸等柵格信息，所以我們要先根據Geometry數組求出這些信息。

• 獲取經緯度范圍

       一個簡單的循環遍歷所有要素比較最大最小值的方法，代碼如下：

var minX = features(0).jtsGeom.getEnvelopeInternal.getMinX
var minY = features(0).jtsGeom.getEnvelopeInternal.getMinY
var maxX = features(0).jtsGeom.getEnvelopeInternal.getMaxX
var maxY = features(0).jtsGeom.getEnvelopeInternal.getMaxY
for (feature <- features) {
if (feature.jtsGeom.getEnvelopeInternal.getMaxX > maxX)
    maxX = feature.jtsGeom.getEnvelopeInternal.getMaxX
if (feature.jtsGeom.getEnvelopeInternal.getMaxY > maxY)
    maxY = feature.jtsGeom.getEnvelopeInternal.getMaxY
if (feature.jtsGeom.getEnvelopeInternal.getMinX < minX)
    minX = feature.jtsGeom.getEnvelopeInternal.getMinX
if (feature.jtsGeom.getEnvelopeInternal.getMinY < minY)
    minY = feature.jtsGeom.getEnvelopeInternal.getMinY
}

• 計算柵格化后的圖像尺寸

       柵格圖像包含分辨率、像素大小、cols、row等要素，在這里我簡單的理解為可以根據矢量數據的經緯度范圍差除以分辨率來得到cols、rows，通過查閱資料可以發現當zoom（表示瓦片的層級）為22時，分辨率為0.037323，所以這里可以簡單的算出其他層級的分辨率如下：

val resolution = 0.037323 * Math.pow(2, 22 - zoom)

       得到了分辨率后即可用范圍差除以分辨率得到圖像尺寸。

此處需要注意圖像的空間參考，若參考不同時需要進行投影轉換：val res1 = Reproject((minX, minY), LatLng, WebMercator)

• 得到RasterExtent

RasterExtent(new Extent(minX, minY, maxX, maxY), cols, rows)

柵格化

       經過查閱Geotrellis的源代碼以及咨詢官方大牛，大概明白了可以使用Rasterizer類進行柵格化操作，其實也很簡單，只需要一句代碼如下：

Rasterizer.rasterizeWithValue(features, re, 100)

       其中features即從shp文件中讀出的Geometry數組，re為上文中得到的RasterExtent，100表示將這些對象在柵格中賦予的像素值。

       柵格化效果如下：

       矢量數據

       柵格化數據

三、總結

       通過以上代碼便完成了柵格化操作，看似沒幾行代碼，確確實實也折騰了很久，主要是對Geotrellis的源代碼還不夠熟悉，對一些基礎的地理空間信息知識掌握還不夠到位。

四、參考鏈接

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