空间信息处理常用工具简介 Part3

大家好，欢迎来到IT知识分享网。

面向空间信息处理的工具

WKT格式

WKT(Well-known text)是一种文本标记语言，用于表示矢量几何对象、 空间参考信息(spatial reference) 。它的二进制表示方式，亦即 WKB (well-known binary)，则以更紧凑、更便于计算机处理的形式来传输和存储相同的信息，虽然这种方式是人类不可读的。

表示几何对象

存储空间参考

地理坐标系

GEOGCS[“WGS 84”, DATUM[“WGS_1984”, SPHEROID[“WGS 84”, , 298., AUTHORITY[“EPSG”, “7030”]], AUTHORITY[“EPSG”, “6326”]], PRIMEM[“Greenwich”, 0, AUTHORITY[“EPSG”, “8901”]], UNIT[“degree”, 0.09433, AUTHORITY[“EPSG”, “9122”]], AUTHORITY[“EPSG”, “4326”]]

投影坐标系

PROJCS[“WGS 84 / UTM zone 50N”, GEOGCS[“WGS 84”, DATUM[“WGS_1984”, SPHEROID[“WGS 84”, , 298., AUTHORITY[“EPSG”, “7030”]], AUTHORITY[“EPSG”, “6326”]], PRIMEM[“Greenwich”, 0, AUTHORITY[“EPSG”, “8901”]], UNIT[“degree”, 0.09433, AUTHORITY[“EPSG”, “9122”]], AUTHORITY[“EPSG”, “4326”]], PROJECTION[“Transverse_Mercator”], PARAMETER[“latitude_of_origin”, 0], PARAMETER[“central_meridian”, 117], PARAMETER[“scale_factor”, 0.9996], PARAMETER[“false_easting”, ], PARAMETER[“false_northing”, 0], UNIT[“metre”, 1, AUTHORITY[“EPSG”, “9001”]], AUTHORITY[“EPSG”, “32650”]]

shapely库

读写WKT格式的最佳的方式是使用Shapely库。Shapely包含一个可以载入和导出WKT格式数据的模块。可以通过pip在线安装，也可以从Unofficial Windows Binaries for Python Extension Packages（https://www.lfd.uci.edu/~gohlke/pythonlibs/）上将其whl文件下载后安装。由于在安装GDAL时已经将Shapely安装过了，所以下面直接调用该库。

# 载入多边形，将其由wkt变为几何对象（多边形），然后计算多边形的面积，最后再将该对象转换为原先的WKT格式 import shapely.wkt wktPoly = "POLYGON((0 0,4 0,4 4,0 4,0 0),(1 1, 2 1, 2 2, 1 2,1 1))" poly = shapely.wkt.loads(wktPoly) print(poly.area) print(poly.wkt) 

输出：

15.0 POLYGON ((0 0, 4 0, 4 4, 0 4, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))

Shapely库还能够处理WKB（Well-known binary， WKB）格式的数据。 Shapely库载入WKB数据使用的是和载入WKT文本类似的wkb模块， 它支持通过调用对象的wkb属性来转换几何体。

在下面的示例中， 将从Shapefile数据集中打开一个包含简单多边形的Shapefile文件-polygon.shp， 然后调用GetLayer()方法， 获取第一个（唯一的） 图层， 然后将其中的特征（Feature）转换为WKT格式。之后就可以使用包含导出结果的wkt变量获取WKT字符串了。 接下来将它重新导入OGR中， 并获取它的多边形边框。

需要注意的是OGR模块须单独读取每个特征然后才能将其导出， 因为ExporttoWkt()方法只能处理特征对象。

from osgeo import ogr shape = ogr.Open("polygon/polygon.shp") layer = shape.GetLayer() feature = layer.GetNextFeature() geom = feature.GetGeometryRef() print('Length:',geom.GetGeometryRef(0).Length()) wkt = geom.ExportToWkt() print('WKT:') print(wkt) poly = ogr.CreateGeometryFromWkt(wkt) print('Envelope:') print(poly.GetEnvelope()) 

输出：

Length: 115.964 WKT: POLYGON ((-99.64 51.51,-75.07 46.08,-75.07 46.08,-75.07 46.08,-76.25 23.69,-76.25 23.69,-76.25 23.69,-114.2 26.37,-114.2 26.37,-99.64 51.51)) Envelope: (-114.2, -75.07, 23.69, 51.51) 

Shapely模块在前面讨论WKT文本导入和导出时已经用到了。 但其实它的主要用途是通用几何库，功能聚焦于几何操作方面。

下面程序首先定义一个简单的WKT多边形， 其次将其导入到Shapely中并测量多边形面积， 接着对该多边形执行缓冲区分析操作， 之后会计算缓冲区的面积。

from shapely import wkt wktPoly = "POLYGON((0 0,4 0,4 4,0 4,0 0))" poly = wkt.loads(wktPoly) print(poly.area) buf = poly.buffer(5.0) print(buf.area) 

输出：

16.0 174.848 

还可以计算该多边形与其缓冲区之间的差异。

print(buf.difference(poly).area) 

输出：

158.848

荟萃知识，滋养你我。