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使用（八）矢量数据栅格化