<<TableOfContents>>

= 创建矢量数据 =
== 一个例子 ==
矢量数据集创建也并没有非常难，就是把读取的时候的顺序重复以便就好了。总体来说这个过程就是构建数据源->构建层->构建要素->构建形状->关闭数据源。下面就通过一个例子来说明如何创建矢量数据集。

{{{
#!python
import ogr,os,math
ds = ogr.Open("e:/gisdata/shp/lines.shp")
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = 'linescopy.shp'
if os.access( outputfile, os.F_OK ):
    driver.DeleteDataSource( outputfile )
newds = driver.CreateDataSource(outputfile)
layernew = newds.CreateLayer('linescopy',None,ogr.wkbLineString)
layer = ds.GetLayer()
extent = layer.GetExtent()
print extent
feature = layer.GetNextFeature()
while feature is not None:
    layernew.CreateFeature(feature)
    feature = layer.GetNextFeature()
newds.Destroy()
#########################################################################
extfile = 'rect.shp'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
width = math.fabs(extent[1]-extent[0])
height = math.fabs(extent[3]-extent[2])
tw = width/2
th = width/2
extnew = extent[0]+tw
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
        extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
        extent[0],extent[3]-th, extent[0],extent[3])
print wkt
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
newds.Destroy()
###########################################################################
outputfile = 'srf.shp'
if os.access( outputfile, os.F_OK ):
    driver.DeleteDataSource( outputfile )
newds = driver.CreateDataSource(outputfile)
layernew = newds.CreateLayer('srf',None,ogr.wkbLineString)
layer.SetSpatialFilter(geom)
feature = layer.GetNextFeature()
while feature is not None:
    layernew.CreateFeature(feature)
    feature = layer.GetNextFeature()
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，还有属性数据。

== 属性数据 ==
我们可以给那个矩形添加属性数据。

{{{
#!python
# -*- encoding:gb2312 -*-
import ogr,os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'rect2.shp'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect2',None,ogr.wkbPolygon)
fieldcnstr = ogr.FieldDefn("关于",ogr.OFTString)
fieldcnstr.SetWidth(32)
layernew.CreateField(fieldcnstr)
fieldf = ogr.FieldDefn("f",ogr.OFTReal)
layernew.CreateField(fieldf)
#laydef = layernew.GetLayerDefn()
#fielde = ogr.FieldDefn('elev',ogr.OFTReal)
#laydef.AddFieldDefn(fielde)
extent = [0,0,300,300]
tw = 100
th = 100
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
extent[0],extent[3]-th, extent[0],extent[3])
print wkt
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
feat.SetField('关于',"啊呀！这孩子呵！您瞧！那么……。阿唷！哈哈！Heh e！")# 测试中文字段名和字段值
feat.SetField('f',127.546)
layernew.CreateFeature(feat)
newds.Destroy()
}}}
这个例子可以看到如何创建一个带属性的数据。要数据带属性表，就要首先定义表头，定义表头就要定义数个字段，而要定义一个字段就要定义一个字段描述描述（第11～12行，第16行），然后把这个字段描述创建到layer中（第14，18行）。有了表头，就可以在输入Feature时添加属性表的内容。方法就用!SetField（第35～36行）。 然后用一个可以开shapefile的软件打开看看，就可以看到结果了。

结果基本符合我们的想像，但是中文不正常。只到了句号，后面截掉了。哈哈，这个就不是ogr的问题了。因为我们设置字符串宽度（长度）是32，就只能显示16个汉字（一个gb2312编码汉字两个字节）。要不信，你把32改大一点，后面的就出来了。

第20～22行有一个添加字段的方法，但是可以在属性表中添加字段，但是这个时刻放在这个位置，python运行就会出错。察看ogr的文档，发现下面的解释：

This method should only be called while there are no[[http://ogr.maptools.org/classOGRFeature.html|OGRFeature]]objects in existance based on this OGRFeatureDefn. The[[http://ogr.maptools.org/classOGRFieldDefn.html|OGRFieldDefn]]passed in is copied, and remains the responsibility of the caller.

这个方法只能在没有任何基于这个OGRFeatureDefn的OGRFeature对象存在的时候才可以调用。传入的OGRFieldDefn会被拷贝，并且保留调用者的责任。

 . 如果剪切到38行处运行，就会发现这回没错，但是在数据中不会有任何效果。下面的例子似乎也没有任何效果。
{{{
#! python
extfile = 'rect3.shp'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect3',None,ogr.wkbPolygon)
fieldcnstr = ogr.FieldDefn("关于",ogr.OFTString)
fieldcnstr.SetWidth(36)
fieldf = ogr.FieldDefn("f",ogr.OFTReal)
laydef = layernew.GetLayerDefn()
laydef.AddFieldDefn(fieldcnstr)
laydef.AddFieldDefn(fieldf)
#print laydef.GetFieldCount()
newds.Destroy()
}}}
输出的图层的属性表中空荡荡没有任何东西。是否我们可以这么认为－－这个函数只是在内存中虚拟得建立了一个字段，而不与磁盘发生关系？或者说这个函数只是用来建立一个属性字段的框架，仅仅是用来定义，不和原来的数据发生联系，而在建立新的数据时用它来进行参考和引用？

== 恶搞 ==
我一直对ogr定义的这层layer的概念理解地很模糊。好吧，我们就对它进行一下恶搞吧，通过恶搞，看看能搞出一些什么惊人的东西。

=== 把layer的名称定义得和datasource不一样 ===
{{{
#! python
# -*- encoding:gb2312 -*-
import ogr,os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'rect2.shp'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rectkkk',None,ogr.wkbPolygon)
extent = [0,0,300,300]
tw = 100
th = 100
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
extent[0],extent[3]-th, extent[0],extent[3])
print wkt
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
feat.SetField('关于',"啊呀！这孩子呵！您瞧！那么……。阿唷！哈哈！Heh e！")
feat.SetField('f',127.546)
layernew.CreateFeature(feat)
newds.Destroy()
ds = ogr.Open("rect2.shp")
print ds.GetLayerCount()
layer = ds.GetLayer(0)
print layer.GetName()
layer = ds.GetLayerByName('rectkkk')
}}}
运行，发现生成后的图层名字还叫rect2，和rectkkk没有什么关系（运行到最后一步出错，说没有图层叫做rectkkk）。 另外我们还发现!SetField在字段不存在的时候不会抛出错误。

=== 生成多层Layer ===
我们都知道，shapefile都是单层的，而且层中都是同种feature。如果我们生成多层layer，每个layer的feature种类不同会怎么样？

{{{
#! python
# -*- encoding:gb2312 -*-
import ogr,os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'rect2.shp'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rectkkk',None,ogr.wkbPolygon)
extent = [0,0,300,300]
tw = 100
th = 100
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
extent[0],extent[3]-th, extent[0],extent[3])
print wkt
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
layernew2 = newds.CreateLayer('rectaaa',None,ogr.wkbPoint)
wktp = 'POINT (%f %f)' % (70.4,220.34)
geom2 = ogr.CreateGeometryFromWkt(wktp)
feat2 = ogr.Feature(layernew2.GetLayerDefn())
feat2.SetGeometry(geom2)
layernew2.CreateFeature(feat2)
newds.Destroy()
}}}
奇迹出现了！！！我们生成了两套shapefile！！一套叫rect2，一套叫rectaaa。就是说我们生成了两层图层。第一套和datasource同名，一套和第二个layer同名。一个是多边形图层一个是点图层！这个结果真……真幽默。 一个偶然的情况下，我把extfile = 'rect2.shp'写成了rec2，没有shp扩展名，终于出现了一个让我舒心的结果。生成了一个rect2的文件夹，然后在文件夹下生成了一套叫rectkkk的shapefile以及叫rectaaa的shapefile。一切都清楚了……

这让我们对!DataSource有了更深入的认识，不但一个文件可以作为!DataSource，一个目录也可以。而目录下的所有shapefile就是一个个的Layer。

换个驱动看看吧，换个mapinfo的怎么样？要知道，一个shapefile中只能放同种类型的feature，而一个tab里面却可以存放很多不同类型的feature啊。只要一个tab在可编辑状态下，加入点，线，面，后保存，最后只形成单一的tab文件。一个图层只存储一种同类的feature，那么是不是意味着一个mapinfo的文件数据源中可以有许多的Layer，而一个layer只存放同种类型的Feature呢？

{{{
#! python
# -*- encoding:gb2312 -*-
import ogr,os,math
ogr.GetDriver
driver = ogr.GetDriverByName("MapInfo File")
extfile = 'mi.tab'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)#,options=['FORMAT=MIF'])
layernew = newds.CreateLayer('mitab',None,ogr.wkbPolygon)
fieldid = ogr.FieldDefn('FID',ogr.OFTInteger)
layernew.CreateField(fieldid)
extent = [0,0,300,300]
tw = 100
th = 100
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
extent[0],extent[3]-th, extent[0],extent[3])
print wkt
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
layernew2 = newds.CreateLayer('mitabaaa',None,ogr.wkbPoint)
fieldid2 = ogr.FieldDefn('FID',ogr.OFTInteger)
layernew2.CreateField(fieldid2)
wktp = 'POINT (%f %f)' % (70.4,220.34)
geom2 = ogr.CreateGeometryFromWkt(wktp)
feat2 = ogr.Feature(layernew2.GetLayerDefn())
feat2.SetGeometry(geom2)
layernew2.CreateFeature(feat2)
newds.Destroy()
}}}
出错了？

