>>> from env_helper import info; info()

页面更新时间： 2023-06-23 11:54:51
运行环境：
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-9-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

4.2. 读取矢量数据¶

矢量数据常用的数据集一种是ESRI的Shapefile，另一种是MapInfo的TAB文件。它们均包含了数个文件，且存放于同一目录，几个同名不同后缀的文件共同组成了一个矢量数据集，当然，不同文件的作用也不同，这点与栅格数据集类似。例如在Shapefile中，至少有 .shp 、.dbf 与 .shx ；在TAB文件中，至少有 .dat 、 .id 、 .map 与 .tab 。通常，OGR处理文件是以数据集为对象的，打开矢量数据时，必须打开以.shp/.tab为后缀的文件。

下面我们就以Shapefile文件为例，运用OGR来操作矢量数据。 Shapefile是最常见且设计上也非常简洁的一种模型，更详细内容可查看ESRI相关说明。我们所使用的数据是标准教程中的Country数据，是一幅最常见的世界地图。

4.2.1. 直接打开数据¶

可使用 ogr.Open() 方法直接打开矢量数据，在这个过程中，ogr会自动根据文件的类型来确定相应的驱动。

>>> inshp = '/gdata/GSHHS_c.shp'
>>> from osgeo import ogr
>>> datasource = ogr.Open(inshp)
>>> inshp = 'shape_towns.shp'
>>> driver = datasource.GetDriver()
>>> driver.name

'ESRI Shapefile'

这样就打开了一个数据源（DataSource），并将其赋值给datasource变量。

使用Python的内省函数dir()查看datasource方法与属性：

>>> dir(datasource)[:6] + ['... ...'] + dir(datasource)[-3:]

['AbortSQL',
 'CommitTransaction',
 'CopyLayer',
 'CreateLayer',
 'DeleteLayer',
 'Dereference',
 '... ...',
 'name',
 'this',
 'thisown']

4.2.2. OGR的数据驱动¶

注意，上述代码是按缺省方式打开的。在实际应用中，应根据打开的数据类型来初始化数据驱动：

>>> driver = ogr.GetDriverByName('ESRI Shapefile')

该函数的原型：

Open(self, char name, int update = 0) -> DataSource

根据数据驱动 driver 的 Open() 方法返回一个数据源对象，update 0为只读，1为可写：

>>> import sys
>>> inshp = '/gdata/GSHHS_c.shp'
>>> driver = ogr.GetDriverByName('ESRI Shapefile')
>>> datasource = driver.Open(inshp,0)
>>> if datasource is None:
>>>     print('could not open')
>>> else:
>>>     print('done')

done

这一节看一下如何使用 ogr 模块来获取图层的信息。

4.2.3. 图层的概念¶

在GIS中，Layer（图层）由同种要素（Feature）（如点、线、多边形等）组在一起的“层”。Layer的描述与ESRI在ArcGIS中定义的模型要素类(Feature Class)一致，其对应的是要素数据集。

在《Modeling our World》中，这两个概念均有明确的阐述：

要素类：具有相同几何形状的要素的集合：点、线或多边形。我们最常见的两种要素类是简单要素类和拓扑要素类。简单要素类是指没有任何拓扑关系的点、线、多边形或注记。也就是说，一个要素类内的点与另一要素类中线要素的终点可以是一致的，但它们是不同的，其特点是可独立编辑。拓扑要素类局限在一定的图形范围内，它是一个由完整拓扑单元组成的一组要素类限定的对象。
要素数据集（要素集）：具有相同坐标系统的要素类的集合。我们可以选择在要素集的内部或外部组织简单要素类，但拓扑要素类只能在要素集内部组织，以确保它们具有相同的坐标系统。

4.2.4. Python中操作Layer方法¶

Layer方法集合如下：

>>> from osgeo import ogr
>>> inshp='/gdata/GSHHS_c.shp'
>>> datasource = ogr.Open(inshp)
>>> layer = datasource.GetLayer(0)
>>> dir(layer)[:6] + ['... ...'] + dir(layer)[-4:]

['AlterFieldDefn',
 'AlterGeomFieldDefn',
 'Clip',
 'CommitTransaction',
 'CreateFeature',
 'CreateField',
 '... ...',
 '__weakref__',
 'schema',
 'this',
 'thisown']

注意：GetLayer的参数是从0开始的。对于Shapefile而言，它只有一个图层，一般情况下这个参数都是0，如果空着，缺省情况下也是0。

