Python读写矢量数据（1）针对读取矢量数据——Python地理数据处理学习分享

Cishazhu_ooook

已于 2022-10-17 20:48:10 修改

阅读量2.3k

点赞数 7

分类专栏： Python地理数据处理学习分享文章标签： python

于 2022-10-16 21:01:18 首次发布

本文链接：https://blog.csdn.net/remote_giser/article/details/127351424

版权

Python地理数据处理学习分享专栏收录该内容

8 篇文章 4 订阅

订阅专栏

1.矢量数据存储格式

常用的矢量存储格式包括两类。Shapefile是存储矢量的通用格式，它至少包括三个文件，而且文件用途各不相同。几何信息存储在.shp和.shx文件中，属性信息存储在.dbf文件中。另一种广泛的格式，尤其针对网络地图是GeoJSON格式，他是纯文本文件，所有信息只存储在这一个文件中。

2.GDAL/OGR

GDAL是地理空间数据抽象库，是Python处理矢量、栅格文件基础的模块。OGR是GDAL的一部分，它并不代表任何东西。GDAL的OGR部分具有读写许多不同矢量数据格式的功能。OGR还允许你创建和操作要素的几何形状、编辑其属性值、对矢量数据进行筛选等操作。如果想使用GDAL/OGR需要安装GDAL库，GDAL库安装参考以下文章。Python环境下GDAL库的安装（Windows系统）_Cishazhu_ooook的博客-CSDN博客https://blog.csdn.net/remote_giser/article/details/127345752?spm=1001.2014.3001.5501 GDAL/OGR支持多种数据格式（矢量数据并不包括以上两种），想了解OGR支持的矢量格式可以浏览以下。OGR/GMT矢量数据格式 — GMT 中文手册 (gmt-china.org)https://docs.gmt-china.org/latest/table/ogrgmt/ 此外，我们还应该知道OGR打开数据源是，会有一个数据源（DataSource）对象。此数据源可以有一个或多个子图层对象，每一个图层代表数据源中的一个数据集。数据集中会包含数据的各个特征。

3.读取矢量数据

3.1简单的读取矢量数据.

在arcgis中打开.shp文件查看，如图。（本文用到的数据在附录下载后globe文件下）

使用简单的Python代码达到读取矢量数据的目的。

# 导入模块
import sys
from osgeo import ogr
# 读取文件路径
fn="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
# 打开文件
ds=ogr.Open(fn,0)
if ds is None:
	sys.exit('Could not open{0}.'.format(fn))
# 打开图层
lyr=ds.GetLayer(0)
# 利用for循环进行打印特征
i=0
for feat in lyr:
	pt=feat.geometry()
	x=pt.GetX()
	y=pt.GetY()
	name=feat.GetField('NAME')
	pop=feat.GetField('POP_MAX')
	print(name,pop,x,y)
	i+=1
	if i==10:
		break
# 关闭文件
del ds

我们成功读取到了文件中前10个点name、pop、x/y信息，运行结果如下：

值得一提的是，这里我们导入了两个模块，sys与ogr。sys是针对与python系统的模块，在读取文件失败后起作用，如果不需要可以自行删除，这里必备的模块ORG在程序运行中起关键作用。（个人不太熟悉，想学习的可以自行搜素。应用在本例修改代码的路径会提示路径错误。）对于以上代码在正常运行情况下可以删掉sys模块和第一个if条件语句。

针对这个程序，我们首先要做的事打开这个程序，令一个变量为文件的路径，在这里我们使用变量fn。调用ds=ogr.Open(fn,0)打开数据源。这里fn我们前面已经叙述过，后面的0代表文件以只读的形式打开。OGR的Open函数会传递文件名和一个可选的更新标记，可选变量标记可以是0或1，0代表只读模式打开，仅能进行读取数据工作；而1代表编辑模式打开，可以对内容进行编辑。如果第二个参数不设定，则会以只读模式打开。例：这里我将代码改为ds=ogr.Open(fn)文件正常运行，因为此程序仅为只读，我们也同样以只读模式打开的文件。

别忘了，上面我们提到过OGR打开数据源后会出现多个图层，因此我们需要获得具体图层。数据源（ds）有一个GetLayer函数，可用用来索引所需要的图层。图层索引是从0开始的，所以第一层的索引为0（类似于列表）。同样如果不给GetLayer设置参数，默认索引第一层。在本例中我们使用lyr=ds.GetLayer(0)索引了第一层。（其实本例就一个图层）

前面只是对数据的打开，接下来才是文件的真正读取。图层构成的要素，可以用for循环来遍历读取。（本例只读取了前10个值）

首先我们要获得几何位置X/Y。先用一个变量（pt）存储他的几何位置。在这里用.geometry（）函数获取了他的几何特征，并用.GetX/Y()来获取x/y。
文本特征这里用.GetField()来获取。但是值得注意的时该函数返回的数据类型与地层数据集中的数据类型一致。本例中name返回的的是字符型，pop返回的为整型。如果想让pop也返回数值型，可使用.GetFieldAsString()函数。
最后使用del ds关闭打开的整个文件。

3.2查看图层要素属性值字段

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# 查看可用操作
use=dir(layer)
print(use)
# 对图层属性查询，包括值字段等
layerdef=layer.GetLayerDefn()
for i in range(layerdef.GetFieldCount()):
 	defn = layerdef.GetFieldDefn(i)
 	print(defn.GetName(),defn.GetWidth(),defn.GetType(),defn.GetPrecision())
 # 关闭文件
