空间数据模型

数据模型

由Shapely提供的几何对象的基本类型包括点, 曲线,和曲面。 每一种都与三套平面的点 (可能无限)相联系。 一个要素的内部, 边界, 外部是相互排斥的,并且它们的并集与整个平面重叠。

  • 点的内部有一个点,边界没有点,点的外部包括其它所有点. 点的拓扑维数是0。
  • 曲线有沿其长度的无穷多个点组成的内部集 (想象一个点在空间拖动), 包括2个端点的边界集, 和包括其它所有点的外部集. 曲线的拓扑维数是1。
  • 曲面有一个包含无穷个点内部集 (想象一条曲线在空间上拖动以覆盖一个区域), 有一个由一个或多个曲面组成的边界, 和包含所有其它点的外部集。包括那些所有在表面可能存在的洞.。曲面的拓扑维数是2。

这可能有一点深奥, 但这将有助于澄清Shapely的空间谓词的含义, 这本手册将和理论一样深刻。 点集理论的推论,其中包括一些 用于不同类时的体现出的缺陷,将会后面进行讨论。

点是由点类实现的; 曲线是由直线和环线类实现的; 曲面是由多边形类实现的。 Shapely不执行平滑的 ( 即有连续的切线)曲线。 所有的曲线必须近似线性样条曲线。 所有的圆斑必须近似线性样条范围内的地区。

点的集合由一多样点类来实现, 曲线的集合由一多样线类来实现, 曲面的集合由一多样多边形来实现. 这些集合在计算上作用不是很大, 但对某些种类的功能建模非常有用。例如,一个Y形线 可以很好地由多样线建模为一个整体。

标准的数据模型对几何对象的特定类型 有特定的附加约束, 这将在接下来的章节里讨论。

一些术语

  • 坐标, 在已经定义好的精确度模型中可以被精确地描绘的空间中的一点。
  • 精确计算, 在操作过程中支持所有数字的数值计算,通常要用到大量复杂的算法。
  • 节点, 相同或不同的几何图形内两条线相交的点。这点不必用坐标表示出来,因为大部分的相交点的输出计算比输入计算需要更加严格的精确度。
  • 节点的计算, 计算出在一个或多个几何体相交之处形成的节点的过程。
  • 非坐标, 不可被描绘成一个坐标的点。
  • 数值稳定性, 算法的稳定性取决于它的输出结果中误差的最大范围。如果算法的误差范围是很小的,那么可以认为它是稳定的。
  • 点, R3中任意一点,一般来说,可以无限地描绘出来。
  • 真交点, 两条线段相交所形成的唯一交点,并且是在这两条线段上的点。
  • 健壮计算, 一种对于所有的输入都必定会输出正确答案的数值计算方法,通常需要有特别处理舍入误差而设的算法。
  • SFS, 开放式地理信息系统协会OGC的简单要素实现规范
  • 分辨率单位, 在已经定义好的精确度模型中的最小的可以表示出来的距离。
  • 顶点, 几何体中的边角点,这些点的坐标都可以确定的话,就可以确定一个几何体的位置。

Shapely的基本技术特征

  1. 得到二维形状的联合分析,交叉分析,或者差异分析的结果;
  2. 决定二维形状是否相交,接触,或者一个中包含另一个。

以及更多。

关系

空间数据模型有一组几何对象之间关系的自然语言 – 包含,相交,叠置,相切等– 和一个理论框架用于理解它们 使用它们的组成点集的相互交叉的3x3矩阵。

运算符

承接JTS技术规格, 本手册将区别构造(缓冲,凸壳)和 理论集操作(相交,联合等)。

坐标系

虽然地球是不平的,但是有很多分析类的问题能够通过 试验的和真正的算法,将地球特征转化到直角平面, 接着将结果转为地理坐标来解决。这种做法与传统的精确纸质地图一样古老。

Shapely不支持坐标系统转换。 所有2个或更多特征的操作假定特征在同一个直角平面。

空间操作与计算

GEOS主要支持的操作和计算:

空间关系计算,主要支持的计算

OGC Filter定义了一些xml标签,从而实现了查询Feature集合的接口。 程序可以将它转换为查询语言,检索Feature集合后返回有效的Feature子集合。

  • 相等(Equals)几何形状拓扑上相等。

  • 脱节(Disjoint)几何形状没有共有的点。

  • 相交(Intersects)几何形状至少有一个共有点(区别于脱节)

  • 接触(Touches)几何形状有至少一个公共的边界点,但是没有内部点。

  • 交叉(Crosses)几何形状共享一些但不是所有的内部点。

  • 内含(Within)几何形状A的线都在几何形状B内部。

  • 包含(Contains)几何形状B的线都在几何形状A内部(区别于内含)

  • 重叠(Overlaps)几何形状共享一部分但不是所有的公共点,而且相交处有他们自己相同的区域。

