我最喜欢的p注浆路由层之一是水层。我感兴趣的是水从哪里来，流向哪里，径流发生在哪里，以及城市发展如何与这种强大的自然力量相互作用。然而，在使用PostGIS和p注浆:Polygons进行路由时，OpenStreetMap水层提出了一个挑战。
为什么多边形是一个挑战?使用p注浆的路由网络是由线(边)构建的。现在，申明一个显而易见的事实:多边形不是直线。
现实世界的水路网络是由直线和多边形组成的。河流、小溪和排水路线主要(但不是全部!)使用线条绘制。这些管线进出池塘、湖泊和水库。下面的动画展示了水多边形对水路网络的影响。当一些非常重要的路径被排除在外时，它们就会消失。
为了使整个水网络可路由，我们需要创建一个组合线路层。组合线层将包括:
-
初始osm。water_line输入
-
来自osm.water_polygon的中轴线
-
连接初始输入和中间轴的线
这篇文章的篇幅很大，因为其中包含了查询示例和截图。
扩展和数据
本文使用了st_approximate atemedialaxis()函数，该函数需要postgis_sfcgal扩展。创建扩展，以及postgis和p注浆，以便它们在需要时可用。
CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;
CREATE EXTENSION postgis_sfcgal;
本文的数据源是使用PgOSM Flex的Docker程序加载的科罗拉多州的OpenStreetMap数据。使用默认层集加载Colorado通常需要不到10分钟的时间来获取.sql，我将其加载到我选择的Postgres实例中。确保使用PgOSM Flex 0.6.0或更高版本。旧版本的PgOSM Flex不包括Standley Lake等水资源关系。
示例湖多边形和中轴
在上面的动画中消失的最大的水体是斯坦德利湖，也被称为3023193关系。可以使用PostGIS函数st_approatemedialaxis()为多个多边形生成合理的移动路径。另一个选项是st_直骨架();这个功能在本文中没有讨论。下面的查询和图像显示了如何计算Standley Lake的中轴。
SELECT osm_id, geom, ST_ApproximateMedialAxis(geom)
FROM osm.water_polygon
WHERE osm_id = -3023193
;
定义小的测试区域
路由的步骤在大型数据集上可能花费很长时间。为了这篇文章的目的，我将科罗拉多的全部数据限制在史坦德利湖大约5公里的重点区域内。我创建了一个名为routing_demo的模式来保存我们的数据，并创建了一个表(routing_demo.focus)来存储焦点的几何形状。在本文的多个查询中使用了焦点区域，以限制使用的数据量。
CREATE SCHEMA routing_demo;
COMMENT ON SCHEMA routing_demo IS 'Contains examples for routing through polygons using ST_ApproximateMedialAxis()';
CREATE TABLE routing_demo.focus AS
SELECT osm_id, geom, ST_Buffer(geom, 5000) AS geom_buffer
FROM osm.water_polygon
WHERE osm_id = -3023193
;
CREATE INDEX ON routing_demo.focus
USING GIST (geom);
COMMENT ON TABLE routing_demo.focus IS 'Focus area for routing';
新创建的routing_demo。焦点表包含建立焦点区域所需的数据。湖泊多边形和计算的缓冲区如下图所示。
SELECT * FROM routing_demo.focus;
我要做的下一步是检查osm。water_polygon数据连接到routing_demo。焦点表。下面的查询显示，我希望从以下步骤中排除两种类型的水多边形：水坝和船坞。虽然水坝对水流有重大影响，但它们不是我想要实现的路线分析的一部分。
SELECT p.osm_type, p.osm_subtype, COUNT(*)
FROM osm.water_polygon p
INNER JOIN routing_demo.focus f ON p.geom && f.geom_buffer
GROUP BY p.osm_type, p.osm_subtype
ORDER BY p.osm_subtype
;
osm_type|osm_subtype|count|
--------+-----------+-----+
waterway|boatyard | 1|
waterway|dam | 2|
natural |water | 172|
natural |wetland | 38|
为多边形创建中间轴
定义了焦点区域后，下一步是计算和持久化这些多边形的中间轴线。下面的查询创建了routing_demo。Polygon_medial_axis表，其中包含构成每个多边形的中间轴的单独线条。st_similatemedialaxis()函数返回MULTILINESTRINGs，因此还使用ST_Dump()来提取每个多边形的单个线组件。
CREATE TABLE routing_demo.polygon_medial_axis AS
SELECT p.osm_id,
(ST_Dump(ST_ApproximateMedialAxis(p.geom))).path[1] AS ma_id,
(ST_Dump(ST_ApproximateMedialAxis(p.geom))).geom AS geom
FROM osm.water_polygon p
INNER JOIN routing_demo.focus f ON p.geom && f.geom_buffer
WHERE p.osm_subtype NOT IN ('dam', 'boatyard')
;
CREATE INDEX ON routing_demo.polygon_medial_axis
USING GIST (geom);
COMMENT ON TABLE routing_demo.polygon_medial_axis IS 'Contains calculated medial axis lines created from the osm.water_polygon table.';
210个多边形的输入创建了4415条中轴线，平均每个多边形21条线。
SELECT COUNT(*) AS rows, COUNT(DISTINCT osm_id) AS input_polygons,
COUNT(*) / COUNT(DISTINCT osm_id) AS avg_medial_axis_per_input
FROM routing_demo.polygon_medial_axis
;
rows|input_polygons|avg_medial_axis_per_input|
----+--------------+-------------------------+
4415| 210| 21|
下面的查询和屏幕截图提供了保存到routing_demo中的中轴数据。polygon_medial_axis表格。
SELECT * FROM routing_demo.polygon_medial_axis;
确定与水多边形相交的水线
以上步骤已经提供了通过我们的水多边形的基本路线。然而，内侧轴线与水道进出多边形的位置不匹配。为了解决这一挑战，我们需要确定osm的水线在哪里。Water_line与水多边形图层相交。为了实现这一点，我重用了我的帖子中关于在OpenStreetMap数据中寻找缺失的交叉标记的技巧。
在下面的查询中，相交CTE查找与多边形相交的线。第一个查询还连接到焦点表，以将范围限制到我们感兴趣的区域。node_id查询使用ST_PointN()和generate_series()从相交于多边形多个点的直线中提取节点id。组合CTE查询将与一个相交节点(来自交集)的相交线和与多个相交节点(来自node_id)的相交线组合在一起。查询后的截图显示斯坦德利湖周围的相交点为蓝点。
CREATE TABLE routing_demo.line_polygon_intersection AS
WITH intersections AS (
SELECT p.osm_id AS water_polygon_osm_id,
l.osm_id,
l.osm_type, l.osm_subtype,
ST_Intersection(p.geom, l.geom) AS geom
FROM osm.water_polygon p
INNER JOIN routing_demo.focus f ON p.geom && f.geom_buffer
INNER JOIN osm.water_line l
ON ST_Intersects(p.geom, l.geom)
AND l.osm_subtype NOT IN ('dam', 'boatyard')
), node_ids AS (
SELECT osm_id, generate_series(1, ST_Npoints(geom)) AS node_id
FROM intersections
WHERE ST_NPoints(geom) > 1
), combined AS (
SELECT DISTINCT i.water_polygon_osm_id, i.osm_id, i.osm_type, i.osm_subtype,
CASE WHEN n.osm_id IS NULL THEN i.geom
ELSE ST_PointN(i.geom, n.node_id)
END AS geom
FROM intersections i
LEFT JOIN node_ids n ON i.osm_id = n.osm_id
)
SELECT *
FROM combined
状态检查
为了查看到目前为止所取得的进展，我从routing_demo中提取了数据。polygon_medial_axis routing_demo。line_polygon_intersection表导入QGIS。下面的截图显示了一条水线与斯丹利湖多边形相交的蓝色点。中间轴有一个分支，它靠近但不接触多边形边界的交点。要创建完整的路由，我们需要将各个部分连接起来。
连接片
下面的查询使用st_closepoint()和ST_MakeLine()将水道与多边形内部的中间轴线相交的点连接起来。INNER JOIN中的子查询使用ST_Collect()将中间轴的各个线折叠为每个多边形的单个几何图形。跳过此步骤将导致计算到每一行的最近点的行。
CREATE TABLE routing_demo.water_line_to_polygon AS
SELECT i.osm_id, i.water_polygon_osm_id,
ST_MakeLine(i.geom, ST_ClosestPoint(ma.geom, i.geom)) AS geom
FROM routing_demo.line_polygon_intersection i
INNER JOIN (
SELECT osm_id, ST_Collect(geom) AS geom
FROM routing_demo.polygon_medial_axis
GROUP BY osm_id
) ma
ON i.water_polygon_osm_id = ma.osm_id
;
把它放在一起上述步骤创建了我们需要缝合在一起的组件。重新查看文章开头的项目列表，下表显示了数据元素当前的位置。
下面的查询将这三（3）个数据集合并为一个行表。这一步的目标是为pgRouting步骤做好一切准备。
DROP TABLE IF EXISTS routing_demo.waterways_to_route;
CREATE TABLE routing_demo.waterways_to_route AS
WITH ww_selected AS (
SELECT l.osm_id, l.osm_type, l.osm_subtype,
NULL::BIGINT AS water_polygon_osm_id,
l.layer,
-- This will simplify some/most...
ST_LineMerge(l.geom) AS the_geom
FROM routing_demo.focus f
INNER JOIN osm.water_line l
ON f.geom_buffer && l.geom
AND l.osm_subtype NOT IN ('dam', 'boatyard')
), ww_extra_cleanup AS (
-- Split the complex lines that made it through ST_LineMerge
SELECT osm_id, osm_type, osm_subtype, water_polygon_osm_id,
layer,
(ST_Dump(the_geom)).geom AS the_geom
FROM ww_selected
WHERE ST_GeometryType(the_geom) = 'ST_MultiLineString'
), combined AS (
SELECT osm_id, osm_type, osm_subtype, water_polygon_osm_id, layer,
the_geom
FROM ww_selected
WHERE ST_GeometryType(the_geom) != 'ST_MultiLineString'
UNION
SELECT osm_id, osm_type, osm_subtype, water_polygon_osm_id, layer,
the_geom
FROM ww_extra_cleanup
WHERE ST_GeometryType(the_geom) != 'ST_MultiLineString'
UNION
SELECT NULL AS osm_id, 'line-to-polygon' AS osm_type, '' AS osm_subtype,
water_polygon_osm_id, 0 AS layer,
geom AS the_geom
FROM routing_demo.water_line_to_polygon
UNION
SELECT NULL AS osm_id, 'medial-axis' AS osm_type, '' AS osm_subtype,
osm_id AS water_polygon_osm_id, 0 AS layer,
geom AS the_geom
FROM routing_demo.polygon_medial_axis
)
SELECT ROW_NUMBER() OVER () AS id, *
FROM combined
WHERE the_geom IS NOT NULL
;
路由网络
以下步骤创建路由所需的结构和数据。在这篇文章中没有足够的篇幅来深入讨论这些步骤。
我的《精通PostGIS和OpenStreetMap》一书的第14、15、16和18章都是关于p注浆的，并详细解释了这个过程。
SELECT pgr_nodeNetwork('routing_demo.waterways_to_route', .1);
-- Establish route costs
ALTER TABLE routing_demo.waterways_to_route_noded
ADD cost NUMERIC
GENERATED ALWAYS AS (
ST_Length(the_geom)
)
STORED
;
ALTER TABLE routing_demo.waterways_to_route_noded
ADD reverse_cost NUMERIC NOT NULL DEFAULT -1
;
UPDATE routing_demo.waterways_to_route_noded n
SET reverse_cost = ST_Length(n.the_geom)
FROM routing_demo.waterways_to_route r
WHERE n.old_id = r.id
AND r.osm_type IN ('line-to-polygon', 'medial-axis')
;
-- Topology and Analyze
SELECT pgr_createTopology('routing_demo.waterways_to_route_noded', 0.1);
SELECT pgr_analyzeGraph('routing_demo.waterways_to_route_noded', 0.1);
创建一条穿过湖泊的路线
唷，我们成功了。是时候创建一条路线了!下面的查询从所选的开始和结束节点创建一条从西到东穿过Standley Lake的路线。
SELECT *
FROM pgr_dijkstra('
SELECT rn.id, rn.source,
rn.target, rn.cost,
rn.reverse_cost
FROM routing_demo.waterways_to_route_noded rn
WHERE source IS NOT NULL
',
235, -- start node
5930, -- end node
True -- Enforce one-way restrictions
) d
INNER JOIN routing_demo.waterways_to_route_noded x
ON d.edge = x.id
;
总结
这篇文章深入介绍了如何使用OpenStreetMap数据创建一个能够穿越多边形的水路路由网络。虽然st_approatemedialaxis()是涉及的一个关键函数，但这篇文章展示了十几个函数的使用!该技术为使用该地层进行更完整、更复杂的水文分析奠定了基础。
在水路层之外，这种方法还可以使用OpenStreetMap中的其他层数据来改进路由。行人路线是另一个需要包含多边形的地方。要关闭这篇文章，再一个查询和一个截图。坚持了下来，值得称赞!
SELECT osm_id, osm_type, route_foot, geom
FROM osm.road_polygon
WHERE osm_id = 323692586
;
原文标题：Routing with Lines through Polygons
原文作者： Ryan Lambert
原文地址：https://blog.rustprooflabs.com/2022/10/pgrouting-lines-through-polygons
