欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  后端开发

python gdal教程之:用ogr读写矢量数据

程序员文章站 2022-04-21 21:29:08
...
为什么用open source?

优点

1. 免费,适合个人和小公司

2. 强大的开发工具,找bug更容易

3. 跨平台,windows和linux都能用

4. 拉风!

缺点

1. 没有内嵌地理处理器

2. 用的人少

Open source RS/GIS模块

1. OGR矢量库:简单的矢量数据读写,是GDAL的一部分

2. GDAL地理空间数据抽象库:

a) 读写栅格数据

b) ArcGIS也是基于GDAL开发的

c) C++库,但是可以用python调用

相关模块

1. Numeric:高速的数组处理,对栅格数据尤其重要

2. NumPy:下一代的Numeric

3. 更强大的gis库 http://www.gispython.org/

导入库:

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()

读取数据层

layer = dataSource.GetLayer(0)

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

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

n = layer.GetFeatureCount()

print 'feature count:', n

读出上下左右边界

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() #复位

提取feature的几何形状

geom = feat.GetGeometryRef()

geom.GetX()

geom.GetY()

print geom.

释放内存

feature.Destroy()

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

dataSource.Destroy()

读完了再说怎么写

创建新文件

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)

以上就是python gdal教程之:用ogr读写矢量数据的内容,更多相关内容请关注PHP中文网(www.php.cn)!