九段线shp数据去哪里下载_【Python】OGR库（1）：读取矢量数据

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

这里推荐两个比较好的OGR学习资源网站：

Welcome to the Python GDAL/OGR Cookbook!​pcjericks.github.io https://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 

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 

也可以直接使用ogr.open函数打开文件。该函数会根据文件后缀名自动选择driver进行数据读取。

from 

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 

2.5.2 查看图形

ospybook 提供了用于绘图的类。它是基于matplotlib的，所以必须安装matplotlib。

交互式图像在创建之后无需使用draw函数调用就可呈现。声明方法是建立一个VectorPlotter类的新实例。

import 

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 

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