del docm

我们定义了一个变量use，采用dir查看可用的函数。如下图我们可以看到对图层layer可以使用的函数，当然该命令对图层的属性数据查询无关，我们仅仅想查看一下可用的函数，

查询图层要输属性值字段操作较为简单，首先使用函数 .GetFieldDefn() 进行定义图层信息。后面使用for循环进行属性的查询。GetFieldCount()的目的是查询获得的属性量。.GetFieldDefn(i)获取第i个属性字段，.GetName().GetType().GetWidth().GetPrecision()分别对应获取字段名、宽度、类型、精度。

如果想要查询shapefile文件某要素属性信息（可能会包含大量的的属性），但我们不知道具体有什么要素，可以通过python进行要素属性字段的检索，然后使用3.1的方法找出自己要寻找的信息。

上述查询比较繁琐，这里提供了一种新方法。（建议采用这种，上方只做了解）

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# 查看属性表
for field in layer.schema:
    print(field.name,field.GetTypeName())
 # 关闭文件
del docm

最简单的方法就是通过图层对象的schema属性来获得FieldDefn对象列表。每一个对象包含诸如属性列表名称和类型的信息。

3.3查看属性（Feature）

除了使用3.1提过的.GetField()方法来查询，Feature的方法也是一种查看属性的方法，但这两种都需要要素的属性值字段，因此3.2的操作来获取字段相当重要。

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# Feature方法查询
num_features=layer.GetFeatureCount()
for i in range(num_features):
	feature=layer.GetFeature(i)
	print(feature.POP_MAX,feature.NAME)
	i+=1
 # 关闭文件
del docm

首先我们要做的是了解一共有多少要素。我们用到的是.GetFeatureCount()方法（觉得更新获得了一个从0至X的数值列表，列表又嵌套了字典，你大可这样理解），在这里num_features=layer.GetFeatureCount()能很快帮我们了解有多少要素，如果我们打印查看一下就可知道与在GIS软件中打开的要素数一致（FID从0计数）。显然这里查询结果为1249，GIS打开属性表的FID也1248，这个实例能让你明白查询的计数到底是什么了。

print(num_features)

然后我们继续使用for循环来实现查询。.GetFeature()用来获取属性，在打印过程中我们需要指出feature要获取什么属性，例如feature.NAME获取了NAME属性。

当然Feature不同于Field的其中之一是可以访问特定数据。

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# Feature指定查询最后一个要素的NAME
num_features=layer.GetFeatureCount()
last_feature=layer.GetFeature(num_features-1)
print(last_feature.NAME)
# Feature指定查询第一个要素的NAME
feature=layer.GetFeature(0)
print(feature.NAME)
 # 关闭文件
del docm

注意在从后面查询的时候一定不要忘记给要素计数，num_features=layer.GetFeatureCount()起了关键作用。在这里我们指定了特定的位置num_feature-1也就是num_feature最后一个值，在打印查询结果时记得不要忘记要查询的属性（NAME）。当然按正序指定要素对于要素的计数不做要求。

3.4获取几何信息

在3.1中我们已经或多或少了解了一部分几何信息的获取，例如一个点在图层范围内的x/y

坐标，下面我们将继续学习一些关于获取几何信息的知识。

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# 获取一个图层的边界范围
extent=layer.GetExtent()
print(extent)
 # 关闭文件
del docm

以上的例子我们对layer使用.GetExtent()函数，另变量得到一个数字元组（x最小值，x最大值，y最小值,y最大值），再使用print语句打印这个变量，能够很直观的得到结果。

当然，我们可能想知道上述文件1249个要素气质一个要素的类型是什么，它到底是点还是线还是面？在这里我们也是有方法解决的。

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# 查询第100个元素的类型
feat=layer.GetFeature(100)
print(feat.geometry().GetGeometryName())
 # 关闭文件
del docm

这里我们采用了一个geometry().GetGeometryName()函数就轻而易举的获取了特定要素的类型，大家也可以尝试改变要素，比如要素10、20、39......甚至可以使用for循环来查找每个要素的类型。或许大家考虑过同一图层下都是点文件呢（使用gis软件时，常常把相同的shape作为一个图层），那我们也可以尝试查询一下整个图层的属性。

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# 获取整个图层的要素类型
lx=layer.GetGeomType()
print(lx)
print(layer.GetGeomType()==ogr.wkbPoint)
print(layer.GetGeomType()==ogr.wkbPolygon)
 # 关闭文件
del docm

运行结果

在这里大家可能会疑惑为什么会出现这个结果，因为我们使用GetGeomType（）函数时返回的时一个字符串（点对应1、线对应2、面对应3），此外layer.GetGeomType()==ogr.wkbPoint表达的意思的获取的layer要素是否与wkpPoint一致，返回一个布尔值（True/False）。OGR有多个常量对应如下：

数据的空间参考系统，是用来描述数据集使用的坐标系统，未经投影的坐标显示的是我们熟知的经纬度坐标，这些坐标系可以转换成各种投影坐标系，以便我们解读。然而如果不知道数据集使用的是哪个坐标系统这将是个很麻烦的事情。

# 导入模块
from osgeo import ogr
# 打开文件
path="D:\\Python_study\\osgeopy-data\\global\\ne_50m_populated_places.shp"
docm=ogr.Open(path,0)
# 加载图层
layer=docm.GetLayer(0)
# 查看坐标系统信息
print(layer.GetSpatialRef())
 # 关闭文件
del docm