Wednesday, July 21, 2010

Network Walking in PostGIS

One of the new features in PostgreSQL 8.4 was the “WITH RECURSIVE” clause available for queries. It allows you to define a subquery based on a recursive term — fancy language for a function that calls itself. One of the favorite uses of recursion is walking a network. Geospatial applications use networks all the time: electrical grids, stream systems, and storm sewers are all directed networks (they have unidirectional flow).

Here’s an example of network walking using a simple collection of segments. As is common in many GIS applications, the segment are implicitly connected — their end points are coincident with the start points of other segments.

CREATE TABLE network ( segment geometry, id integer primary key );

INSERT INTO network VALUES ('LINESTRING(1 1, 0 0)', 1);
INSERT INTO network VALUES ('LINESTRING(2 1, 1 1)', 2);
INSERT INTO network VALUES ('LINESTRING(1 2, 1 1)', 3);
INSERT INTO network VALUES ('LINESTRING(3 1, 2 1)', 4);
INSERT INTO network VALUES ('LINESTRING(3 2, 2 1)', 5);
INSERT INTO network VALUES ('LINESTRING(2 3, 1 2)', 6);
INSERT INTO network VALUES ('LINESTRING(1 3, 1 2)', 7);
INSERT INTO network VALUES ('LINESTRING(4 2, 3 2)', 8);
INSERT INTO network VALUES ('LINESTRING(3 4, 2 3)', 9);
INSERT INTO network VALUES ('LINESTRING(2 4, 2 3)', 10);
INSERT INTO network VALUES ('LINESTRING(1 4, 1 3)', 11);
INSERT INTO network VALUES ('LINESTRING(4 3, 4 2)', 12);
INSERT INTO network VALUES ('LINESTRING(4 4, 3 4)', 13);

CREATE INDEX network_gix ON network USING GIST (segment);

Visually, the network looks like this:

Network

To walk our network, use a WITH clause that starts with one segment, then repeatedly adds the next downstream segment to the collection. In our case, the “next downstream segment” is defined as a segment whose start point is close to the end point of the current segment. We’ll walk down from segment 6.

WITH RECURSIVE walk_network(id, segment) AS (
SELECT id, segment FROM network WHERE id = 6
UNION ALL
SELECT n.id, n.segment
FROM network n, walk_network w
WHERE ST_DWithin(ST_EndPoint(w.segment),ST_StartPoint(n.segment),0.01)
)
SELECT id
FROM walk_network

Which returns:

 id
----
6
3
1
(3 rows)

From 6 to 3 to 1, correct! Once we have our walker producing the results we want, we can wrap more PostGIS and PostgreSQL functions around the walker to produce a finished product. Here’s a function that takes in an edge identifier and outputs a linestring based on the downstream path.

CREATE OR REPLACE FUNCTION downstream(integer)
RETURNS geometry
LANGUAGE sql
AS '
WITH RECURSIVE walk_network(id, segment) AS (
SELECT id, segment FROM network WHERE id = $1
UNION ALL
SELECT n.id, n.segment
FROM network n, walk_network w
WHERE ST_DWithin(ST_EndPoint(w.segment),ST_StartPoint(n.segment),0.01)
)
SELECT ST_MakeLine(ST_EndPoint(segment))
FROM walk_network
' IMMUTABLE;

And here’s the function in action, generating the downstream path from segment 12.

=# SELECT ST_AsText(Downstream(12));

st_astext
---------------------------------
LINESTRING(4 2,3 2,2 1,1 1,0 0)
(1 row)

Check the generated path against our network picture — looks good!

Path 12

9 comments:

Ben Reilly said...

Very cool. I think I asked you how this might be done when you were up in Washington State (WAURISA) a couple of weeks back. Happy to see a post on it.

I haven't gotten much of a chance to work on the problem, but the direction I've been heading to toward a marriage of GDAL/OGR and networkx in Python; something like this.

That might be going about it the hard way, but it's as much for brushing up on the various open source libraries and Python as it is for the end product.

selectoid said...

This is a superb example of the possibilities WITH RECURSIVE gives you. I guess the performance of this approach should be rather good even on bigger networks, right?

nicklas.aven said...

Very nice!
But I guess you will not get the startpoint of edge 12 as the illustration shows in the second example?

Paul Ramsey said...

@nicklas, you're right my graphic doesn't match my result exactly :(

@selectoid, on performance I have no numbers, but am interested in trying this on something larger

snorris said...

Very nice.

I once tried Oracle's recursive queries to select downstream watershed codes in the BC freshwater atlas stream network. I think it worked, but the performance wasn't nearly as good as simply parsing the watershed code in Python. (and trying to code my queries in pl/sql was unpleasant).

I'll have to try this option out - it looks like it will work with side channels and distributatires, where watershed code logic can break down.

William said...

Given the recursive nature of the example in your post I presume that it would not be possible to use WITH RECURSIVE to traverse cyclic graphs as there's no way of determining which edges have already been visited?

I ask as I wish to associate data with nodes and edges that represent a water distribution network then traverse my networks to find the number of entities of one type that lie upstream of entities of another type to assess causality relationships. Water networks are typically modeled as cyclic directed graphs.

Paul Ramsey said...

@William, the PostgreSQL 'with' documentation section describes some methods of dealing with cycles. http://www.postgresql.org/docs/8.4/static/queries-with.html

Marcello Benigno said...

Paul, I'm trying adapted your example to find the longest river inside one basin, so, the result it's exactly the inverse that I need (like your solution), I made some unsuccessfully tests using ST_Reverse() in my geometry, So... Can you help me with this solution? here is my fail logic:

-- Create a View
CREATE OR REPLACE VIEW river_inv AS
SELECT river.gid, st_reverse(river.the_geom) AS the_geom
FROM river;

--Using your query
WITH RECURSIVE walk_network(gid, the_geom) AS (
SELECT gid, the_geom
FROM river_inv
WHERE gid = 212
UNION ALL
SELECT n.gid, n.the_geom
FROM river n, walk_network w
WHERE ST_DWithin(ST_EndPoint(w.the_geom),ST_StartPoint(n.the_geom),0.01)
)
SELECT *
FROM walk_network;

Thanks for shared your knowledge.

pcreso said...

How might this approach be used to identify all the upstream reaches?

About Me

My Photo
Victoria, British Columbia, Canada

Followers

Blog Archive

Labels