Last active
January 11, 2021 16:44
-
-
Save TimMcCauley/064603b6a41517fb60a7cf80862d33b2 to your computer and use it in GitHub Desktop.
Groups parallel lines using DBScan into clusters and unions their buffers
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
-- 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; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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.