Routing with Lines through Polygons
One of my favorite layers to route with pgRouting is the water layer. I am interested in where water comes from, where it goes, where runoff happens, and how urban development interacts with this powerful force of nature. The OpenStreetMap water layer, however, presents a challenge when routing with PostGIS and pgRouting: Polygons.
Why are polygons a challenge? A routing network using pgRouting is built from lines (edges). Now, to state the obvious: polygons are not lines.
Real world waterway networks are made up of both lines and polygons. Rivers, streams, and drainage routes are predominately (but not exclusively!) mapped using lines. These lines feed into and out of ponds, lakes, and reservoirs. The following animation shows how much impact the water polygons can have on a waterway network... some very important paths simply disappear when they are excluded.
To make the full water network route-able we need to create a combined line layer. The combined line layer will include:
- Medial axis lines from
- Lines to connect initial inputs to medial axis
This post is a big one thanks to the query examples and screenshots involved.
This post was inspired by a two-part post written by Bruce Rindahl. The posts, titled "Flow path through a lake along the medial axis," are currently unavailable due to a blog upgrade issue. The original source URL will hopefully be restored in the future: http://berrie.ddns.net/?p=28.
Extensions and Data
This post uses the
which requires the
postgis_sfcgal extension. Create the extension,
pgrouting so they are available as they
CREATE EXTENSION postgis; CREATE EXTENSION pgrouting; CREATE EXTENSION postgis_sfcgal;
The data source for this post is OpenStreetMap data for Colorado
loaded using PgOSM Flex's
Loading Colorado with the default layerset typically takes
less than 10 minutes
to get the
.sql, which I load into my Postgres instance of choice.
Be sure to use PgOSM Flex 0.6.0 or newer. Older versions of PgOSM Flex did not include water relations such as Standley Lake.
My book Mastering PostGIS and OpenStreetMap provides in-depth explanations covering how to load OpenStreetMap data (chapter 8) and prepare it for routing (chapters 14 - 16 and 18). The free preview section shows the book's full table of contents.
Example lake polygon and medial axis
The largest body of water that disappears in the above animation is
Standley Lake, also know as
The PostGIS function
can be used to generate a reasonable path of travel for many polygons.
Another option is ST_StraightSkeleton(); that function is not discussed in this post.
The following query and image shows how to calculate the medial axis
for Standley Lake.
SELECT osm_id, geom, ST_ApproximateMedialAxis(geom) FROM osm.water_polygon WHERE osm_id = -3023193 ;
Define small test area
Routing has steps that can take a long time on large data sets. For the purpose
of this post I limit the full Colorado data down to a roughly 5 kilometer
area of focus are Standley Lake.
I create a schema named
routing_demo to hold our data and a table
routing_demo.focus) to store the geometry of focus. The focus area
is used in multiple queries throughout this post to limit the amount of data
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';
The newly created
routing_demo.focus table contains the data needed to establish
the focus area. The lake polygon and calculated buffer are shown in the following
SELECT * FROM routing_demo.focus;
The next step I take is to check the
osm.water_polygon data joined to the
routing_demo.focus table. The following query shows that there
are two types of water polygons I will want to exclude from the following steps:
boatyard. While dams have a significant impact to the flow of water,
they are not part of the routing analysis I want to achieve.
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|
Create medial axis for polygons
With the focus area defined, the next step is to calculate and persist the
medial axis lines for these polygons. The following query creates the
routing_demo.polygon_medial_axis table with the individual lines that make
up the medial axis for each polygon. The
ST_Dump() is also used
in order to extract the individual line components for each polygon.
CREATE TABLE routing_demo.polygon_medial_axis AS SELECT p.osm_id, (ST_Dump(ST_ApproximateMedialAxis(p.geom))).path 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.';
The input of 210 polygons created 4,415 medial axis lines, an average of 21 lines per polygon.
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|
The following query and screenshot provide a look at the medial axis data
saved to the
SELECT * FROM routing_demo.polygon_medial_axis;
Identify water lines intersecting water polygons
The above steps have provided basic line routes through our water polygons.
The medial axis lines, however, do not match up to where waterways enter and
exit the polygons. To solve that challenge, we need to identify where
the water lines from
osm.water_line intersect with the water polygon layer.
To accomplish this I'm reusing a trick from my post on
finding missing crossing tags
in OpenStreetMap data.
In the following query, the
intersections CTE finds lines that intersect
with polygons. This first query also joins to the
to limit the scope to our region of interest.
node_ids query uses
to extract node ids from lines that intersect the polygon at multiple points.
combined CTE query combines the intersecting lines
with one intersecting node (from
intersections) and those
with multiple intersecting nodes (from
The screenshot after the query shows the intersecting points
around Standley Lake as blue dots.
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 ;
In case you were curious, the points that are contained within polygons are where the waterway line has been extended into (or through) the polygon itself. The wetland area to the west of Standley Lake is an example of that. The query above could be easily modified to exclude these points if the user desired that change. I have decided to leave them in, they won't hurt anything beyond taking up a bit more space.
To see the progress made so far I pull the data from
tables into QGIS.
The following screenshot shows the blue point where a water line intersects
with the Standley Lake polygon. The medial axis has a branch that comes close
to, but does not touch, the intersection at the boundary of the polygon.
To create complete routes, we need to connect the pieces.
Connect the pieces
The following query uses
ST_MakeLine() to bridge
the point where a waterway intersects a polygon and the medial axis line(s)
inside the polygon. The subquery in the
INNER JOIN uses
to collapse the individual lines of the medial axis into a single geometry
per polygon. Skipping this step would result in lines calculated to the closest
point of each line.
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 ;
Bring it together
The steps above created the components we need to stitch together. Revisiting the list of items from the beginning of the post, the following table shows where the data elements are currently located.
|Data needed||Available in|
|Medial axis lines from
|Lines to connect initial inputs to medial axis||
The following query combines these three (3) data sets into a single table
The goal in this step is to make everything ready for the
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 ;
The following steps create the structures and data needed for routing.
There is not room in this post to go into these steps in depth.
Chapters 14, 15, 16 and 18 of my book Mastering PostGIS and OpenStreetMap are all about
pgRouting and explain this process in detail.
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);
Create route through a lake
Phew, we made it. It's time to create a route! The following query creates a route from a chosen start and end node that goes through Standley Lake from west to east.
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 ;
The start and end nodes that I used in this post (235 and 5930) will almost certainly be different for you. Any changes to the data or exact steps above will result in different IDs for the same point from the data.
This overstuffed post gave an in-depth guide to using OpenStreetMap data
to create a waterway routing network capable of traversing polygons.
ST_ApproximateMedialAxis() was a key function involved,
this post shows the use of more than a dozen functions!
This technique sets the stage for more complete and complex hydrology analysis
using this layer.
Beyond the waterway layer, this approach can be used to improve routing using other layers of data available in OpenStreetMap. Pedestrian routes are another place where polygons need to be included. To close this post, one more query and one more screenshot. Kudos for hanging in there!
SELECT osm_id, osm_type, route_foot, geom FROM osm.road_polygon WHERE osm_id = 323692586 ;
Need help with your PostgreSQL servers or databases? Contact us to start the conversation!
Published October 23, 2022
Last Updated November 20, 2022