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