>>> from env_helper import info; info()
页面更新时间: 2023-04-15 20:01:58
运行环境:
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-7-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

8.8. 评价几何对象之间的关系

在前面一节中,我们使用了 MBR 来评价几何对象之间的空间关系。 但是这种比较只是一种算法上的优化方法,不能进行完全可信的比较,也就是说结果并不一定是对的。 对于要进行精确结果的数据,可以使用 MBR 对要进行空间关系比较的对象进行快速的筛选, 然后再进一步对筛选剩下的结果进行评价。

总体而言,任何一组的几何对象, 都可以建立起下面的一种(或多种)关系:

  • 如果两个几何的交点产生的尺寸小于两个几何的最大尺寸的值, 并且交点值包括两个几何的内部的点,并且交点值不等于第一或第二几何。

  • 如果两个几何的交点都不为空。

  • 如果第二个几何id完全包含在第一个。

  • 允许根据给定的模式进行任意比较矩阵。

  • 是否是空间一致的(空间上相等)。

  • 两个几何对象的交集是否为空。

  • 如果在的两个几何对象之间的共同点在于它们的边界,而没有内部的点。

  • 如果第一个几何对象完全被第二个包含。

  • 如果两个几何对象的“相交”操作生成的结果的维度与这两个是一样的(没有降维),并且此结果与第一个或第二个都不相同。

一些实用(图示)的例子可能会帮助您更好地了解:

  • Equal(g1, g2)

  • Disjoint(g1, g2)

  • Touches(g1, g2)

  • Within(g1, g2)

  • Overlaps(g1, g2)

  • Crosses(g1, g2)

  • Contains(g1, g2)

SpatiaLite包装相应的GEOS库函数,以评估几何之间的空间关系。

SpatiaLite支持以下SQL函数:

  • Equals(geom1,geom2)

  • Disjoint(geom1,geom2)

  • Touches(geom1,geom2)

  • Within ( geom1 , geom2 )

  • Overlaps ( geom1 , geom2 )

  • Crosses ( geom1 , geom2 )

  • Intersects ( geom1 , geom2 )

  • Contains ( geom1 , geom2 )

  • Relate ( geom1 , geom2 , patternMatrix )

所有这些功能可能会返回以下值之一:

  • 1 如果空间关系是 TRUE

  • 0 如果空间关系是 FALSE

  • -1 如何这里有一些误差 [即一个或两个几何形状为NULL]

relate()功能是非常普遍的,而且需要patternMatrix定义的空间关系进行检查:

patternMatrix是一个9个字符的字符串代表一个3 x 3矩阵,每个字符代表内部之间的交集,和每一个几何外观。

  • T 意味着 TRUE

  • F 意味着 FALSE

  • * 意味着不用在乎

  • 0、1、2表示交点必须存在,且必须具有尺寸

例如:T*T***T** 字符串对应于重叠矩阵,而 FF*FF**** 字符串对应不相交矩阵。

在开始之前,先初始化环境:

>>> import os,shutil,stat
>>> import sqlite3 as sqlite
>>> tmp_db = '/tmp/xx_new_db.sqlite'
>>> if os.path.exists(tmp_db):
>>>     os.remove(tmp_db)
>>> shutil.copy("spalite.db", tmp_db)
>>> os.chmod(tmp_db, stat.S_IRUSR + stat.S_IWUSR)
>>> conn = sqlite.connect(tmp_db)
>>> conn.enable_load_extension(True)
>>> conn.execute('SELECT load_extension("mod_spatialite.so.7")')
>>> cursor = conn.cursor()
>>> sql = '''SELECT Regions.Name, COUNT(*) FROM Towns, Regions
>>>  WHERE Regions.Name IN (
>>>  'VALLED','AOSTA', 'PIEMONTE', 'UMBRIA', 'LOMBARDIA',
>>>  'CALABRIA', 'MOLISE', 'MARCHE', 'BASILICATA') AND
>>>  Within(Towns.geometry, Regions.Geometry)
>>>  GROUP BY Regions.Name;'''
>>> cursor.execute(sql)
>>> [print(rec) for rec in cursor]
---------------------------------------------------------------------------

OperationalError                          Traceback (most recent call last)

