Skip to content

Instantly share code, notes, and snippets.

@TimMcCauley
Last active January 11, 2021 16:44
Show Gist options
  • Save TimMcCauley/064603b6a41517fb60a7cf80862d33b2 to your computer and use it in GitHub Desktop.
Save TimMcCauley/064603b6a41517fb60a7cf80862d33b2 to your computer and use it in GitHub Desktop.
Groups parallel lines using DBScan into clusters and unions their buffers
-- This set of queries derives a perpendicular line through all input lines
DROP TABLE IF EXISTS intersection_points;
CREATE TABLE intersection_points
AS
-- 2 lines will suffice to compute a perpendicular
WITH nbr_parallel_lines AS (
SELECT A.geom AS A_geom, B.geom AS B_geom
FROM parallel_lines AS A, parallel_lines AS B
WHERE A.fid_1 <> B.fid_1
ORDER BY A.geom <-> B.geom
LIMIT 1
),
-- Take the first and find the nearest point on the second which will be the perpendicular
perpendicular_line AS (
SELECT ST_MakeLine(ST_Centroid(A_geom), ST_ClosestPoint(B_geom, ST_Centroid(A_geom))) as geom
FROM nbr_parallel_lines
),
-- Extend the line to make it intersect all input lines which are parallel
perpendicular_line_extended AS (
SELECT ST_MakeLine(ST_TRANSLATE(a, sin(az1) * len, cos(az1) *
len),ST_TRANSLATE(b,sin(az2) * len, cos(az2) * len)) as the_geom
FROM (
SELECT a, b, ST_Azimuth(a,b) AS az1, ST_Azimuth(b, a) AS az2, ST_Distance(a,b) + 100 AS len
FROM (
SELECT ST_StartPoint(geom) AS a, ST_EndPoint(geom) AS b
FROM perpendicular_line
) AS sub
)
AS sub2)
-- Find all intersecting points
SELECT all_lines.fid_1, ST_Intersection(perp_line.the_geom, all_lines.geom) geom
FROM perpendicular_line_extended AS perp_line, parallel_lines AS all_lines;
-- This set of queries uses DBScan to find the cluster of intersecting points
-- and joins them with the input data which then are buffered
DROP TABLE IF EXISTS buffers;
CREATE TABLE buffers
AS
-- Determine clusters using a user provided epsilon and min points
WITH clusters AS (
SELECT pl.geom as linestring_geom,
ipts.geom as pt_geom, ipts.fid_1,
ST_ClusterDBSCAN(ipts.geom, eps := 0.3, minpoints := 5) OVER() AS clst_id
FROM intersection_points AS ipts
JOIN parallel_lines AS pl
ON ipts.fid_1 = pl.fid_1
),
-- For each cluster found determine the maximum distance between the points and divide it by
-- a user given input;
-- this we will use as input for our buffering
clusters_distances AS (
SELECT
a.clst_id, MAX(ST_Distance(a.pt_geom, b.pt_geom))/2 AS dist
FROM
clusters AS a
CROSS JOIN
clusters AS b
WHERE a.clst_id = b.clst_id
GROUP BY a.clst_id
),
clusters_enriched AS (
SELECT a.*, b.dist FROM clusters AS a
JOIN clusters_distances AS b
ON a.clst_id = b.clst_id
)
-- Union all clustered linestrings after buffering them
SELECT ST_Union(ST_Buffer(cl.linestring_geom, cl.dist)) AS geom
FROM clusters_enriched AS cl
GROUP BY cl.clst_id;
-- Last but not least compute the approximate medial axis which
-- will yield the center line
DROP TABLE IF EXISTS buffers_medial_axes;
CREATE TABLE buffers_medial_axes AS
SELECT ST_ApproximateMedialAxis(geom) as medial_axis FROM buffers;
@TimMcCauley
Copy link
Author

TimMcCauley commented Jan 11, 2021

This query is a way to find what the OP asked for in https://gis.stackexchange.com/questions/383999/how-to-calculate-buffer-distance-for-different-groups-of-parallel-lines-in-qgis - however using PostGIS.

  1. It first determines the clusters using DBScan given a perpendicular between 2 parallel lines which intersects all input lines creating points which can be used as an input for unsupervised learning:

Screenshot 2021-01-11 at 15 06 42

  1. Afterwards it buffers the clusters and computes the medial axis:

Screenshot 2021-01-11 at 17 36 56

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment