python读取raw数据文件_【Python】OGR库（1）：读取矢量数据

OGR库是一个非常流行的处理地理空间矢量数据的开源库。它可以读取丰富的数据格式，允许用户进行几何处理、属性表操作、数据分析，是个非常强大的开源GIS库。目前OGR已集成在GDAL库中，可以说是GIS的本源之一了，有大量的软件用到了这个库。本篇文章是关于OGR库的一些基础用法汇总，将会持续更新~

这里推荐两个比较好的OGR学习资源网站：Welcome to the Python GDAL/OGR Cookbook!​pcjericks.github.iohttps://gdal.org/python/​gdal.org

1. 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。整理一下可见下图：OGR类结构

2.3 读取矢量数据

通过对driver对象调用Open方法，可以打开一个文件，打开的结果就是一个Data Souce对象。通过print可以查看到这个信息shp_test是个osgeo.ogr.DataSource对象。

from osgeo import ogr

driver = ogr.GetDriverByName('ESRI Shapefile') #这个函数是不区分大小写的

filename = r'F:\Zhihu\DATA'

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 ogr

filename = r'F:\Zhihu\DATA'

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 pb

pb.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:\Zhihu\DATA') #设置工作空间

from ospybook.vectorplotter import VectorPlotter

vp = 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:

break

x=feat.geometry().GetX()

y=feat.geometry().GetY()

fid = feat.ID

area = 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')