{{{
#! python
C:\WINDOWS\system32\cmd.exe /c D:\Python24\python.exe create.py
POLYGON ((0.000000 300.000000,100.000000 300.000000,100.000000 200.000000,0.0000
00 200.000000,0.000000 300.000000))
ERROR 1: Unable to create new layers in this single file dataset.
Traceback (most recent call last):
File "create.py", line 24, in ?
layernew2 = newds.CreateLayer('mitabaaa',None,ogr.wkbPoint)
File "D:\Python24\lib\site-packages\ogr.py", line 532, in CreateLayer
raise OGRError, gdal.GetLastErrorMsg()
ogr.OGRError: Unable to create new layers in this single file dataset.
shell returned 1
Hit any key to close this window...
}}}
而把mi后面的tab扩展名去掉，也生成了和shapefile类似的目录结构。在ogr的库下，真的不能在一个tab文件中塞入很多不同类型的Layer吗？但是mapinfo确实可以做到一套tab文件中同时塞入点线面呀……

后来我试验了一下，打开生成的mitab.tab文件，加入了几个点后保存，这样图层就既有点要素又有多边形要素，然后用ogr打开，

{{{
#! python
>>> import ogr
>>> ds = ogr.Open('e:/shp2shp/mi/mitab.tab')
>>> ds.GetLayerCount()
1
}}}
还是1？哦，看来多层layer多种要素类型是过于理想化了。那么如何存储多种要素呢？

