创建矢量数据

矢量数据集创建也并没有非常难,就是把读取的时候的顺序重复以便就好了。总体来说这个过程就是构建数据源->构建层->构建要素->构建形状->关闭数据源。下面就通过一个例子来说明如何创建矢量数据集。

   1 import ogr,os,math
   2 ds = ogr.Open("e:/gisdata/shp/lines.shp")
   3 driver = ogr.GetDriverByName("ESRI Shapefile")
   4 outputfile = 'linescopy.shp'
   5 if os.access( outputfile, os.F_OK ):
   6     driver.DeleteDataSource( outputfile )
   7 newds = driver.CreateDataSource(outputfile)
   8 layernew = newds.CreateLayer('linescopy',None,ogr.wkbLineString)
   9 layer = ds.GetLayer()
  10 extent = layer.GetExtent()
  11 print extent
  12 feature = layer.GetNextFeature()
  13 while feature is not None:
  14     layernew.CreateFeature(feature)
  15     feature = layer.GetNextFeature()
  16 newds.Destroy()
  17 #########################################################################
  18 extfile = 'rect.shp'
  19 if os.access( extfile, os.F_OK ):
  20     driver.DeleteDataSource( extfile )
  21 newds = driver.CreateDataSource(extfile)
  22 layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
  23 width = math.fabs(extent[1]-extent[0])
  24 height = math.fabs(extent[3]-extent[2])
  25 tw = width/2
  26 th = width/2
  27 extnew = extent[0]+tw
  28 wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
  29         extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
  30         extent[0],extent[3]-th, extent[0],extent[3])
  31 print wkt
  32 geom = ogr.CreateGeometryFromWkt(wkt)
  33 feat = ogr.Feature(layernew.GetLayerDefn())
  34 feat.SetGeometry(geom)
  35 layernew.CreateFeature(feat)
  36 newds.Destroy()
  37 ###########################################################################
  38 outputfile = 'srf.shp'
  39 if os.access( outputfile, os.F_OK ):
  40     driver.DeleteDataSource( outputfile )
  41 newds = driver.CreateDataSource(outputfile)
  42 layernew = newds.CreateLayer('srf',None,ogr.wkbLineString)
  43 layer.SetSpatialFilter(geom)
  44 feature = layer.GetNextFeature()
  45 while feature is not None:
  46     layernew.CreateFeature(feature)
  47     feature = layer.GetNextFeature()
  48 newds.Destroy()

这段代码其实有三段,干了3件不同的事。这些事除了构建数据源外,还实践了SetSpatialFilter这个方法的使用,以及拷贝数据源。这三件事分别被两行######号隔开了。

第一段代码,就是进行数据拷贝,把一个叫lines.shp的数据拷贝成名叫linescopy.shp的数据。可以看到,拷贝一个数据源是很容易的(我这里已经把它做得复杂了。实际上直接在DataSource和layer层上都可以进行拷贝操作了,之所以做得这么复杂,只是为了说明如何层层拨开矢量datasource神秘的外衣)。第三行首先创建一个Shapefile的驱动,然后用这个驱动创建一个数据源(第7行),再用数据源创建一个layer(第8行),再用layer创建Feature(第12~15行)。最后关掉创建好的数据集(第16行)。

特别注意:16行的Destroy不能省略,Destroy除了销毁数据还担负了数据flush到磁盘的任务。如果没有这一句,那么刚才创建的一系列Feature都不会写入磁盘的文件中。

第二段的代码解释了如何单独创建一个矢量数据集,这里我创建了一个矩形,用的是wkt创建,当然对于人工创建Geometry,wkt是最直观的,因为字符串是人为可以外部控制的。

代码看起来和第一段也很像,也是创建Driver,创建dataSource,创建Layer,接下来就创建Feature了,但是feature需要geometry。于是我们就用wkt创建了一个矩形。这个矩形是在上面的lines.shp数据集的左上角的一半位置(1/4的总面积),把矩形四个角点写成wkt,通过CreateGeometryFromWkt创建了一个矩形形状,然后放入Feature中,最后用Destroy写入磁盘。就构成了一个矩形的数据集。

第三段代码如果看懂了上面的两端代码,就很容易理解,另外,代码附带解释了SetSpatialFilter的用法和效果(第43行)。这里可以看到我们用第二段代码绘制成的矩形作为空间过滤的窗口去过滤第一段代码拷贝的数据。然后把过滤后的数据源输出成一个新的数据源。以验证SetSpatialFilter的结果。这个结果如下:

attachment:setspatialfilter.jpg

sr图层(红色的线)是第三段代码生成的,rect(蓝色矩形)是第二段代码生成的,linescopy(绿色的线,一部分被sr图层覆盖)是原始数据。可以看出效果了。ogr的文档中解释SetSpatialFilter的过滤条件是源数据的feature的外包矩形和过滤形状的外包矩形互相有重叠的就被选入,而没有重叠的就被排除(而且是整根同时被选)。看看结果,果然如此,左数第五根线虽然和矩形没有搭界,但是外包和矩形相交,所以也被选中了,输出到sr图层中。左下角那根线其实和左数第2根线其实是一根,所以也被选入了。

当然,创建一个矢量数据源不是那么容易的,除了Feature,还有属性数据。