2. GDAL python教程（1）——用OGR读写矢量数据¶

本教程的讲义和源码都是取自Utah State University的openGIS课程

相关资料，包括讲义、源码、数据样例，请从此处下载

本人只是做点翻译，写写学习体会而已，版权属于原作者。

欢迎转载，不过别忘了上面这段话。

2.1. 为什么用open source？¶

优点

免费，适合个人和小公司
强大的开发工具，找bug更容易
跨平台，windows和linux都能用
拉风！

缺点

没有内嵌地理处理器
用的人少

2.2. Open source RS/GIS模块¶

OGR矢量库：简单的矢量数据读写，是GDAL的一部分
GDAL地理空间数据抽象库： a) 读写栅格数据 b) ArcGIS也是基于GDAL开发的 c) C++库，但是可以用python调用

2.3. 相关模块¶

Numeric：高速的数组处理，对栅格数据尤其重要
NumPy：下一代的Numeric
更强大的gis库 http://www.gispython.org/

2.4. 导入库：¶

import ogr

或者：

from osgeo import ogr

万能的方法是：

try:
    from osgeo import ogr
except:
    import ogr

要读取某种类型的数据，必须要先载入数据驱动，也就是初始化一个对象，让它“知道”某种数据结构。

import ogr
driver = ogr.GetDriverByName(‘ESRI Shapefile’)

数据驱动driver的open()方法返回一个数据源对象

open(<filename>, <update>)

其中update为0是只读，为1是可写

例如：

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
filename = 'C:/Users/gongwei/Documents/My eBooks/python_and_sage/GDAL python/test/ospy_data1/sites.shp'
dataSource = driver.Open(filename,0)
if dataSource is None:
    print 'could not open'
    sys.exit(1)
print 'done!'

注意filename一定要写绝对路径！

因为一定要用绝对路径，为了简化代码，经常会使用到os.chdir()

2.5. 读取数据层¶

layer = dataSource.GetLayer(0)

一般ESRI的shapefile都是填0的，如果不填的话默认也是0.

再看看这个数据层里面有几个点呢？

n = layer.GetFeatureCount()
print 'feature count:', n

2.5.1. 读出上下左右边界¶

extent = layer.GetExtent()
print 'extent:', extent
print 'ul:', extent[0], extent[3]
print 'lr:', extent[1], extent[2]

读取某一要素feature（总算切入正题了），这里读取的是一个点

feat = layer.GetFeature(41)
fid = feat.GetField('id')
print fid
feat = layer.GetFeature(0)
fid = feat.GetField('id') #should be a different id
print fid

另外还有按顺序读取feature，循环遍历所有的feature

feat = layer.GetNextFeature()  #读取下一个
while feat:
    feat = layer.GetNextFeature()
later.ResetReading()  #复位

2.5.2. 提取feature的几何形状¶

geom = feat.GetGeometryRef()
geom.GetX()
geom.GetY()
print geom.

2.6. 释放内存¶

feature.Destroy()

关闭数据源，相当于文件系统操作中的关闭文件

dataSource.Destroy()

2.7. 读完了再说怎么写¶

2.7.1. 创建新文件¶

driver.CreateDataSource(<filename>)

但是这个文件不能已经存在了，否则会出错

创建新的layer

dataSource.CreateLayer(<name>,CreateLayer(<name>, geom_type=<OGRwkbGeometryType>, [srs])

举个例子：

ds2 = driver.CreateDataSource('test.shp')
layer2 = ds2.CreateLayer('test', geom_type=ogr.wkbPoint)

要删除一个shp文件

driver.DeleteDataSource(‘test.shp’)

要添加一个新字段，只能在layer里面加，而且还不能有数据

添加的字段如果是字符串，还要设定宽度

fieldDefn = ogr.FieldDefn('id', ogr.OFTString)
fieldDefn.SetWidth(4)
layer.CreateField(fieldDefn)

添加一个新的feature，首先得完成上一步，把字段field都添加齐了

然后从layer中读取相应的feature类型，并创建feature

featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)

设定几何形状

feature.SetGeometry(point)

设定某字段的数值

feature.SetField('id', 23)

将feature写入layer

layer.CreateFeature(feature)