矢量数据常用的数据集一种是ESRI的Shapefile,另一种是MapInfo的TAB文件。
它们均包含了数个文件,且存放于同一目录,
几个同名不同后缀的文件共同组成了一个矢量数据集,当然,不同文件的作用也不同,这点与栅格数据集类似。
例如在Shapefile中,至少有 .shp
、.dbf
与 .shx
;
在TAB文件中,至少有 .dat
、 .id
、 .map
与 .tab
。
通常,OGR处理文件是以数据集为对象的,打开矢量数据时,必须打开以.shp
/.tab
为后缀的文件。
下面我们就以Shapefile文件为例,运用OGR来操作矢量数据。 Shapefile是最常见且设计上也非常简洁的一种模型,更详细内容可查看ESRI相关说明。 我们所使用的数据是标准教程中的Country数据,是一幅最常见的世界地图。
直接打开数据
可使用 ogr.Open()
方法直接打开矢量数据,在这个过程中,ogr会自动根据文件的类型来确定相应的驱动。
inshp = '/data/gdata/GSHHS_c.shp'
from osgeo import ogr
datasource = ogr.Open(inshp)
driver = datasource.GetDriver()
# driver.name
这样就打开了一个数据源(DataSource),并将其赋值给datasource
变量。
使用Python的内省函数dir()
查看datasource
方法与属性:
dir(datasource)[:10]
['AbortSQL', 'AddBand', 'AddFieldDomain', 'AddRelationship', 'AdviseRead', 'BeginAsyncReader', 'BuildOverviews', 'ClearStatistics', 'Close', 'CommitTransaction']
driver = ogr.GetDriverByName('ESRI Shapefile')
该函数的原型:
Open(self, char name, int update = 0) -> DataSource
根据数据驱动 driver
的 Open()
方法返回一个数据源对象,update
0为只读,1为可写:
import sys
inshp = '/data/gdata/GSHHS_c.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
datasource = driver.Open(inshp,0)
if datasource is None:
print('could not open')
else:
print('QED.')
QED.
这一节看一下如何使用 ogr
模块来获取图层的信息。
图层的概念
在GIS中,Layer(图层)由同种要素(Feature)(如点、线、多边形等)组在一起的“层”。Layer的描述与ESRI在ArcGIS中定义的模型要素类(Feature Class)一致,其对应的是要素数据集。
在《Modeling our World》中,这两个概念均有明确的阐述:
- 要素类:具有相同几何形状的要素的集合:点、线或多边形。我们最常见的两种要素类是简单要素类和拓扑要素类。简单要素类是指没有任何拓扑关系的点、线、多边形或注记。也就是说,一个要素类内的点与另一要素类中线要素的终点可以是一致的,但它们是不同的,其特点是可独立编辑。拓扑要素类局限在一定的图形范围内,它是一个由完整拓扑单元组成的一组要素类限定的对象。
- 要素数据集(要素集):具有相同坐标系统的要素类的集合。我们可以选择在要素集的内部或外部组织简单要素类,但拓扑要素类只能在要素集内部组织,以确保它们具有相同的坐标系统。
Python中操作Layer方法
Layer方法集合如下:
from osgeo import ogr
inshp='/data/gdata/GSHHS_c.shp'
datasource = ogr.Open(inshp)
layer = datasource.GetLayer(0)
同样看下 layer
对象支持的方法:
dir(layer)
['AlterFieldDefn', 'AlterGeomFieldDefn', 'Clip', 'CommitTransaction', 'CreateFeature', 'CreateField', 'CreateFieldFromArrowSchema', 'CreateFieldFromPyArrowSchema', 'CreateFields', 'CreateGeomField', 'DeleteFeature', 'DeleteField', 'Dereference', 'Erase', 'ExportArrowArrayStreamPyCapsule', 'FindFieldIndex', 'GetArrowArrayStreamInterface', 'GetArrowStream', 'GetArrowStreamAsNumPy', 'GetArrowStreamAsPyArrow', 'GetDataset', 'GetDescription', 'GetExtent', 'GetExtent3D', 'GetFIDColumn', 'GetFeature', 'GetFeatureCount', 'GetFeaturesRead', 'GetGeomType', 'GetGeometryColumn', 'GetGeometryTypes', 'GetLayerDefn', 'GetMetadata', 'GetMetadataDomainList', 'GetMetadataItem', 'GetMetadata_Dict', 'GetMetadata_List', 'GetName', 'GetNextFeature', 'GetRefCount', 'GetSpatialFilter', 'GetSpatialRef', 'GetStyleTable', 'GetSupportedSRSList', 'Identity', 'Intersection', 'IsArrowSchemaSupported', 'IsPyArrowSchemaSupported', 'Reference', 'Rename', 'ReorderField', 'ReorderFields', 'ResetReading', 'RollbackTransaction', 'SetActiveSRS', 'SetAttributeFilter', 'SetDescription', 'SetFeature', 'SetIgnoredFields', 'SetMetadata', 'SetMetadataItem', 'SetNextByIndex', 'SetSpatialFilter', 'SetSpatialFilterRect', 'SetStyleTable', 'StartTransaction', 'SymDifference', 'SyncToDisk', 'TestCapability', 'Union', 'Update', 'UpdateFeature', 'UpsertFeature', 'WriteArrow', 'WriteArrowBatch', 'WriteArrowSchemaAndArrowArrayCapsule', 'WriteArrowStreamCapsule', 'WritePyArrow', '__arrow_c_stream__', '__bool__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__iter__', '__le__', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__nonzero__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_parent_ds', 'schema', 'this', 'thisown']
注意:GetLayer
的参数是从0开始的。
对于Shapefile而言,它只有一个图层,一般情况下这个参数都是0,如果空着,缺省情况下也是0。
查看图层要素:
layer.GetFeatureCount()
1765
还可以查看图层范围:
layer.GetExtent()
(-180.0, 180.00000000000023, -89.99999999999994, 83.53036100000006)
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
这就是整个表的结构。 当然类型是枚举类型,若要查看数据类型名,可对应地建立一个数据类型字典。
获取要素(Feature)信息
地理要素是地图的地理内容,包括自然地理要素与社会经济要素。其中,自然地理要素是指地球表面自然形态所包含的要素,如地貌、水系、植被和土壤等;社会经济要素是指人类在生产活动中改造自然界所形成的要素,如居民地、道路网、通讯设备、工农业设施、经济文化和行政标志等。
如果单从计算机角度看,要素就是一些几何形状,其实这种观点是欠妥的。几何形状,具体来说包括点、线、面、多边形、弧段、向量、控制点等多种类型,其实就是要素的抽象模型。
在Shapefile中,要素模型由点、线、面三种类型构成。要素类自带属性信息,一般对应属性表中的一行。
Layer的核心是针对Feature的操作。GetFeatureCount()
和 GetFeature()
或 GetNextFeature()
可获取层内所有的要素;GetFeaturesRead()
可查询目前已读取了多少条Feature
;生成器可读取要素值;若要再次访问获取过的要素,可访问ResetReading()
从头读取。
from osgeo import ogr
inshp ='/data/gdata/GSHHS_c.shp'
datasource = ogr.Open(inshp)
layer = datasource.GetLayer(0)
feature = layer.GetFeature(0)
dir(feature)
['Clone', 'Dereference', 'Destroy', 'DumpReadable', 'DumpReadableAsString', 'Equal', 'ExportToJson', 'FillUnsetWithDefault', 'GetDefnRef', 'GetFID', 'GetField', 'GetFieldAsBinary', 'GetFieldAsDateTime', 'GetFieldAsDouble', 'GetFieldAsDoubleList', 'GetFieldAsISO8601DateTime', 'GetFieldAsInteger', 'GetFieldAsInteger64', 'GetFieldAsInteger64List', 'GetFieldAsIntegerList', 'GetFieldAsString', 'GetFieldAsStringList', 'GetFieldCount', 'GetFieldDefnRef', 'GetFieldIndex', 'GetFieldType', 'GetGeomFieldCount', 'GetGeomFieldDefnRef', 'GetGeomFieldIndex', 'GetGeomFieldRef', 'GetGeometryRef', 'GetNativeData', 'GetNativeMediaType', 'GetStyleString', 'IsFieldNull', 'IsFieldSet', 'IsFieldSetAndNotNull', 'Reference', 'SetFID', 'SetField', 'SetFieldBinary', 'SetFieldBinaryFromHexString', 'SetFieldDoubleList', 'SetFieldInteger64', 'SetFieldInteger64List', 'SetFieldIntegerList', 'SetFieldNull', 'SetFieldString', 'SetFieldStringList', 'SetFrom', 'SetFromWithMap', 'SetGeomField', 'SetGeomFieldDirectly', 'SetGeometry', 'SetGeometryDirectly', 'SetNativeData', 'SetNativeMediaType', 'SetStyleString', 'UnsetField', 'Validate', '_SetField2', '_SetFieldBinary', '__class__', '__cmp__', '__copy__', '__del__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__getitem__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setitem__', '__sizeof__', '__str__', '__subclasshook__', '__swig_destroy__', '__weakref__', '_add_geom_ref', '_getfieldindex', '_invalidate_geom_refs', 'geometry', 'items', 'keys', 'this', 'thisown']
其中名称中含有 Field
是属性操作, 含有 Geometry
是图形操作。
按顺序读取feature,循环遍历所有的feature:
layer.ResetReading()
feature = layer.GetNextFeature()
while feature:
feature = layer.GetNextFeature()
layer.ResetReading()
feat=layer.GetFeature(0)
feat.keys()
['id', 'level', 'source', 'parent_id', 'sibling_id', 'area', 'Shape_Leng', 'Shape_Area']
feat.GetField('level')
1
以上是使用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
geom = feat.GetGeometryRef()
geom.GetGeometryName()
'POLYGON'
geom.GetGeometryCount()
326
geom.GetPointCount()
0
geom.ExportToWkt()[:80]
'POLYGON ((107.544833 76.914028,106.533222 76.499556,111.093417 76.7653890000001,'
这是Geometry的主要操作方法。这些操作通过循环已经可以把图形画出来。
需要注意的是如果几何类型是多边形,那么就有两层 Geometry ;获取多边形后,还需要获取构成多边形边缘的折线,才可以读取构成多边形的点。
如果类型是单线或者点,就可以直接用GetX()
, GetY()
,而不用再多花一步GetGeometryRef()
来获取子形状了。
polygon=geom.GetGeometryRef(0)
polygon.GetGeometryName()
'LINEARRING'
polygon.GetGeometryCount()
0
polygon.GetPointCount()
1004
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'
需要注意的是这些点都是地理意义上的坐标点,而不是数据意义的坐标。单位是根据空间参考决定的(此数据的单位应该是经纬度,而不是投影坐标的单位米)。获取了这些坐标后,不需要像栅格数据那样用四周边界点来计算地理坐标,矢量数据可直接铺展到一个实际平面已经定义好的投影或者地理坐标系上。
栅格数据有像元长宽,并且有数据意义上的行列坐标,而矢量在没有贴到某个实际存在的物体表面-如计算机屏幕或者打印的纸张之前,是无法表示数据意义上的长宽,所以比例尺对于矢量数据来说,没有绘制操作也就无根本意义。
from osgeo import ogr
inshp ='/data/gdata/GSHHS_c.shp'
datasource = ogr.Open(inshp)
layer = datasource.GetLayer(0)
feature = layer.GetFeature(0)
feature.Destroy()
关闭数据源,相当于文件系统操作中的关闭文件。
from osgeo import ogr
inshp ='/data/gdata/GSHHS_c.shp'
datasource = ogr.Open(inshp)
datasource.Destroy()
from osgeo import ogr
inshp = '/data/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)[:10]
['AddGeometry', 'AddGeometryDirectly', 'AddPoint', 'AddPointM', 'AddPointZM', 'AddPoint_2D', 'Area', 'AssignSpatialReference', 'Boundary', 'Buffer']
这里获取的外包矩形只针对于这个要素,而不是整个图层。