Wednesday, 10 April 2019

postgis - How to find linestrings which aren't joined


Given a PostGIS table with a column of linestrings, which are expected to all be joined (they represent a network), how can I query the table to determine whether or not they are all joined (i.e. it should be possible to reach any linestring from any other, by traversing the network)?


I was thinking of creating some kind of aggregate geometry and comparing some measurement of that aggregate geometry (e.g. total "length") with the sum of all lengths, but haven't been able to get anywhere.



Answer



Assuming your dataset is small enough to fit in memory, you can cluster all of your geometries and then examine the number of clusters and the number of lines in each cluster.


If you want to use intersection as the basis of testing, you can use ST_ClusterIntersecting.


If you want to use use a distance tolerance instead, you can use ST_CusterWithin or ST_ClusterDBSCAN (which has a better interface).


Here's an example using ST_ClusterDBSCAN:



SELECT
cluster_id,
count(*) AS num_lines,
array_agg(feature_id) AS features
FROM
(SELECT
feature_id,
ST_ClusterDBSCAN(geom, eps := 1e-8) OVER() AS cluster_id
FROM lines) sq
ORDER BY num_lines DESC

OFFSET 1

In this case, we're linking our line geometries with a tolerance of 1e-8, and then assigning them to a cluster_id. If our network is entirely connected, then we'll only have one cluster ID. Assuming this is not the case, we count the number of features in each cluster, and then list out the feature IDs of every feature in each non-dominant cluster. For each non-dominant cluster, one linkage needs to be established to the dominant cluster in order to build a fully connected network.


One last thing: if you wanted to verify that your network was entirely connected by noded intersections (that is, lines that cross without a node would not be considered to intersect), you could do this by running your cluster function on ST_Points(geom) instead of geom.


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...