Cell In [3], line 7
      1 sql = '''SELECT Regions.Name, COUNT(*) FROM Towns, Regions
      2  WHERE Regions.Name IN (
      3  'VALLED','AOSTA', 'PIEMONTE', 'UMBRIA', 'LOMBARDIA',
      4  'CALABRIA', 'MOLISE', 'MARCHE', 'BASILICATA') AND
      5  Within(Towns.geometry, Regions.Geometry)
      6  GROUP BY Regions.Name;'''
----> 7 cursor.execute(sql)
      8 [print(rec) for rec in cursor]


OperationalError: no such table: Towns

这个查询计算了几个城镇落在区域边界内,我们只想扫描有限的地区列表, 因此我们使用IN(..)条款, Within()函数将测试每个城镇和区域之间存在的空间关系边界, 你可以注意到,IN()子句内有一个奇怪的符号:

在字符串中的“VALLE D”是真的需要[它的意图,这不是一个typo]

这是因为有一个apostrophe名词和SQL语法。

解释这种行为的原因很容易识别:

within()功能意味着一个巨大的计算量(任何其他功能)

NewTowns表包含8000行。

我们要求选择8行从newregions表

所以,这个查询评估有64000次within()关系;这需要很长时间。

我们可以从空间索引的适时使用中获得利润, 以优化此查询,我们将更好地解释这一步, 远离分析空间索引的东西。

>>> sql = '''SELECT t2.Name, t2.Peoples,
>>> Distance(t1.geometry, t2.geometry) AS Distance
>>> FROM Towns AS t1, Towns AS t2
>>> WHERE t1.Name = 'Firenze' AND
>>> Distance(t1.geometry, t2.geometry) < 10000'''
>>> cursor.execute(sql)
>>> [print(rec) for rec in cursor]
---------------------------------------------------------------------------

OperationalError                          Traceback (most recent call last)

Cell In [2], line 6
      1 sql = '''SELECT t2.Name, t2.Peoples,
      2 Distance(t1.geometry, t2.geometry) AS Distance
      3 FROM Towns AS t1, Towns AS t2
      4 WHERE t1.Name = 'Firenze' AND
      5 Distance(t1.geometry, t2.geometry) < 10000'''
----> 6 cursor.execute(sql)
      7 [print(rec) for rec in cursor]


OperationalError: no such table: Towns

SpatiaLite Distance()函数测量两个几何之间的最小距离, 这也是一种有用的空间关系。 这个查询确定了布置在佛罗伦萨镇附近的任何城镇。

如果两个几何体没有被分离,我们假定10Km的范围 它们的相对距离恰好为0。

8.8.1. 几何对象之间的逻辑运算

SpatiaLite支持以下SQL函数:

  • geom3 = Intersection ( geom1 , geom2 )

  • geom3 = Difference ( geom1 , geom2 )

  • geom3 = GUnion ( geom1 , geom2 )

  • geom3 = SymDifference ( geom1 , geom2 )

1.所有函数都返回一个表示布尔运算结果的GEOMETRY,请注意返回的GEOMETRY可能是一个函数。

2.OpenGIS标准定义了Union()函数; 但是在SQLite自己的语法中是一个 *  *,所以SpatiaLite将此函数重命名为GUnion()

我们可以在两个几何上应用各种布尔运算,来获得第三个几何。

下面是一些实例,来增进对集合对象的一些了解。

* g3 = Intersection(g1, g2)
* g3 = Difference(g1, g2)
* g3 = GUnion(g1, g2)
* g3 = SymDifference(g1, g2)

SpatiaLite 封装了 GEOS 库中对应的功能,在几何对象上执行逻辑操作。

8.8.2. 几何对象的更多操作

SpatiaLite支持以下SQL函数:

  • geom2 = ConvexHull ( geom1  )

  • geom2 = Buffer ( geom1 , radius )

  • geom2 = Simplify ( geom1 , tolerance )

  • geom2 = SimplifyPreserveTopology ( geom1 , tolerance )

所有这些函数返回表示所请求操作结果的GEOMETRY, 请注意返回的GEOMETRY可能是OpenGis标准定义了ConvexHull()函数和Buffer()函数;

现在我们来看一下空间操作。在一个几何对象上进行操作,然后生成另外一个几何对象。

  • ConvexHull : 最小凸包,给定一组点 X,返回包含这些点 X 的最小凸包。

  • Buffer : “缓冲区操作”,返回任何几何对象一定距离之内的特定区域,返回的是多边形要素。

  • Simplify : “图形概化操作”,返回一个简化的几何图形。这个图形包含了较少的“顶点”,但是仍然尽量保持了令人满意的近似形状。