一些情况下,确切定义有些微妙。 你可以参考JTS技术说明书来完全确定在任何情况下的返回情况。

以上的运算返回的都是 true 或者 false

另外一种是空间叠加分析操作。主要有下面几个操作:

  • 缓冲区分析(Buffer)包含所有的点在一个指定距离内的多边形和多多边形。

  • 凸壳分析(ConvexHull)包含几何形体的所有点的最小凸壳多边形,(就是外包多边形啦)

  • 交叉分析(Intersection)交叉操作就是多边形AB中所有共同点的集合。

  • 联合分析(Union)AB的联合操作就是AB所有点的集合。

  • 差异分析(Difference)AB形状的差异分析就是A里有B里没有的所有点的集合。

  • 对称差异分析(SymDifference)AB形状的对称差异分析就是位于A中或者B中但不同时在AB中的所有点的集合

另外还支持多边形化,连接有向线段,压出节点等等操作。

快速上手

用法

Shapely版本,GEOS库版本和GEOS C API版本均可通过shapely访问。version,

shapely.geos.geos_version_string, and shapely.geos.geos_capi_version.

In [1]:
from osgeo import ogr
import shapely

print(shapely.__version__)
import shapely.geos

print(shapely.geos.geos_version)
print(shapely.geos.geos_version_string)
1.6.3
(3, 6, 2)
3.6.2-CAPI-1.10.2 4d2925d6

缓冲一个点

In [2]:
from shapely.geometry import Point
point = Point(10, 10)
pt_buf = point.buffer(5)
print(pt_buf.wkt)
POLYGON ((15 10, 14.97592363336098 9.509914298352198, 14.90392640201615 9.02454838991936, 14.78470167866104 8.54857661372769, 14.61939766255643 8.086582838174554, 14.40960632174178 7.643016315870014, 14.15734806151273 7.222148834901992, 13.86505226681369 6.828033579181776, 13.53553390593274 6.464466094067266, 13.17196642081823 6.134947733186318, 12.77785116509802 5.842651938487276, 12.35698368412999 5.590393678258227, 11.91341716182545 5.380602337443569, 11.45142338627232 5.215298321338958, 10.97545161008065 5.096073597983849, 10.49008570164781 5.024076366639017, 10.00000000000001 5, 9.509914298352205 5.024076366639015, 9.024548389919367 5.096073597983846, 8.548576613727697 5.215298321338953, 8.086582838174561 5.380602337443563, 7.643016315870021 5.59039367825822, 7.222148834901997 5.842651938487268, 6.82803357918178 6.134947733186308, 6.464466094067269 6.464466094067256, 6.134947733186321 6.828033579181766, 5.842651938487278 7.222148834901982, 5.590393678258229 7.643016315870005, 5.380602337443569 8.086582838174545, 5.215298321338958 8.548576613727683, 5.096073597983849 9.024548389919353, 5.024076366639016 9.509914298352191, 5 9.999999999999995, 5.024076366639015 10.4900857016478, 5.096073597983847 10.97545161008064, 5.215298321338954 11.45142338627231, 5.380602337443565 11.91341716182545, 5.590393678258224 12.35698368412999, 5.842651938487273 12.77785116509801, 6.134947733186315 13.17196642081823, 6.464466094067261 13.53553390593274, 6.828033579181771 13.86505226681368, 7.222148834901985 14.15734806151272, 7.643016315870007 14.40960632174177, 8.086582838174545 14.61939766255643, 8.548576613727679 14.78470167866104, 9.024548389919348 14.90392640201615, 9.509914298352184 14.97592363336098, 9.999999999999986 15, 10.49008570164779 14.97592363336098, 10.97545161008062 14.90392640201616, 11.45142338627229 14.78470167866105, 11.91341716182543 14.61939766255644, 12.35698368412997 14.40960632174179, 12.77785116509799 14.15734806151274, 13.17196642081821 13.8650522668137, 13.53553390593272 13.53553390593276, 13.86505226681367 13.17196642081825, 14.15734806151271 12.77785116509804, 14.40960632174176 12.35698368413002, 14.61939766255642 11.91341716182548, 14.78470167866103 11.45142338627235, 14.90392640201615 10.97545161008068, 14.97592363336098 10.49008570164784, 15 10.00000000000004, 15 10))