查看图层要素：

>>> layer.GetFeatureCount()

还可以查看图层范围：

>>> layer.GetExtent()

(-180.0, 180.00000000000023, -89.99999999999994, 83.53036100000006)

注意这些坐标的顺序，四个值分别是“西、东、南、北”四至，在其他软件中，这些顺序可能会不一致。

4.2.5. 图层的属性¶

若想查看整个表的结构及各字段名称等信息，可在layer的附加信息里查找。

>>> layerdef = layer.GetLayerDefn()
>>> for i in range(layerdef.GetFieldCount()):
>>>     defn = layerdef.GetFieldDefn(i)
>>>     print(defn.GetName(),defn.GetWidth(),defn.GetType(),defn.GetPrecision())

id 80 4 0
level 10 12 0
source 80 4 0
parent_id 10 12 0
sibling_id 10 12 0
area 19 2 11
Shape_Leng 19 2 11
Shape_Area 19 2 11

这就是整个表的结构。当然类型是枚举类型，若要查看数据类型名，可对应地建立一个数据类型字典。

4.2.6. 获取要素(Feature)信息¶

地理要素是地图的地理内容，包括自然地理要素与社会经济要素。其中，自然地理要素是指地球表面自然形态所包含的要素，如地貌、水系、植被和土壤等；社会经济要素是指人类在生产活动中改造自然界所形成的要素，如居民地、道路网、通讯设备、工农业设施、经济文化和行政标志等。

如果单从计算机角度看，要素就是一些几何形状，其实这种观点是欠妥的。几何形状，具体来说包括点、线、面、多边形、弧段、向量、控制点等多种类型，其实就是要素的抽象模型。

在Shapefile中，要素模型由点、线、面三种类型构成。要素类自带属性信息，一般对应属性表中的一行。

Layer的核心是针对Feature的操作。GetFeatureCount() 和 GetFeature() 或 GetNextFeature()可获取层内所有的要素；GetFeaturesRead()可查询目前已读取了多少条Feature；生成器可读取要素值；若要再次访问获取过的要素，可访问ResetReading()从头读取。

读取图层中的要素¶

下面看一下如何获取图层中的要素（Feature）。

>>> from osgeo import ogr
>>> inshp ='/gdata/GSHHS_c.shp'
>>> datasource = ogr.Open(inshp)
>>> layer = datasource.GetLayer(0)
>>> feature = layer.GetFeature(0)
>>> dir(feature)[:6] + ['... ...'] +  dir(feature)[-3:]

['Clone',
 'Dereference',
 'Destroy',
 'DumpReadable',
 'Equal',
 'ExportToJson',
 '... ...',
 'keys',
 'this',
 'thisown']

按顺序读取feature，循环遍历所有的feature：

>>> layer.ResetReading()
>>> feature = layer.GetNextFeature()
>>> while feature:
>>>     feature = layer.GetNextFeature()
>>> dir(feature)

['__bool__',
 '__class__',
 '__delattr__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getstate__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__']

ResetReading() 为重置，以便再次获取。GetNextFeature()为获取下一个要素。

feature的有关方法：’Clone’, ’Destroy’, ’DumpReadable’, ’Equal’, ’GetDefnRef’, ’GetFID’, ’GetField’ , ’GetFieldAsDouble’, ’GetFieldAsInteger’, ’GetFieldAsString’, ’GetFieldCount’, ’GetFieldDefnRef’, ’GetFieldIndex’, ’GetGeometryRef’, ’GetStyleString’, ’IsFieldSet’, ’SetFID’, ’SetField’, ’SetFrom’, ’SetGeometry’, ’SetGeometryDirectly’, ’SetStyleString’, ’UnsetField’。

其中，Field是属性操作。Geometry是形状操作。

读取要素属性¶

下面是一个简单的例子：

>>> layer.ResetReading()
>>> feat=layer.GetFeature(0)
>>> feat.keys()

['id',
 'level',
 'source',
 'parent_id',
 'sibling_id',
 'area',
 'Shape_Leng',
 'Shape_Area']

>>> feat.GetField('level')

以上是使用feature.GetField('CNTRY_NAME')来获取字段的值。字段的名称，在上节图层属性中获取，也可以使用要素keys()方法获取。

除了使用字段名称之外，还可以使用索引值来获取要素属性。下面是对所有的属性值进行遍历：

>>> for i in range(feat.GetFieldCount()):
>>>      print(feat.GetField(i))