{{{
#! python
>>> layer = ds.GetLayer()
>>> f = layer.GetNextFeature()
>>> geom = f.GetGeometryRef()
>>> geom.GetGeometryType()
3
>>> f = layer.GetNextFeature()
>>> geom = f.GetGeometryRef()
>>> geom.GetGeometryType()
1
>>>
}}}
现在真相大白了。多层layer的单一文件是我们一厢情愿。在一个文件数据源中塞入多种feature不是在layer中动手脚而是在Feature的geometry中做的文章！

但是这样一来，layerdefn的形状类型是什么？

{{{
#! python
>>> laydef = layer.GetLayerDefn()
>>> laydef.GetGeomType()
0
>>> ogr.wkbUnknown
0
>>>
}}}
吼吼，当然什么也不是咯。

那么这么一来，能不能在shapefile里面也像mapinfo那样乱塞一些什么东西进去呢？

{{{
#! python
# -*- encoding:gb2312 -*-
import ogr,os,math
ogr.GetDriver
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'mi'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)#,options=['FORMAT=MIF'])
layernew = newds.CreateLayer('mitab',None,ogr.wkbUnknown)
fieldid = ogr.FieldDefn('FID',ogr.OFTInteger)
layernew.CreateField(fieldid)
extent = [0,0,300,300]
tw = 100
th = 100
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
        extent[0]+tw,extent[3], extent[0]+tw,extent[3]-th,
        extent[0],extent[3]-th, extent[0],extent[3])
print wkt
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
wktp = 'POINT (%f %f)' % (70.4,220.34)
geom2 = ogr.CreateGeometryFromWkt(wktp)
feat2 = ogr.Feature(layernew.GetLayerDefn())
feat2.SetGeometry(geom2)
layernew.CreateFeature(feat2)
newds.Destroy()
}}}
结果

{{{
#! python
C:\WINDOWS\system32\cmd.exe /c D:\Python24\python.exe create.py
POLYGON ((0.000000 300.000000,100.000000 300.000000,100.000000 200.000000,0.0000
00 200.000000,0.000000 300.000000))
ERROR 1: Attempt to write non-polygon (POINT) geometry to  type shapefile.
Hit any key to close this window...
}}}
看来不是什么都可以侥幸的。