下面是一些实际的例子, 来更好地理解这些概念:

最小凸包

g2 = ConvexHull(g1)

最小凸包的概念,可以比较容易地想象, 可以通过想象伸展的弹性带包围给定物体来容易地看出凸包; 当释放时,它将呈现所需凸包的形状。

ConvexHull()可以应用于任何类型的几何,因为我们可以假设点。

关于凸包的概念,中译「凸包」或「凸壳」。在多维空间中有一群散佈各处的点,「凸包」是包覆这群点的所有外壳当中,表面积暨容积最小的一个外壳,而最小的外壳一定是凸的。 至于「凸」的定义是:图形内任意两点的连线不会经过图形外部。「凸」并不是指表面呈弧状隆起,事实上凸包是由许多平坦表面组成的。

缓冲区操作

g2 = Buffer(g1, radius)

在GIS中, 缓冲区操作会在点、线,或面的周围,并根据指定的宽度创建一个区域。

这个缓冲区也指在更普遍的几何对象 (也包括复杂对象,如复合点、复合线等)的指定宽度的区域。

这个实例展示了在同一个几何对象上,随着“半径”的增加,而产生的不同的缓冲区。

Buffer()可以应用于任何类型的几何。 请求的半径表示参考几何使用的相同单位。

几何概化操作

g2 = Simplify(g1, tolerance)

Douglas Peuker 算法在这里用来概化几何对象。 在 SpatiaLite 中,封装了 GEOS 库的功能,用来在几何对象上执行这样的操作。

Douglas Peuker算法试图保持线的方向趋势,根据所需的简化量变化,可以使用一个公差因子。

Simplify() 函数和 SimplifyPreserveTopology() 函数并没有在 PostGIS 中定义, 但是在很多空间数据库中,大都已经实现了。

Simplify() 只能应用在 线性 几何要素上(这也包括多边形的边界), 而不能在点或复合点上应用。

公差因子确定了所需的简化量。[即: 最大的容差值比最小容差值产生最强的简化]

公差必须以引用几何所使用的相同单位表示 SimplifyPreserveTopology()表示Simplify()的专用版本, 需要特别小心才能完全保留原始拓扑。

例如:最适合内孔在内的多边形。

8.8.3. 评价MBR关系

在开始之前,先运行如下程序:

>>> import os,shutil,stat
>>> import sqlite3 as sqlite
>>> tmp_db = '/tmp/xx_new_db.sqlite'
>>> if os.path.exists(tmp_db):
>>>     os.remove(tmp_db)
>>> shutil.copy("spalite.db", tmp_db)
>>> os.chmod(tmp_db, stat.S_IRUSR + stat.S_IWUSR)
>>>
>>> conn = sqlite.connect(tmp_db)
>>> conn.enable_load_extension(True)
>>> conn.execute('SELECT load_extension("mod_spatialite.so.7")')
>>> cursor = conn.cursor()

每一个几何要素都有其最小边界矩形(MBR)。 需要了解 MBR 的概念和 Envelope() 函数。

下面是一组 MBR 中可能会存在的空间关系:

通过比较两个几何要素的 MBR 关系来计算这两个元素的实际空间关系,这是一种不完全的估计的方式。 另一方面,MBR 的比较是一种粗略但又快速的比较空间关系的方式。 比较任意两个复杂多边形的实际空间关系是非常耗时的任务。MBR 是一种非常快速的方法。 通过比较 MBR可以最大程度地限制比较的范围。 MBR是不相交的两个要素,则这两个要素必然不会相交,这个结果设计到一些集合论的知识, 但是理解起来还是很简单的。一旦已经对级和要素进行了大量的不完全比较, 你可以对可能相互的要素进行进一步的检查。

如果两个 MBR 不相关,则相关的几何要素也必须不相交。 但是如果两个 MBR 重合或相交,则不能肯定这两个要素也重合或相交, 在这种情况下,需要进行进一步的空间关系计算。

在GIS应用中一个非常普通的问题就是如何快速有效地获取属于某个特定图框的所有实体,如画一个多边形地图。 针对上面的问题,需要获取图框内的所有实体,而放弃图框外的所有实现。

