OGR库是一个非常流行的处理地理空间矢量数据的开源库。它可以读取丰富的数据格式,允许用户进行几何处理、属性表操作、数据分析,是个非常强大的开源GIS库。目前OGR已集成在GDAL库中,可以说是GIS的本源之一了,有大量的软件用到了这个库。本篇文章是关于OGR库的一些基础用法汇总,将会持续更新~
这里推荐两个比较好的OGR学习资源网站:
Welcome to the Python GDAL/OGR Cookbook!pcjericks.github.iohttps://gdal.org/python/gdal.org1. GDAL库的安装
首先要去下面的网站
https://www.lfd.uci.edu/~gohlke/pythonlibs/www.lfd.uci.edu找到GDAL对应的地址
下载gdal对应python版本的whl文件,这里因为我是python 3.7 64位的,所以下载相应版本就行了:
下载到本地后,对文件存放的文件夹按住shift+鼠标右键,选择“在此处打开Powershell窗口”
然后在弹出来的对话框中键入“pip install 文件名”按回车,等待安装即可
在安装好之后,会成功安装osgeo这个包。OGR模块是在osgeo包中的。这个包里的所有模块都是以小写字母命名。
2. 读取矢量数据
2.1 定义Driver
在正式读写之前得先理解一下OGR的对象组织形式,弄清楚OGR类之间的关系。
在打开一个矢量数据前,首先需要定义一个driver,用来告诉OGR你想要读取何种格式的数据。每个driver就是对应一种数据格式,比如‘ESRI Shapefile’就是对应的shp格式。OGR可以读取70多个数据格式,基本常用的都包含了。下面是定义driver的一个示例:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
2.2 OGR类结构
在讲如何打开数据之前,最好先了解一下OGR的对象组织逻辑。
使用OGR打开一个矢量数据,如shp文件或GeoJSON文件,会产生一个DataSource对象。该对象包含若干个Layer,每个Layer就是一个要素集。值得注意的是,大多数矢量数据格式一般只有一个Layer,少数有多个(如SpatiaLite格式)。既然Layer是个要素集,可想而知它包含的就是一个个的feature了。Feature的概念稍微学过GIS原理的朋友应该都知道,Arcmap里也经常提到,就是几何对象和属性表的集和。所以,feature包含的就是geometry和attribute。整理一下可见下图:
2.3 读取矢量数据
通过对driver对象调用Open方法,可以打开一个文件,打开的结果就是一个Data Souce对象。通过print可以查看到这个信息shp_test是个osgeo.ogr.DataSource对象。
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile') #这个函数是不区分大小写的filename = r'F:ZhihuDATA'
shp_test = driver.Open(filename,0) #0是只读,1是可写
if shp_test is None:print ('could not open')sys.exit(1)
print(shp_test)
也可以直接使用ogr.open函数打开文件。该函数会根据文件后缀名自动选择driver进行数据读取。
from osgeo import ogrfilename = r'F:ZhihuDATA'
shp_test = ogr.Open(filename)
2.4 读取要素个数
OGR的Layer概念类似于ArcGIS里的FeatureClass,就是多个同类要素(点、线、多边形)的集和。
可以通过dir函数获取Layer可以使用的方法。这里使用GetLayer方法获取shp数据的图层,再对其使用GetFeatureCount方法可以获取shp数据中的元素个数。
layer = shp_test.GetLayer()
dir(layer)
n = layer.GetFeatureCount()
print ('feature count:', n)
2.5 查看数据
2.5.1 查看属性
可在cmd中使用pip install https://github.com/cgarrard/osgeopy-code/raw/master/ospybook-latest.zip安装ospybook包。这个包是书《Geoprocessing with python》作者开发的。可以使用print_attributes来打印出数据属性表,可以直接输入文件名也可以输入layer。
不建议使用这个函数来打印大型数据的全部属性表,可使用数字来控制打印几行并指定打印的字段名。
import ospybook as pbpb.print_attributes(filename,3,['Shape_Area','ID','long','lat']) pb.print_attributes(lyr,3,['Shape_Area','ID','long','lat']) #既可以调用文件名,也可以调用图层,效果等价
2.5.2 查看图形
ospybook 提供了用于绘图的类。它是基于matplotlib的,所以必须安装matplotlib。
交互式图像在创建之后无需使用draw函数调用就可呈现。声明方法是建立一个VectorPlotter类的新实例。
import os
os.chdir(r'F:ZhihuDATA') #设置工作空间
from ospybook.vectorplotter import VectorPlottervp = VectorPlotter(True) #声明图像为交互式图像
vp.plot(filename,fill=False) #fill是用来控制是否填充要素,默认是true
2.6 读取四至范围
使用GetExtent方法获得四至信息,结果是一个元组,顺序为左右下上。若shp数据本身含有投影坐标,则输出的也是投影坐标系的值。
extent = layer.GetExtent()
print ('extent:', extent)
print ('x range:', extent[0], extent[1])
print ('y range:', extent[2], extent[3])
2.7 读取单个要素
- 使用GetFeature方法按照FID读取要素,这里读取的第二个要素,即FID=1的那个要素。
- 通过GetField方法可读取要素指定列信息,值得注意的是这里需要输入的列名不分大小写,同shp格式的要求一致。
feat = layer.GetFeature(1)
fid = feat.GetField('id')
area = feat.GetField('shape_area') print (fid)
print (area)
dir(feat)
2.8 遍历要素
- 使用GetNextFeature方法可以省去使用For循环按ID读取的低效率;
- 要使用个try except机制,不然再最后一个要素读完之后,GetNextFeature方法仍然会读下一个空要素,这时输出面积会报错。
- ResetReading函数是用来复位的,不然下次使用GetNextFeature程序接着上次读的位置继续读。
feat = layer.GetNextFeature() #读取下一个
while feat:feat = layer.GetNextFeature()try:area = feat.GetField('shape_area') print(area)except:print('Done!')
layer.ResetReading() #复位
也可以通过for循环直接遍历每个要素。但值得注意的一点是当前要素问题。比如我们用For循环遍历了一遍Layer中的所有feature并输出了它们一个字段的值。若想再通过for循环遍历一遍输出另一个字段的值,你会发现得不到任何结果。这是因为在第一次for循环后,指针停在了最后一个feature上,必须使用Layer.ResetReading()函数来进行重置。
for i,feat in enumerate(layer):if i >= 5:breakx=feat.geometry().GetX()y=feat.geometry().GetY()fid = feat.IDarea = feat.GetField('shape_area') print (fid,area,x,y)
layer.ResetReading() #复位
2.9 提取要素几何信息
- 使用GetGeometryRef方法读取要素几何信息,通过dir函数可以查看geom可以使用的方法;
- GetX和GetY可以直接打印一个个点的坐标;
- 使用geom.Area()可以读feature的面积,默认单位为㎡。
feat = layer.GetFeature(1)
geom = feat.GetGeometryRef()print(geom)
print(geom.Area())
2.10 释放内存
- 要素.Destory是先关闭单个要素,后面的Destory是关闭整个DataSource;
- 关闭数据源,相当于文件系统操作中的关闭文件。
feat.Destroy()
shp_test.Destroy()
2.11 删除文件
使用DeleteDataSource可以删除shp文件及其附属文件(如dbf,poj等文件)。
import os
filename = 'F:/Zhihu/DATA/testCopy.shp'
if os.path.exists(filename):driver.DeleteDataSource(filename)print('File was deleted!')
else:print('File not exist')