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)
You will need:
1) A table with LineString
geometries:
CREATE TABLE lin (
id serial PRIMARY KEY,
geom geometry(LineString, 31370)
);
CREATE INDEX ON lin (id);
CREATE INDEX ON lin USING gist (geom);
2) A table with Point
geometries where you want to split your overlapping lines. They can represent train/metro stations, intersections, bifurcations etc.
(edit: if you do not have a table with intersection points, use the query at the end of this answer to generate them)
CREATE TABLE pt (
id serial PRIMARY KEY,
geom geometry(Point, 31370)
);
CREATE INDEX ON pt (id);
CREATE INDEX ON pt USING gist (geom);
For this example I've drawn some lines with overlapping sections and some points. The goal is to cut our lines in segments between the points, and then to count how many segments exists between pt_from
and pt_to
.
Here is the query:
WITH
a AS
(
-- Calculate the linear position (loc) of points on the line
SELECT
lin.id lin_id,
pt.id pt_id,
ST_LineLocatePoint(lin.geom, pt.geom) loc
FROM
lin
JOIN
pt ON ST_DWithin(lin.geom, pt.geom, 100.0)
),
b AS
(
-- Calculate from/to segments on lines by retrieving the next segment
-- on the line with the LEAD window function, ordered by location (loc)
SELECT
lin_id,
pt_id pt_from,
loc loc_from,
LEAD(pt_id) OVER w_b pt_to,
LEAD(loc) OVER w_b loc_to
FROM
a
WINDOW w_b AS (PARTITION BY lin_id ORDER BY loc)
)
-- Use row_number to count the number of the overlapping lines between pt_from and pt_to
-- It can be used to parametrize the offset to be displayed
SELECT
row_number() OVER() fid, -- generate unique fid
pt_from,
loc_from,
pt_to,
loc_to,
lin_id,
row_number() OVER w_c offset_nr, -- iterate on overlapping lines
ST_LineSubstring(lin.geom, loc_from, loc_to)::geometry(LineString, 31370) geom -- cut lines to make segments
FROM
b
JOIN
lin ON b.lin_id = lin.id
WHERE
pt_to IS NOT NULL
WINDOW w_c AS (PARTITION BY pt_from, pt_to ORDER BY lin_id, loc_from)
;
You can load this query in QGIS via DB Manager or by creating a VIEW
in PostgreSQL.
Then you simply need to use offset_nr
as a data-defined offset in QGIS.
The result looks like that:
Please note that the lines may not always connect correctly at intersection points. In depends on the order used to calculate offset_nr
.
Edit: to generate the points to be used for intersections, use a query like that:
WITH ix AS
(
SELECT DISTINCT ST_Intersection(ST_SnapToGrid(a.geom, 50.0), ST_SnapToGrid(b.geom, 50.0)) geom
FROM lin a JOIN lin b ON ST_Intersects(a.geom,b.geom)
),
ix_simple_lines AS
(
SELECT
(ST_Dump(ST_LineMerge(ST_CollectionExtract(geom, 2)))).geom geom
FROM
ix
),
ix_points AS
(
SELECT
(ST_Dump(ST_CollectionExtract(geom, 1))).geom geom
FROM
ix
)
SELECT
row_number() OVER() id,
geom
FROM
(
SELECT ST_StartPoint(geom) geom FROM ix_simple_lines
UNION
SELECT ST_EndPoint(geom) FROM ix_simple_lines
UNION
SELECT geom FROM ix_points
) points_union
Best Answer
You should not bother about using &&. That functionality is built in for example ST_Contains.
You should also avoid using buffer in this way.
What is it really that you want to get. Now you are selecting the linestrings where s1 is totally contained in the buffered line from r1. But from you write, what you want is all lines in s1 that is intersecting with the line(s) in r1.
What I don't understand is also why the lines is in different tables. What is the difference between the lines in s1 and r1.
aha, now I think I understand. You want the parts of the linestrings in s1 tat is shared with the lines in r1. Then you want the intersecting part of the linestring, not just all linestrings that is intersecting. s that what you want.
Ok.
If you, as I thought first, want all lines in s1 that intersects with lines in r1 or comes closer than 0.2*(length on line) than you can do:
That will use any gist indexes present. But why do you want a longer line to catch lines wider than a short line?
If, what you want is the intersecting part of the line in s1 and r1, then you should use ST_Intersection together with ST_Intersects to only put work on lines intersecting. In this case you will also have to use ST_Buffer if you want the calculation to be forgiving about inaccuracy.
HTH Nicklas