0-E
1
WVS
-1
0
50654050.6945
2284.13578022
6280.55928902

此操作列出的是这个要素的所有属性值。

要素的几何形状(Geometry)¶

>>> geom = feat.GetGeometryRef()
>>> geom.GetGeometryName()

'POLYGON'

>>> geom.GetGeometryCount()

>>> geom.GetPointCount()

>>> geom.ExportToWkt()[:80]

'POLYGON ((107.544833 76.914028,106.533222 76.499556,111.093417 76.7653890000001,'

这是Geometry的主要操作方法。这些操作通过循环已经可以把Aruba的国家轮廓画出来了。需要注意的是如果Geometry类型是多边形，那么就有两层Geometry，获取多边形后，还需要获取构成多边形边缘的折线，才可以读取构成多边形的点。如果类型是单线或者点，就可以直接用GetX()， GetY()，而不用再多花一步GetGeometryRef()来获取子形状了。

>>> polygon=geom.GetGeometryRef(0)
>>> polygon.GetGeometryName()

'LINEARRING'

>>> polygon.GetGeometryCount()

>>> polygon.GetPointCount()

>>> polygon.GetX(0)

107.54483300000004

>>> polygon.GetY(0)

76.91402800000003

>>> polygon.GetZ(0)

0.0

>>> polygon.ExportToWkt()[:80]

'LINEARRING (107.544833 76.914028,106.533222 76.499556,111.093417 76.765389000000'

需要注意的是这些点都是地理意义上的坐标点，而不是数据意义的坐标。单位是根据空间参考决定的（此数据的单位应该是经纬度，而不是投影坐标的单位米）。获取了这些坐标后，不需要像栅格数据那样用四周边界点来计算地理坐标，矢量数据可直接铺展到一个实际平面已经定义好的投影或者地理坐标系上。

栅格数据有像元长宽，并且有数据意义上的行列坐标，而矢量在没有贴到某个实际存在的物体表面-如计算机屏幕或者打印的纸张之前，是无法表示数据意义上的长宽，所以比例尺对于矢量数据来说，没有绘制操作也就无根本意义。

4.2.7. 关闭数据与释放内存¶

在读取完数据之后，应保持良好习惯，释放数据，及时清理操作所占内存。

释放内存需要将要素的资源释放，使用Destroy()函数。

>>> from osgeo import ogr
>>> inshp ='/gdata/GSHHS_c.shp'
>>> datasource = ogr.Open(inshp)
>>> layer = datasource.GetLayer(0)
>>> feature = layer.GetFeature(0)
>>> feature.Destroy()

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

>>> from osgeo import ogr
>>> inshp ='/gdata/GSHHS_c.shp'
>>> datasource = ogr.Open(inshp)
>>> datasource.Destroy()

4.2.8. 矢量数据的空间参考¶

矢量数据的空间参考信息可从每个要素中获取。ShapeFile的坐标系统和仿射参数是分开的，通过.prj附加在shp文件同目录下的方式来获取坐标系统，而仿射参数是直接放在ShapeFile文件内部的，只不过它和GeoTiff的tfw不一样，它表示一个外包矩形，记录的是整个layer的地理范围。

>>> from osgeo import ogr
>>> inshp = '/gdata/GSHHS_c.shp'
>>> datasource = ogr.Open(inshp)
>>> layer = datasource.GetLayer(0)
>>> layer.GetSpatialRef()
>>> layer.GetExtent()

(-180.0, 180.00000000000023, -89.99999999999994, 83.53036100000006)

除了整个图层有地理范围外，每个要素（feature）也有地理范围，可通过GeometryDef来获取这个地理范围：

>>> feature = layer.GetFeature(0)
>>> geom = feature.GetGeometryRef()
>>> geom.GetEnvelope()

(-9.500388999999927, 180.0000000000001, 1.2695000000000505, 77.71625000000006)

>>> geom.GetSpatialReference()
>>> dir(geom)[:6] + ['... ...'] + dir(geom)[-3:]

['AddGeometry',
 'AddGeometryDirectly',
 'AddPoint',
 'AddPointM',
 'AddPointZM',
 'AddPoint_2D',
 '... ...',
 '__weakref__',
 'this',
 'thisown']

这里获取的外包矩形只针对于这个要素，而不是整个图层。

4.1. 使用OGR模块打开矢量数据

4.3. 使用 OGR 创建矢量数据

Python与开源GIS