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

目錄

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

1、前言

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

       前幾天碰到一個任務，須要將矢量數據導入到Accumulo中，而後經過geotrellis進行調用。這一下又犯難了，以前處理的全是raster數據，經過ETL類能夠直接進行導入生成金字塔等，如何將矢量數據導入平臺以前不曾碰到，可是大體分析首先須要進行柵格化，由於柵格化以後就能夠直接使用Geotrellis進行處理，矢量數據柵格化以前也未遇到過，解決問題就要一步步來，一步步分析，下面就爲你們講解我本次實現的過程。

2、柵格化處理

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

讀取矢量數據

本文中主要講解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表示將這些對象在柵格中賦予的像素值。

       柵格化效果以下：

       矢量數據

       柵格化數據

3、總結

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

4、參考連接

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