后来又仔细想想，大悟！的确不可以把Layer想像成形状类型分层，这从设计逻辑上可以得到结果的；既然一个数据源可以是一个文件夹，也可以是一个文件，而文件夹下的文件被表示成数据源中的layer，那么就不可能再在单个文件的数据源中存在按形状分类型的layer，如果可以分，那么在一个文件夹数据源下打开一个多类型的文件Layer中的形状类型层应该怎么表示？

= 创建MapInfo的文件时的一些问题 =
[[http://ogr.maptools.org/drv_mitab.html|这里]]提到一些在创建!MapInfo过程中的问题： TAB格式在一个新的文件中写入第一个要素前需要先创建一个地理范围。

当时现在没有清楚的机制来通过OGRDataSource来为新文件设置缺省的框架。

我们应该在适合的驱动上为每个投影设置合理的范围，但是现在暂时用下面的默认框架来为mapinfo驱动在创建新Layer时设置框架范围。

 .* 一个经纬坐标系的框架范围BOUNDS (-180, -90) (180, 90)
 .* 其他的用BOUNDS (-30000000, -15000000) (30000000, 15000000)
如果在创建Layer时没有提供坐标系，使用了投影实例，如果坐标系真的是地理坐标，就没有可以导致非常低精度地理坐标系。你可以添加"-a_srs wgs84"到 ogr2ogr 命令行在一个转换中来强制为地理模型

Mapinfo有很多属性限制：

 .* 只可以创建Integer，Real，String字段，其他的类型不能创建。
 .* 字符串字段，字段宽被用来创建dat文件中的存储大小。这意味着比字段宽长得字符串会被阶段。
 .* 没有分配宽的字符串字段被对待为254个字符。
但是BOUNDS (-30000000, -15000000) (30000000, 15000000)这个框架的限制太不好了，由于中国的西安北京坐标系的投影的关系，还要东移500公里，再加个带号什么的，于是通常会超过30000000（福建在3度分带的西安或者北京坐标系上带号差不多在39－40，加上东移的500000m，最后的坐标值变成39XXXXXX，铁定超出(-30000000,30000000)这个范围框架）。那么用OGR创建的Mapinfo福建部分都是错的。这很不好。不过，我试了一下，发现mif可以正常创建，没有问题。还好还好……虽然离mapinfo能用还差一步之遥，虽然是个ACSII字符文件，稍微大了点……

= 添加投影 =
上面的!CreateLayer的第二个参数都是None，其实它是一个空间参考。可以添加它来让那些乱七八糟的数字有实际的含义。添加了空间参考，那些数字才能在地图上找到自己的位置。同一个Layer需要拥有相同的空间参考。

创建空间参考最简单的办法是找一个wkt来，然后用!ImportFromWkt来导入。Wkt到处都是啊。最简便好用的方法就是安装一个PostgreSQL，然后安装PostGIS(windows里面有个安装PostGIS的选项，勾起来就好;Linux下比较麻烦，一般需要自己编译，好像还没有很好用二进制安装包源，rpm和deb的我到现在都没有发现有哪个官方的源中有包含postgis的，可能有一些零散的二进制包，但是要解决依赖问题，算了，还是编译来得清爽！[[http://www.3snews.net/blog/user1/linux23/archives/2006/3584.html|这里]]有我写的编译和启动方法，虽然是在FC4下，但是其他系统也可以参考)。里面有spatial_reference的表（表名可能不叫这个，但和这个差不多，好像有个前缀，找一找吧，其实刚建好没有任何数据的PostGIS就两张表，一张登记数据集一张就是空间参考信息）里面有全套的空间参考，几乎所有的投影坐标信息的wkt表达都在里面。要哪个就拷贝出来，直接就能用。

假设已经有了一个wkt了，创建一个空间参考就像这样：

{{{
#! python
>>> import osr
>>> spatial = osr.SpatialReference()
>>> spatial.ImportFromWkt(wkt)
>>>
}}}
当然还有其他创建空间参考的方法，但是这就涉及很多专业知识，哪个下面要跟哪个，都要考虑半天的。[[http://wiki.woodpecker.org.cn/moin/lilin/gdal-coordinatetrans|这里]]有一些简单的解释，可以看看。

然后将spatial来代替!CreateLayer的第二个参数None，就可以建立矢量数据的空间参考了。

= 反馈 =
如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆[[http://groups.google.com/group/geosings?hl=zh-CN|这个论坛]]。