Thursday, 7 March 2019

Grouping connected linestrings in PostGIS?


I have a table of streets I've selected based on a set of attributes (let's say it's speed_limit < 25). There are groups of streets that are locally contiguous; I'd like to group these sets of connected linestrings into GeometryCollections. In the image below, there would be two GeometryCollections: one with the red lines and one with the blue lines.


enter image description here


I tried running a couple of "dissolve, deaggregate" queries along the lines of:


SELECT (ST_Dump(st_union)).geom
FROM
(SELECT ST_Union(geom) FROM roads) sq

With everything I've tried, I either end up with a single feature (ST_Union) or my original geometry (ST_Dump of ST_Union).


Maybe it's possible to do this with some kind of WITH RECURSIVE magic?




Answer



So, by example. Here's a simple table with two connected groups of edges:


drop table lines;
create table lines ( id integer primary key, geom geometry(linestring) );
insert into lines (id, geom) values ( 1, 'LINESTRING(0 0, 0 1)');
insert into lines (id, geom) values ( 2, 'LINESTRING(0 1, 1 1)');
insert into lines (id, geom) values ( 3, 'LINESTRING(1 1, 1 2)');
insert into lines (id, geom) values ( 4, 'LINESTRING(1 2, 2 2)');
insert into lines (id, geom) values ( 11, 'LINESTRING(10 10, 10 11)');
insert into lines (id, geom) values ( 12, 'LINESTRING(10 11, 11 11)');

insert into lines (id, geom) values ( 13, 'LINESTRING(11 11, 11 12)');
insert into lines (id, geom) values ( 14, 'LINESTRING(11 12, 12 12)');
create index lines_gix on lines using gist(geom);

Now, here's a recursive function that, given the id of an edge, accumulates all the edges that touch:


CREATE OR REPLACE FUNCTION find_connected(integer) returns integer[] AS
$$
WITH RECURSIVE lines_r AS (
SELECT ARRAY[id] AS idlist, geom, id
FROM lines

WHERE id = $1
UNION ALL
SELECT array_append(lines_r.idlist, lines.id) AS idlist,
lines.geom AS geom,
lines.id AS id
FROM lines, lines_r
WHERE ST_Touches(lines.geom, lines_r.geom)
AND NOT lines_r.idlist @> ARRAY[lines.id]
)
SELECT

array_agg(id) AS idlist
FROM lines_r
$$
LANGUAGE 'sql';

That just leaves us needing to find, after each group is accumulated, the id of an edge that is not already part of a group. Which, tragically, requires a second recursive query.


WITH RECURSIVE groups_r AS (
(SELECT find_connected(id) AS idlist,
find_connected(id) AS grouplist,
id FROM lines WHERE id = 1)

UNION ALL
(SELECT array_cat(groups_r.idlist,find_connected(lines.id)) AS idlist,
find_connected(lines.id) AS grouplist,
lines.id
FROM lines, groups_r
WHERE NOT idlist @> ARRAY[lines.id]
LIMIT 1)
)
SELECT id, grouplist
FROM groups_r;


Which taken together return a nice set with the seed id and each group it accumulated. I leave it as an exercise to the reader to turn the arrays of id back into a query to create geometry for mapping.


 id |   grouplist   
----+---------------
1 | {1,2,3,4}
11 | {11,12,13,14}
(2 rows)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...