OGR + Python 要素的描述
目录
需要解释解释Layer了。这里的Layer指是一个由同种要素(Feature)组合在一起的”层”。相当于在ESRI定义的模型中的要素类(FeatureClass)，也可能相当于一个要素数据集(FeatureDataSet)。总之是要素的集合。 关于要素数据集，Modeling our world这样解释：
要素数据集（要素集）是具有相同坐标系统的要素类的集合。我们可以选择在要素集的内部或外部 组织简单要素类，但拓扑要素类只能在要素集内部组织，以确保它们具有相同的坐标系统。
对于要素类这样解释：
要素类是具有相同几何形状的要素的集合：点、线或多边形。我们最关心的两种要素类是简单要素 类和拓扑要素类。
简单要素类包括没有任何拓扑关系的点、线、多边形或注记。也就是说，一个要素类内的点与另一 要素类中的线的终点可以是一致的，但它们是不同的。这些要素可以彼此独立地编辑。
拓扑要素类局限在一定的图形范围内，它是一个由完整拓扑单元组成的一组要素类限定的对象。
可以看到，Layer核心操作都是针对的是要素的操作(Feature)，并且通过GetFeatureCount和GetFeature或者GetNextFeature可以获取层内所有的要素，通过GetFeaturesRead可以知道现在已经读取了多少条Feature了(获取一次值少1，如果要从头开始，用ResetReading)。那么什么是要素呢？
其实，要素就是一些形状。具体得分会分成点，线，多边形，弧段，向量，控制点等等非常多种。不管那么多，先把它理解成点线面吧。要素类一般要带一个属性表 （要素+属性才组成GIS，所谓的GIS空间查询就是要素和属性之间的对应，根据条件选择的要素对应到属性，选择的属性对应到要素）。一个要素一般对应属 性表中的一行。
1. OGR简介
OGR是一个读取和处理GSI矢量数据的库。这个库可以读取和处理多种流行的矢量数据(关于能处理什么数据，请看这里，注意一些数据只能读取不能写入)。
2. 安装
OGR是GDAL库的一个部分，只要你安装了GDAL库，就已经拥有了OGR库。至于GDAL如何安装，请察看这里。
3. 快速开始
3.1. 打开一个矢量数据源
第一步还是看看如何打开一 个数据。矢量数据常用的数据集一种是ESRI的shapefile，一种是MapInfo的TAB。它们都分别包含了数个文件，单一得拷贝某个文件是无法 打开的。但是在OGR中要打开它们只要打开其中最主要的文件就可以。在shapefile里，这个文件是以.shp结尾的。在TAB里，这个文件是 以.tab结尾的。
下面我们打开一个ShapeFile文件。这个文件是标准教程数据中的Country数据。它其实就是一幅最常见的世界地图。
3.2. 数据源, DataSource
1 >>> import ogr
2 >>> dir(ogr)
3 ['BuildPolygonFromEdges', 'CreateGeometryFromGML', 'CreateGeometryFromWkb', 'Cre
4 ateGeometryFromWkt', 'DataSource', 'Driver', 'Feature', 'FeatureDefn', 'FieldDef
5 n', 'Geometry', 'GetDriver', 'GetDriverByName', 'GetDriverCount', 'GetOpenDS', '
6 GetOpenDSCount', 'Layer', 'ODrCCreateDataSource', 'ODrCDeleteDataSource', 'ODsCC
7 reateLayer', 'ODsCDeleteLayer', 'OFTBinary', 'OFTInteger', 'OFTIntegerList', 'OF
8 TReal', 'OFTRealList', 'OFTString', 'OFTStringList', 'OFTWideString', 'OFTWideSt
9 ringList', 'OGRError', 'OJLeft', 'OJRight', 'OJUndefined', 'OLCCreateField', 'OL
10 CDeleteFeature', 'OLCFastFeatureCount', 'OLCFastGetExtent', 'OLCFastSetNextByInd
11 ex', 'OLCFastSpatialFilter', 'OLCRandomRead', 'OLCRandomWrite', 'OLCSequentialWr
12 ite', 'OLCTransactions', 'Open', 'OpenShared', 'SetGenerate_DB2_V72_BYTE_ORDER',
13 '__builtins__', '__doc__', '__file__', '__name__', '_gdal', 'gdal', 'osr', 'ptr
14 add', 'ptrcast', 'ptrcreate', 'ptrfree', 'ptrmap', 'ptrptrcreate', 'ptrptrset',
15 'ptrptrvalue', 'ptrset', 'ptrvalue', 'sys', 'types', 'wkb25Bit', 'wkbGeometryCol
16 lection', 'wkbGeometryCollection25D', 'wkbLineString', 'wkbLineString25D', 'wkbL
17 inearRing', 'wkbMultiLineString', 'wkbMultiLineString25D', 'wkbMultiPoint', 'wkb
18 MultiPoint25D', 'wkbMultiPolygon', 'wkbMultiPolygon25D', 'wkbNDR', 'wkbNone', 'w
19 kbPoint', 'wkbPoint25D', 'wkbPolygon', 'wkbPolygon25D', 'wkbUnknown', 'wkbXDR']
20 >>> datasource = ogr.Open("e:/gisdata/data/country.shp")
这样就打开了一个数据源(DataSource)。
让我们看看这个dataSource里有什么吧
1 >>> dir(datasource) 2 ['CopyLayer', 'CreateLayer', 'DeleteLayer', 'Dereference', 'Destroy', 'ExecuteSQ 3 L', 'GetDriver', 'GetLayer', 'GetLayerByName', 'GetLayerCount', 'GetName', 'GetR 4 efCount', 'GetSummaryRefCount', 'Reference', 'Release', 'ReleaseResultSet', 'Tes 5 tCapability', '__doc__', '__getitem__', '__init__', '__len__', '__module__', '_o 6 ', 'driver_o'] 7 >>>
还是那个规矩，在快速开始的时候我们只看Get类型的函数。这里看来可以操作的概念的只有Layer了。
3.3. 层, Layer
需要解释解释Layer了。这里的Layer指是一个由同种要素(Feature)组合在一起的”层”。相当于在ESRI定义的模型中的要素类(FeatureClass)，也可能相当于一个要素数据集(FeatureDataSet)。总之是要素的集合。 关于要素数据集，Modeling our world这样解释：
要素数据集（要素集）是具有相同坐标系统的要素类的集合。我们可以选择在要素集的内部或外部 组织简单要素类，但拓扑要素类只能在要素集内部组织，以确保它们具有相同的坐标系统。
对于要素类这样解释：
要素类是具有相同几何形状的要素的集合：点、线或多边形。我们最关心的两种要素类是简单要素 类和拓扑要素类。
简单要素类包括没有任何拓扑关系的点、线、多边形或注记。也就是说，一个要素类内的点与另一 要素类中的线的终点可以是一致的，但它们是不同的。这些要素可以彼此独立地编辑。
拓扑要素类局限在一定的图形范围内，它是一个由完整拓扑单元组成的一组要素类限定的对象。
下面是Layer的方法集合
1 >>> layer = datasource.GetLayer(0) 2 >>> dir(layer) 3 ['CommitTransaction', 'CreateFeature', 'CreateField', 'DeleteFeature', 'Derefere 4 nce', 'GetExtent', 'GetFeature', 'GetFeatureCount', 'GetFeaturesRead', 'GetLayer 5 Defn', 'GetName', 'GetNextFeature', 'GetRefCount', 'GetSpatialFilter', 'GetSpati 6 alRef', 'Reference', 'ResetReading', 'RollbackTransaction', 'SetAttributeFilter' 7 , 'SetFeature', 'SetNextByIndex', 'SetSpatialFilter', 'SetSpatialFilterRect', 'S 8 tartTransaction', 'SyncToDisk', 'TestCapability', '__doc__', '__init__', '__len_ 9 _', '__module__', '_o'] 10 >>>
可以看到，Layer核心操作都是针对的是要素的操作(Feature)，并且通过GetFeatureCount和GetFeature或者GetNextFeature可以获取层内所有的要素，通过GetFeaturesRead可以知道现在已经读取了多少条Feature了(获取一次值少1，如果要从头开始，用ResetReading)。那么什么是要素呢？
3.4. 要素，Feature
其实，要素就是一些形状。 具体得分会分成点，线，多边形，弧段，向量，控制点等等非常多种。不管那么多，先把它理解成点线面吧。要素类一般要带一个属性表（要素+属性才组成 GIS，所谓的GIS空间查询就是要素和属性之间的对应，根据条件选择的要素对应到属性，选择的属性对应到要素）。一个要素一般对应属性表中的一行。
抓一个feature出来看看……
1 >>> feature = layer.GetFeature(0) 2 >>> dir(feature) 3 ['Clone', 'Destroy', 'DumpReadable', 'Equal', 'GetDefnRef', 'GetFID', 'GetField' 4 , 'GetFieldAsDouble', 'GetFieldAsInteger', 'GetFieldAsString', 'GetFieldCount', 5 'GetFieldDefnRef', 'GetFieldIndex', 'GetGeometryRef', 'GetStyleString', 'IsField 6 Set', 'SetFID', 'SetField', 'SetFrom', 'SetGeometry', 'SetGeometryDirectly', 'Se 7 tStyleString', 'UnsetField', '__cmp__', '__copy__', '__del__', '__doc__', '__get 8 attr__', '__init__', '__module__', '__setattr__', '_o', 'thisown']
里面可以看到有关Field和Geometry的操作，没错，这就是Feature的本质。Field有关的是属性的操作。Geometry有关的是形状的操作。
1 >>> for i in range(feature.GetFieldCount()): 2 ... print feature.GetField(i) 3 ... 4 AA 5 ABW 6 Aruba 7 Netherlands 8 67074 9 182.926 10 70.628 11 Florin 12 AWG 13 N 14 1
这是这个feature的所有属性值 如果你要看整个表的结构，各个字段的名称等等信息，那就不能在feature里看了。要看得到layer的附加信息里看。
1 >>> layerdef = layer.GetLayerDefn() 2 >>> for i in range(layerdef.GetFieldCount()): 3 ... defn = layerdef.GetFieldDefn(i) 4 ... print defn.GetName(),defn.GetWidth(),defn.GetType(),defn.GetPrecision() 5 ... 6 FIPS_CNTRY 2 4 0 7 GMI_CNTRY 3 4 0 8 CNTRY_NAME 40 4 0 9 SOVEREIGN 40 4 0 10 POP_CNTRY 10 0 0 11 SQKM_CNTRY 12 2 3 12 SQMI_CNTRY 12 2 3 13 CURR_TYPE 16 4 0 14 CURR_CODE 4 4 0 15 LANDLOCKED 1 4 0 16 COLOR_MAP 1 4 0 17 >>>
这就是整个表结构。当然类型是枚举，要看数据类型名，自己去建立一个数据类型的字典去对应，我这里就不写了吧。
3.5. 形状，Geometry
1 >>> geom = feature.GetGeometryRef() 2 >>> dir(geom) 3 ['AddGeometry', 'AddGeometryDirectly', 'AddPoint', 'AddPoint_2D', 'AssignSpatial 4 Reference', 'Buffer', 'Centroid', 'Clone', 'CloseRings', 'Contains', 'ConvexHull 5 ', 'Crosses', 'Destroy', 'Difference', 'Disjoint', 'Empty', 'Equal', 'ExportToGM 6 L', 'ExportToWkb', 'ExportToWkt', 'FlattenTo2D', 'GetArea', 'GetBoundary', 'GetC 7 oordinateDimension', 'GetDimension', 'GetEnvelope', 'GetGeometryCount', 'GetGeom 8 etryName', 'GetGeometryRef', 'GetGeometryType', 'GetPointCount', 'GetSpatialRefe 9 rence', 'GetX', 'GetY', 'GetZ', 'Intersect', 'Intersection', 'Overlaps', 'SetCoo 10 rdinateDimension', 'SetPoint', 'SetPoint_2D', 'SymmetricDifference', 'Touches', 11 'Transform', 'TransformTo', 'Union', 'Within', 'WkbSize', '__copy__', '__del__', 12 '__doc__', '__init__', '__module__', '__str__', '_o', 'thisown'] 13 >>> geom.GetGeometryName() 14 'POLYGON' 15 >>> geom.GetGeometryCount() 16 1 17 >>> geom.GetPointCount() 18 0 19 >>> polygon = geom.GetGeometryRef(0) 20 >>> polygon.GetGeometryName() 21 'LINEARRING' 22 >>> polygon.GetGeometryCount() 23 0 24 >>> polygon.GetPointCount() 25 11 26 >>> polygon.GetX(0) 27 -69.882232666015625 28 >>> polygon.GetY(0) 29 12.411109924316406 30 >>> polygon.GetZ(0) 31 0.0 32 >>>
这是有关Geometry的主要操作(你 还可以试试ExportToWkt这个方法，可以得到整个多边形的Wkt表示，可以看到所有的点坐标)。这些操作通过循环已经可以把Aruba的国家轮廓
画出来了(把数据贴在这里应该不犯法吧～，这个可是标准Demo数据哦)。需要注意的是Geometry的类型如果是多边形，那么就有两层 Geometry，获取多边形后，还需要获取构成多边形边缘的折线，才可以获得构成多边形的点。如果类型是单线或者点，那么就可以直接用 GetX，GetY，而不用再多花一步GetGeometryRef获取子形状了。
需 要注意的是这些点都是地理意义上的坐标，而不是数据意义坐标。单位是根据空间参考决定的(这个数据所用的应该是经纬度，而不是投影坐标单位米)。获取了这 些坐标，不需要像栅格数据那样根据四周边界点来计算地理坐标(栅格数据有象元长宽，并且有数据意义上的行列坐标，而矢量在没有贴到某个实际存在的物体表面 －如计算机屏幕或者打印的纸张－前实际上无法表示数据意义上的长宽，所以比例尺对于矢量数据来说，没有绘制操作就根本没有意义)，应该直接铺展到一个实际 平面上已经定义好的投影或者地理坐标系上。
矢量和栅格数据集一样，有空间参考信息的。而且空间参考信息直接渗入每个Feature中。 因为ShapeFile的坐标系统和仿射参数是分开的。坐标系统是用.prj文件附加在shp文件的同目录下获取的。但是仿射参数是直接放在 ShapeFile文件内部的，只不过它和GeoTiff的tfw不一样，它的表示是一个外包矩形。记录整个layer图层的地理范围。因为我找不到那个
prj文件了，所以这里的空间参考也没法读出。但是地理范围还是可以的。
1 >>> layer.GetSpatialRef() 2 >>> layer.GetExtent() 3 (-180.0, 180.0, -90.0, 83.62359619140625) 4 >>>
除了整个图层有地理范围，每个feature也有地理范围。不过这个地理范围是放在GeometryDef里的。
1 >>> geom.GetEnvelope() 2 (-70.059661865234375, -69.874862670898437, 12.411109924316406, 12.627776145935059) 3 >>> geom.GetSpatialReference() 4 >>>
这时获取的外包矩形只是对于这个要素的，而不是整个图层的。而这里的坐标系统也没法获取。获取出来也应当与Layer的一样。
我们还看到Geometry里面有很多奇怪的方法，如Buffer，Contains，Crosses，Disjoint，Overlaps，Union等等。这些都是一些几何关系运算，已经涉及空间分析的范畴了，关于几何关系运算可以单独开一篇来讨论了(几何关系运算一般需要另外一个库GEOS的支持，不然OGR无法正常运算，这个话题留到后面去慢慢研究吧)。
有一个方法还是不错的。那就是GetArea()这个方法可以获得整个多边形包围的面积。这个方法我竟然在C++的文档中没有发现，而在Python中发现了(哦～～上帝爱Python～～)。
1 >>> geom.GetArea() 2 0.016598101236013463 3 >>>
不过这个是用“度”为单位的。已经失去它作为面积的意义了。如果要准确的面积，可以根据空间参考来进行从度到米的转换，然后算出平方米。不过我的prj投影文件不见了，没有办法了。
4. 总结
我们可以看到，整个 ogr的结构很明确，每个部分都很标准。首先是数据源。打开数据源里面有一系列的图层(不管它是矢量数据集还是要素类)，图层中就是许多相同空间参考下的 表示地物的要素(要素是地物的抽象，至于具体的抽象表达形式，可以找本导论看看)。要素包括了相同空间参考下的几何形状和关联属性。几何形状又由几个子形 状或者一系列点组成，这些点正好代表它们各自的地理位置。这样整个矢量模型就被完整得拆分开了。而创建矢量数据的过程就是上面这个过程的反过程。
转载自：https://blog.csdn.net/yht_roy/article/details/38515325
