31. 拓扑¶
PostGIS 通过名为 postgis_topology 的扩展支持 SQL/MM SQL-MM 3 Topo-Geo 和 Topo-Net 3 规范。您可以在 手册:PostGIS 拓扑 中了解此扩展提供的所有函数和类型。postgis_topology 扩展包含另一种类型的核心空间类型,称为 拓扑几何。除了 拓扑几何 空间类型之外,您还会发现用于构建 拓扑 和填充拓扑的函数。
在开始使用拓扑之前,您必须安装 postgis_topology 扩展,如下所示
CREATE EXTENSION postgis_topology;
安装扩展后,您将在数据库中看到一个名为 topology 的新模式。topology 模式将对数据库中的所有拓扑进行编目。
topology 模式包含两个表和所有用于拓扑的辅助函数。
topology - 列出数据库中的所有拓扑,以及它们存储在哪个模式中
layer - 列出数据库中所有包含拓扑几何的表列
layer 表与我们之前了解的 raster_columns、geometry_columns 和 geography_columns 目录非常相似,但专门针对拓扑几何。
31.1. 创建拓扑¶
拓扑和拓扑几何到底是什么,它们之间有什么关系?在解释之前,让我们先创建一个拓扑来容纳我们的纽约市拓扑完美数据,使用 CreateTopology 函数并将容差设置为 0.5 米。请注意,0.5 是以米为单位,因为我们的空间参考系统是纽约州平面米。
SELECT topology.CreateTopology('nyc_topo', 26918, 0.5);
哪些输出
1
这是它分配给新拓扑的 ID。运行上述命令后,您将在数据库中看到一个名为 nyc_topo 的新模式。您可以将拓扑命名为任何您想要的名字。我的惯例是在末尾添加 _topo 以将其与数据库中的其他模式区分开来。
如果您探索 topology.topology 表格,
SELECT * FROM topology.topology;
您将看到
id | name | srid | precision | hasz
----+----------+-------+-----------+------
1 | nyc_topo | 26918 | 0 | f
(1 row)
31.2. 拓扑和拓扑几何的存储¶
拓扑在 PostgreSQL 数据库中实现为一个模式。如果您探索 nyc_topo 模式,您将看到这些表格和视图
- edge - 这是针对 edge_data 构建的视图,主要用于 SQL/MM 兼容性。
它包含 edge_data 表格中的一组列。
edge_data - 包含构成拓扑的所有线字符串。
- face - 包含所有可以从 edge_data 形成的封闭表面的列表。
它不包含实际的几何图形,而是只包含几何图形的边界框。
node - 包含所有边的起点和终点,以及未连接到任何东西的点(孤立节点)。
relation - 这定义了拓扑中的哪些元素构成拓扑几何。
那么什么是拓扑几何?拓扑几何是通过拓扑中的边、面、节点和其他拓扑几何形成的几何图形的表示。
拓扑几何位于哪里?它位于其他地方,通过 relation 表格引用拓扑的元素。虽然我们可以将拓扑几何放在我们的 nyc_topo 模式中,但一般惯例是定义其他模式中的其他表格,这些表格具有拓扑几何,并且还具有您可能感兴趣的任何其他类型的数据。
31.3. 为什么要使用拓扑几何?¶
使用拓扑几何可以使您的数据井井有条并保持连接。拓扑几何对于地籍工作非常有用,在这些工作中,您希望确保两块土地不会重叠,即使您更改了其中一块的边界,或者您希望确保道路在您更改形成它们的几何图形时保持连接。几何图形独立存在,您可以复制它们,变形它们。几何图形无忧无虑,不关心与它们共享空间的其他几何图形。相反,拓扑几何遵循其拓扑规则;除非存在定义它们的边、节点、面或其他拓扑几何,否则它们无法存在。拓扑几何属于一个且仅属于一个拓扑。拓扑几何是几何图形的关系模型,因此,当每个组件(边/面/节点)移动、添加等时,它们不会改变一个拓扑几何形状,而是改变所有具有共同组件的拓扑几何形状。
我们有一个名为 nyc_topo 的拓扑结构,它没有任何数据。让我们用我们的纽约市数据填充它。拓扑结构的边、面和节点可以通过两种主要方式创建。
可以使用拓扑结构基元函数直接创建边、面和节点。
可以通过创建拓扑几何来形成边、面和节点。当从几何图形创建拓扑几何图形时,如果缺少与坐标匹配的边、面或节点,则在创建过程中会创建新的边、面和节点。
31.4. 定义拓扑几何列并创建拓扑几何¶
填充拓扑结构最常见的方法是创建拓扑几何。让我们从创建一个用于保存街区的表开始,然后使用 AddTopoGeometryColumn 函数添加一个拓扑几何列。
CREATE TABLE nyc_neighborhoods_t(boroname varchar(43), name varchar(67),
CONSTRAINT pk_nyc_neighborhoods_t PRIMARY KEY(boroname,name) );
SELECT topology.AddTopoGeometryColumn('nyc_topo', 'public', 'nyc_neighborhoods_t',
'topo', 'POLYGON') As layer_id;
以上操作的输出为
layer_id
--------
1
现在我们准备填充我们的表。最好在添加之前确保您的几何图形有效,否则您会遇到错误,例如 SQLMM 几何图形不简单。
因此,让我们从添加有效的几何图形开始。这里使用的 1 指的是从上一个查询生成的 layer_id。如果您不知道 layer_id,您可以使用 FindLayer 函数查找它,我们将在后面的示例中使用它。
对于这些示例,您将使用 toTopoGeom 函数将几何图形转换为其拓扑几何等效项。toTopoGeom 函数为您处理了许多簿记工作。
toTopoGeom 函数会检查传入的几何图形,并在您的拓扑结构中根据需要注入节点、边和面,以形成几何图形的形状。然后,它会将关系添加到 relation 表中,以定义此新的拓扑几何与这些新的和现有的拓扑元素之间的关系。在某些情况下,几何图形的部分可能存在,或者需要拆分现有部分以形成新的几何图形。
INSERT INTO nyc_neighborhoods_t(boroname,name, topo)
SELECT boroname, name, topology.toTopoGeom(geom, 'nyc_topo', 1)
FROM nyc_neighborhoods
WHERE ST_ISvalid(geom);
上述步骤应该需要 3-4 秒。现在让我们添加无效的几何图形
INSERT INTO nyc_neighborhoods_t(boroname,name, topo)
SELECT boroname, name, topology.toTopoGeom(
ST_UnaryUnion(
ST_CollectionExtract(
ST_MakeValid(geom), 3)
), 'nyc_topo', 1)
FROM nyc_neighborhoods
WHERE NOT ST_ISvalid(geom);
上述步骤应该需要大约 300-400 毫秒。
现在我们的拓扑结构中有数据了。快速检查将显示 nyc_topo.edge、nyc_topo.node 和 nyc_topo.face 都有数据
SELECT 'edge' AS name, count(*)
FROM nyc_topo.edge
UNION ALL
SELECT 'node' AS name, count(*)
FROM nyc_topo.node
UNION ALL
SELECT 'face' AS name, count(*)
FROM nyc_topo.face;
输出
name | count
------+-------
edge | 580
node | 396
face | 218
(3 rows)
现在我们可以通过在 nyc_boros_t 表中定义一个名为 topo 的列来声明性地表达行政区是由多个街区组成的。该列的类型为 POLYGON,并且是来自 nyc_neighborhoods_t.topo 列的其他拓扑几何的集合。
CREATE TABLE nyc_boros_t(boroname varchar(43),
CONSTRAINT pk_nyc_boros_t PRIMARY KEY(boroname) );
SELECT topology.AddTopoGeometryColumn('nyc_topo', 'public', 'nyc_boros_t',
'topo', 'POLYGON',
(topology.FindLayer('public', 'nyc_neighborhoods_t', 'topo')).layer_id
) AS layer_id;
哪些输出
为了填充这个新表,我们将使用 CreateTopoGeom 函数。CreateTopoGeom 函数不是从几何图形开始构建新的拓扑几何,而是从现有的拓扑元素开始,这些元素可能是基本几何图形或其他拓扑几何图形,以定义新的拓扑几何。
INSERT INTO nyc_boros_t(boroname, topo)
SELECT n.boroname,
topology.CreateTopoGeom('nyc_topo',
3, (topology.FindLayer('public', 'nyc_boros_t', 'topo')).layer_id ,
topology.TopoElementArray_Agg( ARRAY[ (n.topo).id, (n.topo).layer_id ]::topoelement ) )
FROM nyc_neighborhoods_t AS n
GROUP BY n.boroname;
这将插入 5 条记录,对应于纽约的五个行政区。
注意
如果您使用的是 PostGIS 3.4 或更高版本,您可以使用新的强制转换将拓扑几何转换为拓扑元素,并将上面示例中的 topology.TopoElementArray_Agg( ARRAY[ (n.topo).id, (n.topo).layer_id ]::topoelement ) ) 替换为更短的 topology.TopoElementArray_Agg( n.topo::topoelement )
要在 pgAdmin 中查看这些内容,您可以将拓扑几何转换为几何图形,如下所示
SELECT boroname, topo::geometry AS geom
FROM nyc_boros_t;
输出将如下所示
如果您认为这太乱了,是的,确实很乱。这是在多次简化和其他几何处理之后发生的,其中每个几何图形都被视为一个独立的单元。您会得到间隙、悬空岛屿以及街区相互侵犯领土的情况。
幸运的是,我们可以使用拓扑来清理这种混乱,并帮助我们维护干净、连接良好的数据。
让我们戴上土地测量员的帽子,并提出一个问题,如果我们将土地划分为区域(行政区或街区),使得每个区域可能与其他区域接壤,但不能共享任何共同区域,那么区域共享共同区域是否有意义?没有意义。而我们的数据指出了某些区域属于多个街区或多个行政区。
首先让我们看看行政区,并寻找共享共同元素的街区
SELECT te, array_agg(DISTINCT b.boroname)
FROM nyc_boros_t AS b, topology.GetTopoGeomelements(topo) AS te
GROUP BY te
HAVING count(DISTINCT b.boroname) > 1;
输出是
te | array_agg
--------+-------------------
{44,3} | {Brooklyn,Queens}
{51,3} | {Brooklyn,Queens}
{76,3} | {Brooklyn,Queens}
{114,3} | {Brooklyn,Queens}
{117,3} | {Brooklyn,Queens}
(5 rows)
这告诉我们皇后区和布鲁克林正在进行边界战争。在这个查询中,我们使用 GetTopoGeomElements 函数来声明性地检查跨行政区共享的组件。
返回的是一组拓扑元素。拓扑元素表示为一个包含两个整数的数组,第一个数字是元素的 ID,第二个数字是元素的层级(或基本类型)。PostGIS GetTopoElements 返回拓扑几何体的基本类型,类型号 1-3 分别对应于(1:节点,2:边,3:面)。所有街区和行政区对应的拓扑元素类型都是 3,对应于拓扑面。我们可以使用 ST_GetFaceGeometry 来获取这些共享面的可视化表示,如下所示
SELECT te, t.geom, ST_Area(t.geom) AS area, array_agg(DISTINCT d.boroname) AS shared_boros
FROM nyc_boros_t AS d, topology.GetTopoGeomelements(topo) AS te
, topology.ST_GetFaceGeometry('nyc_topo',te[1]) AS t(geom)
GROUP BY te, t.geom
HAVING count(DISTINCT d.boroname) > 1
ORDER BY area;
结果将是 5 行,对应于皇后区和布鲁克林之间的边界争议。
如果我们查看我们的街区,我们会看到类似的情况,但有 44 个边界争议。
SELECT te, t.geom, ST_Area(t.geom) AS area, array_agg(DISTINCT d.name) AS shared_d
FROM nyc_neighborhoods_t AS d, topology.GetTopoGeomelements(d.topo) AS te
, topology.ST_GetFaceGeometry('nyc_topo',te[1]) AS t(geom)
GROUP BY te, t.geom
HAVING count(DISTINCT d.name) > 1
ORDER BY area;
由于行政区是街区的聚合,我们可以通过解决街区边界争议来解决行政区问题。
我们可以通过多种方式解决这个问题。我们可以进行实地调查,询问人们他们认为自己站在哪个街区。或者,我们可以简单地将土地碎片分配给面积最小的街区或出价最高的街区。
从拓扑几何体中删除元素使用 TopoGeom_remElement 函数处理。因此,让我们开始吧,从面积最大的街区中删除重复的元素,如下所示
WITH to_remove AS (SELECT te, MAX( ST_Area(d.topo::geometry) ) AS max_area, array_agg(DISTINCT d.name) AS shared_d
FROM nyc_neighborhoods_t AS d, topology.GetTopoGeomelements(d.topo) AS te
, topology.ST_GetFaceGeometry('nyc_topo',te[1]) AS t(geom)
GROUP BY te
HAVING count(DISTINCT d.name) > 1)
UPDATE nyc_neighborhoods_t AS d SET topo = TopoGeom_remElement(topo, te)
FROM to_remove
WHERE d.name = ANY(to_remove.shared_d)
AND ST_Area(d.topo::geometry) = to_remove.max_area;
上述操作的结果是更新了 29 个街区。如果你重新运行街区和行政区的边界争议查询,你会发现不再有边界争议。
我们仍然有由于过度简化造成的街区之间的空隙。可以通过使用 拓扑编辑器系列函数 直接编辑拓扑,以及/或者填充空洞并将它们分配给街区来解决此类问题。