>>> recs = cursor.execute("SELECT count (*) FROM prov_capital;")
>>> [print(rec) for rec in cursor]
(34,)
[None]

这一个是NewTowns的“未过滤”数据集记住:对于这个表,GEOMETRY是一个POINT类

>>> sql = '''SELECT count(*) FROM prov_capital WHERE MBRContains(
>>>     GeomFromText('POLYGON((554000 4692000, 770000 4692000, 770000 4925000,
>>>     554000 4925000, 554000 4692000))'), geometry)'''
>>> cursor.execute(sql)
<sqlite3.Cursor at 0x7f80f8c676c0>
>>> [print(rec) for rec in cursor]
(0,)
[None]

SpatiaLite 的 MBRContains() 函数允许限制检索的范围。

GeomFromText('POLYGON(...)') 从文本中构建对应现在图框的 MBR, 如屏幕、栅格数据、纸张或打印机等。

下面看一下图中的示意,然后再进行说明。

使用 BuildMBR 函数

GeomFromText() 在构建几何对象时非常有用的, 但在构建 MBR 时,显得有些复杂,因为重复的坐标输入了两遍。 在 SpatiaLite 中,实现了一种简化的方法,你可以使用 BuildMBR(X1 , Y1 , X2 , Y2) 函数在定义 MBR 。

>>> sql = '''SELECT count(*) FROM prov_capital WHERE MBRContains( BuildMBR(554000, 4692000, 770000, 4925000), geometry);'''
>>> cursor.execute(sql)
<sqlite3.Cursor at 0x7f80f8c676c0>
>>> [print(rec) for rec in cursor]
(0,)
[None]
  • [X1, Y1] 坐标值声明了第一个点;

  • [X2 ,Y2] 坐标值声明了第二个点;

  • 生成的 MBR ,是以这两个点为对角线的矩形。

BuildCircleMBR 函数

还有一个与 BuildMBR 类似的函数,即 BuildCircleMBR(X , Y, radius) 函数。

  • [X, Y] 坐标值声明了一个点;

  • 以这个点为圆心,使用给定半径的圆来构建 MBR 。注意,这是矩形,不要以为会是圆形。

下面是实例:

>>> sql = '''SELECT count (*) FROM prov_capital WHERE MBRContains(BuildMBR (654000,4692000, 770000, 4924000), geometry);'''
>>> cursor.execute(sql)
>>> [print(rec) for rec in cursor]
(0,)
[None]

为了获取一些“缩放”的序列,你可以轻松地使用不同的框作为参数来调用 MBRContains() 函数。 根据需要,每次可以设置宽一点或窄一点的 MBR。 对于自动化数据处理,这样做是很常见的。

要进行空间上的滑动窗口遍历,则需要修改 MBR 的起始点的值。 这种空间上的遍历也是很多空间分析需要用到的,对于一些专业的数据分析, 更需要掌握这个函数的思想与使用。

MBRWithin()函数

>>> sql = '''SELECT count (*) FROM prov_capital WHERE MBRWithin( geometry , BuildMBR (754000, 4692000, 770000, 4924000))'''
>>> res = cursor.execute(sql)
>>> [print(rec) for rec in cursor]
(0,)
[None]

MBRWithin() 函数的作用与 MBRContains() 函数是一样, 但是 MBR 的顺序是相反的。

你可以根据需要,使用更符合上下文语义的函数。其结果都是一样的。

MBRIntersects()函数

>>> sql = '''SELECT count (*) FROM HighWays WHERE MBRIntersects(BuildMBR(754000, 4692000, 770000, 4924000), geometry)'''
>>> cursor.execute(sql)
>>> [print(rec) for rec in cursor]
---------------------------------------------------------------------------

OperationalError                          Traceback (most recent call last)

Cell In [9], line 2
      1 sql = '''SELECT count (*) FROM HighWays WHERE MBRIntersects(BuildMBR(754000, 4692000, 770000, 4924000), geometry)'''
----> 2 cursor.execute(sql)
      3 [print(rec) for rec in cursor]


OperationalError: no such table: HighWays

当访问LINESTRING或POLYGONS或其MULTI集合时, HighWay表具有类型为LINESTRING / MULTILINESTRING的GEOMETRY类, SpatiaLite 的MBRIntersects()函数允许您选择给定帧中的实体。 对于MBRIntersects()函数,引用几何的顺序没有任何影响。