有关Shapely提供的空间操作和谓词的更多示例, 请参阅tests.txt和Predicates.txt。 或参见Point.txt,LineString?.txt等几何API的例子。

与shapfile

In [11]:
# Generated in chapter 1.
datasource = ogr.Open('/gdata/GSHHS_c.shp')
layer = datasource.GetLayer(0)
feature = layer.GetFeature(0)
# geom = feature.GetGeometryRef()
polygon = feature.GetGeometryRef()

# print(polygon.ExportToWkt())

buf = polygon.Buffer(2)
buf.ExportToWkt()[:40]
Out[11]:
'POLYGON ((-10.9650568460992 37.422419294'

Numpy接口

所有的Shapely几何实例都提供了Numpy数组接口:

In [4]:
from numpy import asarray
a = asarray(point)
a.size
a.shape
Out[4]:
(2,)

Numpy数组也可以被调整为Shapely的点和线串(linestrings)。

In [5]:
from shapely.geometry import asLineString
a = asarray([[1.0, 2.0], [3.0, 4.0]])
line = asLineString(a)
print(line.wkt)
LINESTRING (1 2, 3 4)

性能

Shapely 使用 GEOS 库来完成所有的操作。 GEOS 是使用 C++ 来完成的并且得到广泛的应用, 可以相信它已经得到了高度的优化。

shapely.speedups模块包含用C编写的性能增强功能。 当Python在安装过程中访问编译器时,它们会自动安装。

您可以检查加速功能是否与可用属性一起安装。 默认情况下禁用构造函数加速。 你可以启用加速调用enable()。 你也可以使用disable()恢复默认的实现。

In [6]:
from shapely import speedups
speedups.available
    
Out[6]:
True
In [12]:
speedups.enable()
In [ ]:
 

JTS、GEOS与Shapely

Shapely是一个用来处理二维地理空间几何要素的Python包。 它基于GEOS库来获得不同几何类型、操作与谓词的定义与实例。

JTS、GEOS与Shapely关系密切。 简单地说,JTS使用Java语言实现了拓扑操作;GEOS则使用C++语言重新进行了实现; 而Shapely则提供了GEOS的Python接口。

对于 Python 开发来讲,使用 Shapely 是最方便的。 但是理解这些概念,还是需要查看JTS的网站,而不是GEOS。

下面分别对各个进行介绍,重点对JTS进行一些说明。

JTS Topology Suite的理解

