Clustering similar lines in PostGIS?

Question

In PostGIS I have a large database of linestrings which I have obtained using ST_makeline function.

I would like to cluster similar lines based on where the lines start and end.

Here is an example of what I would like to obtain:

From this post Grouping lines together to form one line after a distance threshold is met I know that in ArcGIS Desktop there is a tool that identifies the mean direction, length, and geographic center for a set of lines.

Do you know if there is an equivalent tool in QGIS?

Do you think it would be easier to cluster start and end coordinates instead of lines?

---

Basing on this post How to merge connected lines with same direction (PostGis), I wrote the following query that seems to work although it is very slow:

select a.routeid as route1, b.routeid as route2,a.the_geom as geom1, b.the_geom as geom2 from prova a, prova b 
WHERE (ST_DWithin(st_startpoint(a.the_geom),st_startpoint(b.the_geom),0.1) and ST_DWithin(st_endpoint(a.the_geom),st_endpoint(b.the_geom),0.1) and (a.routeid <> b.routeid))

---

To prevent a line from being included in more than one separate cluster, I used distinct on as follows where line1 and line2 are the same tables as prova:

 SELECT * FROM(
 SELECT DISTINCT ON (line2.routeid) line1.routeid as route1, line2.routeid as route2,line1.the_geom, line2.the_geom  
 FROM line2 INNER JOIN line1 ON
 ST_DWithin(st_startpoint(line1.the_geom),st_startpoint(line2.the_geom),0.1) 
 AND ST_DWithin(st_endpoint(line1.the_geom),  st_endpoint(line2.the_geom),  0.1)
 ORDER BY line2.routeid) as clusters
 ORDER BY route1;

Answer 1

An easy way to cut this query's work in half is to use a.routeid < b.routeid instead of a.routeid <> b.routeid. This prevents PostGIS from making the same line comparisons twice.

Another change that may be helpful is to add a spatial index operator. With the query as is, I get the following plan:

Nested Loop  (cost=0.00..58.73 rows=1 width=72)
  Join Filter: ((a.routeid <> b.routeid) AND (st_startpoint(a.the_geom) && st_expand(st_startpoint(b.the_geom), 0.1::double precision)) AND (st_startpoint(b.the_geom) && st_expand(st_startpoint(a.the_geom), 0.1::double precision)) AND (st_endpoint(a.the_ge (...)
  ->  Seq Scan on prova a  (cost=0.00..1.10 rows=10 width=36)
  ->  Materialize  (cost=0.00..1.15 rows=10 width=36)
        ->  Seq Scan on prova b  (cost=0.00..1.10 rows=10 width=36)

I'm using dummy data, so your plan may differ. But the nested loop of sequential scans can means lots of comparisons, and lots of disk reads, especially as the number of rows in the table increases.

If no lines are vertical or horizontal, you could use the && operator to force a bounding box intersection:

EXPLAIN
SELECT a.routeid as route1, b.routeid as route2,a.the_geom as geom1, b.the_geom as geom2       from prova a, prova b 
WHERE ST_DWithin(st_startpoint(a.the_geom),st_startpoint(b.the_geom),0.1) 
  AND ST_DWithin(st_endpoint(a.the_geom),  st_endpoint(b.the_geom),  0.1)
  AND a.the_geom && b.the_geom
  AND a.routeid <> b.routeid

On my test set with only 10 rows, this still causes a switch to an inner-loop index scan:

Nested Loop  (cost=0.00..17.43 rows=1 width=72)
  Join Filter: ((a.routeid <> b.routeid) AND (st_startpoint(a.the_geom) && st_expand(st_startpoint(b.the_geom), 0.1::double precision)) AND (st_startpoint(b.the_geom) && st_expand(st_startpoint(a.the_geom), 0.1::double precision)) AND (st_endpoint(a.the_ge (...)
  ->  Seq Scan on prova a  (cost=0.00..1.10 rows=10 width=36)
  ->  Index Scan using prova_the_geom_idx on prova b  (cost=0.00..1.07 rows=1 width=36)
        Index Cond: (a.the_geom && the_geom)

If there's a possibility that the lines could be similar without their bounding boxes intersecting (vertical or horizontal lines), you can expand the original bounding box by your threshold distance using ST_Expand:

SELECT a.routeid as route1, b.routeid as route2,a.the_geom as geom1, b.the_geom as geom2       from prova a, prova b 
WHERE ST_DWithin(st_startpoint(a.the_geom),st_startpoint(b.the_geom),0.1) 
  AND ST_DWithin(st_endpoint(a.the_geom),  st_endpoint(b.the_geom),  0.1)
  AND ST_Expand(a.the_geom, 0.1) && b.the_geom
  AND a.routeid <> b.routeid

An alternative, using the bounding box distance operator <#>, causes a switch back to the original plan on my test data. You might get different results on your actual data.

If these queries are still too slow, you might try extracting the endpoints up front and save those coordinates in their own geometry columns, which could then be indexed.

Nothing about this approach prevents a line from being included in more than one separate cluster. If that is a problem, you'll need to add in some additional logic to handle that (see PostgreSQL's DISTINCT ON keyword for some ideas).