JTS(JTS Topology Suite [1] 由加拿大的 Vivid Solutions [2]有限公司开发的一套Java API开放源码。 它提供了一套空间数据操作的核心算法, 并为兼容OGC标准的空间对象模型中进行基础的几何操作提供2D空间谓词API。 JTS是一个用Java语言描述的几何拓扑套件,它遵循OpenGIS的Simple Feature Specification, 封装了2D几何类型和非常多的空间分析操作, 而且包含了不少常见的计算几何算法。 JTS解决了对象与对象之间的拓扑关系的判定和计算, 并提供很多有用的算法来解决对象的面积和长度等问题。

JTS拓扑结构程序组是一个通过使用精确的精度模型和明确的几何算法 来执行空间数据操作的Java应用编程接口。 JTS是可以用来支持空间数据集的确认、处理、综合及质疑的应用程序进行改善。

JTS试图尽可能精确地执行开放式地理信息系统GIS中的简单要素实现规范(SFS)。 在一些情况下,SFS是不清晰的,或者是遗漏了某个规范。 在这种情况下,JTS试图选择一个合理的相容的替代物。 SFS的不同之处与其细节之处都将会在这个说明书中论述到。 有关等级与方法的详细文献将会以Java文件的形式出现,作为原始资料代码。

JTS Topology Suite从根本上而言其实并不是很复杂, 它主要是完成了Java对几何对象、空间拓扑得核心操作算法。 个人感觉如果如果把它简单的认为是一个类似于java.utils.*之类的开发包 可能不能真正的体现它的意义, 实际上它除了集成了Java对几何对象(点、线、面等)的对象管理外, 更大一部分工作是在完成对各种几何对象的bufferanalyze以及空间索引。 它尽可能实现了OpenGIS Simple Features Specification规范 [3],

还有OpenGIS Simple Features Specification就用SFS代替

空间数据模型和方法的定义都与SFS规范要求一致; JTS所提供的几何算法可以支持自定义几何对象以及自定义的几何精度, 当然,方法的返回值也是与自定义的几何对象以及几何精度对应的; JTS对算法的正确性要求永远放在第一位, 而其他的效率之类的问题到是放在了其次的地位, 所以说不管是Geotools也好还是GeoServer也好, 只要是核心采用了JTS技术, 那么就需要自己去搞定效率问题了 ; 当然,整个JTS为了能让更多的人使用,所以结构相当的不错, 借用Crespo小兄弟的一句话就是“相当的优雅”,这又是JTS被广泛采用的缘由之一吧。

对JTS而言所有几何数字运算都依托于相应的精度格式。 而精度格式一般会有这么几种。 Fixed,这种情况下坐标相当于在一个由N个正方形构成的离散固定矩阵上; Floaing,这种情况下计算坐标值必须依靠参数来形成一个正方形构成的矩阵然后进行计算, 所以坐标值这个时候并不是一个固定的值,而且会出现对小数点的四舍五入等带来的误差; Exact,这个嘛,呵呵,反正我现在几乎没有遇到这么高的要求, 这种模式要求坐标点是一个分子式来表示(那是相当的精确呀!)。

JTS设计的目的

JTS的设计目的:

  1. 空间模型和方法定义能够达到与开放GIS简单要素实现规范尽可能相符,并且能够正确地执行。
  2. 应用编程孔API设计在任何情况下都符合Java的规范。例如:

    1. 接入口函数要采用Java的getX 和setX的规范。
    2. 谓词要采用isX的规范。
    3. 方法要以小写字母开头。
  3. JTS函数可以支持用户定义的精确度模型。JTS算法对于该精确度模型仍然保持其健壮性。

  4. 在可能的时候,方法会在之前定义的精确度模型中返回拓扑结构和几何图形正确的结果。
  5. 正确性是最重要的,空间与时间的有效性的重要性是次于正确性的。
  6. JTS的速度足够快,可以适用于生产环境中。
  7. 为了让其他的开发者更好地理解,在JTS中使用的算法和代码应该是条理清晰并且结构优化的。

JTS应用情况

关于GIS的开源世界里如Geotools、Udig等, JTS Topology Suite都得到了大量的应用, 甚至可以说没有JTS Topology Suite的话, Geotools等的实现会很复杂,不对,是相当的复杂。 从这些角度而言,我更愿意把JTS Topology Suite看作一个几何对象中间件。

JTS被广泛地应用在开源GIS软件中, JTS被为GeoTools和基于GeoTools的GeoServer和uDig的底层库。

Geotools2使用了JTS所提供的几何对象模型(点、线等)来实现它的矢量化几何模型。 使用JTS 可以将postgis(或者shapefile等其他格式)地图数据转化为JTS几何对象。 就拿昨天发的文章《终于搞定GeoTools对PostGis的操作》中的代码而言, 当featurereader被定义的同时会有一个JTS几何对象数组同时产生, 通过迭代我们可以取得这个数组中每个几何对象的属性 (当然,我们可以先申明一批JTS几何对象然后将此数组赋予featurewriter, 执行write()方法,将几何对象插入到postgisdatastore对象中)。

GEOS

在官方声明中[4]指出,GEOS 是 JTS 的 C++ 端口(port), 这个声明有些容易混淆。 实际上GEOS (Geometry Engine - Open Source)是 JTS的 C++ 实现(独立于JTS运行)。 它的目的是使用 C++ 完成所有 JTS 功能。 包括所有的 OpenGIS 简单要素模型的 SQL 空间谓词函数与空间运算符,以及一些JTS增强的拓扑函数。

GEOS是一个集合形状的拓扑关系操作实用库(可能这么说不太准确), 简单地说,GEOS就是判断两个几何形状之间关系和 对两个几何形状进行操作以形成新的几何形状的库。

Shapely

Shapely是对流行的GEOS库的一个Python封装。 Python 使用 ctypes 实现了使用 Python 直接调用 GEOS 的功能。

Shapely 1.0忽略了坐标与空间参考系统的处理,投影处理的功能留给特别的应用再处理。

Shapely 需要 Python 2.4以上版本。

已知问题:

GEOS的2.2.3和3.0.0版本有几个上游问题。 GEOS的WKT和WKB写入程序泄漏, 因此访问几何的 wktwkb 属性和一些几何类型的创建相关联的漏洞较小。