Selecting the longest connected lines from PostGIS ST_ApproximateMedialAxis

Question

I used PostGIS ST_ApproximateMedialAxis to transform polygons representing mountain ranges into a series of lines. Now, I'd like to select only the "longest" line segment for each group, so that I can dissolve the lines and use it for "curved" labeling. Can I extend the query in PostGIS to end up with the single series of arcs?

The end result should look something like the lines highlighted in yellow in the picture below.

the ultimate solution to this (and sweet sugar to many more of PostGIS' functions and functionality) is building and enabling the SFCGAL backend for PostGIS (with default GEOS, the medial axis is directly derived off the ST_StraightSkeleton)...but that can be a massive pain... — geozelot, Commented Sep 7, 2018 at 7:22
I'm using Postgres 9.5 and PostGIS 2.3, and have installed the postgis_sfcgal extension, thus allowing me to use ST_ApproximateMedialAxis. I chose Medial Axis results because the StraightSkeleton had so many more "fishbones" coming in from the sides of the polygons — RyanKDalton, Commented Sep 8, 2018 at 5:25
What do you mean by each group? Do you have a groupid or somesuch? Perhaps you should post the SQL you have run so far. — John Powell, Commented Sep 8, 2018 at 20:54
Instead of"group", perhaps I could have used "set", as in "set of connected arcs". For example, in the picture there are 3 "sets", and the longest axis of arcs are highlighted in yellow. — RyanKDalton, Commented Sep 10, 2018 at 4:36
In that case, you probably want to use one of the clustering functions, DBScan, for example, to create the sets/groups, and then do a max(length) ... group by set kind of thing. It is a window function, so you can get a group ID out of it quite naturally. — John Powell, Commented Sep 10, 2018 at 9:22

Cyril Mikhalchenko · Accepted Answer · 2022-05-19 19:37:39Z

I am publishing my answer for those who visit this page hoping to find an answer to the problem you have formulated!

So, I invite everyone to fly with me in a plane with the engine I have assembled with the name: "Free and Magic geoSQL", in order to admire the beautiful mountain scenery of the Cordillera mountain system and from a height of flight to see their centerlines :-)! We have a difficult and interesting task ahead of us, but we can do it even without SFCGAL library and ST_ApproximateMedialAxis and ST_StraightSkeleton functions, because we are pilots and navigators of high class :-), I strongly hope so! I command: "Fasten your seatbelts", "From the screw", "Fly"! and of course I start the engine, which works so far for only one polygon:

WITH
tbla AS (SELECT (ST_Dump(geom)).geom FROM <name_polygon_table>),
tblb AS (SELECT (ST_DumpPoints(geom)).geom FROM tbla),
tblc AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tblb),
tbld AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom FROM tblc a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
tble AS (SELECT line1, line2, CASE WHEN line1>line2 THEN line1 WHEN line1<=line2 THEN line2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
(ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS line1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS line2 FROM tbla) AS foo),
tblf AS (SELECT (ST_Buffer(line2, ST_Length(line2), 'endcap=flat join=round')) geom FROM tble),
tblg AS (SELECT DISTINCT (ST_Dump(ST_Difference(a.geom, b.geom))).geom FROM tblf a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
tblh AS (SELECT (ST_Dump(ST_Buffer(geom, -0.000001))).geom geom FROM tblg),
tbli AS (SELECT (ST_Buffer(geom, 0.000002)) geom FROM tblh a JOIN tble b ON ST_Intersects(a.geom, b.line2)),
tblj AS (SELECT ST_ExteriorRing(ST_Union(a.geom)) geom FROM tbld a JOIN tbli b ON ST_Intersects(a.geom, b.geom)),
tblk AS (SELECT ST_MakeLine(pt1, pt2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) pt1, ST_PointN(geom, generate_series(2, ST_NPoints(geom))) pt2 FROM tblj) AS geom),
tbll AS (SELECT (ST_Dump(ST_Intersection(ST_Buffer(a.geom, -0.0001), b.geom))).geom FROM tbla a JOIN tblk b ON ST_Contains(a.geom, b.geom))
SELECT ST_Union(geom) geom FROM tbll

(If necessary, adjust it to the size of your polygons (function ST_Buffer), I played with data in EPSG:4326).

From now on, everyone will be able to see and create the centerline of the polygon.

Original solutions.

The script is called: ST_LongestAxisPolygonsFromVoronoiDiagrams.

(!Not to use for simple shapes like rectangle!)

"I want to thank all those who contribute to the creation and development of promising geoinstruments by sharing their ideas and codes..."

First editing:

CREATE OR REPLACE FUNCTION ST_LongestAxisPolygonsFromVoronoiDiagrams(
    geom GEOMETRY
)
RETURNS GEOMETRY AS
$BODY$
WITH
tbla AS (SELECT (ST_Dump(geom)).geom),
tblb AS (SELECT (ST_DumpPoints(geom)).geom FROM tbla),
tblc AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tblb),
tbld AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom FROM tblc a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
tble AS (SELECT line1, line2, CASE WHEN line1>line2 THEN line1 WHEN line1<=line2 THEN line2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
(ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS line1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS line2 FROM tbla) AS foo),
tblf AS (SELECT (ST_Buffer(line2, ST_Length(line2), 'endcap=flat join=round')) geom FROM tble),
tblg AS (SELECT DISTINCT (ST_Dump(ST_Difference(a.geom, b.geom))).geom FROM tblf a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
tblh AS (SELECT (ST_Dump(ST_Buffer(geom, -0.000001))).geom geom FROM tblg),
tbli AS (SELECT (ST_Buffer(geom, 0.000002)) geom FROM tblh a JOIN tble b ON ST_Intersects(a.geom, b.line2)),
tblj AS (SELECT ST_ExteriorRing(ST_Union(a.geom)) geom FROM tbld a JOIN tbli b ON ST_Intersects(a.geom, b.geom)),
tblk AS (SELECT ST_MakeLine(pt1, pt2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) pt1, ST_PointN(geom, generate_series(2, ST_NPoints(geom))) pt2 FROM tblj) AS geom),
tbll AS (SELECT (ST_Dump(ST_Intersection(ST_Buffer(a.geom, -0.0001), b.geom))).geom FROM tbla a JOIN tblk b ON ST_Contains(a.geom, b.geom))
SELECT ST_Union(geom) geom FROM tbll
$BODY$
LANGUAGE SQL;

RUN

SELECT ST_LongestAxisPolygonsFromVoronoiDiagrams(geom) geom FROM <source_table>

If necessary, make the line smoother, add the function generate_series...

Second edition 13.03.2020

CREATE OR REPLACE FUNCTION ST_LongestAxisPolygonsFromVoronoiDiagrams_gs(
    geom GEOMETRY,
    n integer)
RETURNS GEOMETRY AS
$BODY$
WITH
    tbla AS (SELECT (ST_Dump(geom)).geom),
    tblb AS (WITH btbl AS (SELECT (ST_ExteriorRing(((ST_Dump(geom)).geom))) geom),
            intervals AS (SELECT generate_series (0, n) as steps)
            SELECT steps AS stp, ST_LineInterpolatePoint(geom, steps/(SELECT count(steps)::float-1 FROM intervals)) geom FROM btbl, intervals GROUP BY intervals.steps, geom),
    tblc AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tblb),
    tbld AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom FROM tblc a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
    tble AS (SELECT line1, line2 FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
        (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS line1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS line2 FROM tbla) AS foo),
    tblf AS (SELECT (ST_Buffer(line2, ST_Length(line2), 'endcap=flat join=round')) geom FROM tble),
    tblg AS (SELECT DISTINCT (ST_Dump(ST_Difference(a.geom, b.geom))).geom FROM tblf a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
    tblh AS (SELECT (ST_Dump(ST_Buffer(geom, -0.000001))).geom geom FROM tblg),
    tbli AS (SELECT (ST_Buffer(geom, 0.000002)) geom FROM tblh a JOIN tble b ON ST_Intersects(a.geom, b.line2)),
    tblj AS (SELECT ST_ExteriorRing(ST_Union(a.geom)) geom FROM tbld a JOIN tbli b ON ST_Intersects(a.geom, b.geom)),
    tblk AS (SELECT ST_MakeLine(pt1, pt2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) pt1, ST_PointN(geom, generate_series(2, ST_NPoints(geom))) pt2 FROM tblj) AS geom),
    tbll AS (SELECT (ST_Dump(ST_Intersection(ST_Buffer(a.geom, -0.0001), b.geom))).geom FROM tbla a JOIN tblk b ON ST_Contains(a.geom, b.geom))
            SELECT ST_MakeLine(ST_Centroid(geom)) geom FROM tbll
$BODY$
LANGUAGE SQL;

RUN

SELECT ST_LongestAxisPolygonsFromVoronoiDiagrams_gs(geom, 300) geom FROM <source_table>

Third edition 30.03.2020

As an option, for 2 or more polygons run the following script:

SELECT DISTINCT b.geom
FROM <source_table> AS a
JOIN LATERAL (SELECT ST_LongestAxisPolygonsFromVoronoiDiagrams_gs(geom) AS geom
FROM <source_table>) AS b ON true;

Fourth edition 13.05.2021

CREATE OR REPLACE FUNCTION ST_CentraLAxisLongestLineFromPolygons_VD(
geom GEOMETRY,
n integer)
RETURNS GEOMETRY AS
$BODY$
WITH
    tbla AS (SELECT (ST_Dump(geom)).geom geom),
    tblb AS (SELECT ST_Buffer(ST_Boundary(ST_OrientedEnvelope(geom)),0.00001) geom FROM tbla),
    tblc AS (SELECT (ST_Dump(ST_Difference(ST_Boundary(a.geom), b.geom))).geom geom FROM tbla a JOIN tblb b ON true),
    tbld AS (SELECT (ST_Dump(ST_Difference(ST_Buffer(ST_OrientedEnvelope(a.geom), -0.00001), b.geom))).geom geom FROM tbla a JOIN tbla b ON true),
    tble AS (SELECT ST_Buffer(geom, 0.00001) geom  FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
        (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),4))) AS geom FROM tbla) AS foo),
    tblf AS (SELECT generate_series (0, n) as steps),
    tblg AS (SELECT steps AS stp, ST_LineInterpolatePoint(geom, steps/(SELECT count(steps)::float-1 FROM tblf)) geom FROM tblc, tblf GROUP BY tblf.steps, geom),
    tblh AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tblg),
    tbli AS (SELECT (a.geom) geom FROM tbld a JOIN tble b ON ST_Intersects(a.geom,  ST_Buffer(b.geom, 0.00001))),
    tblj AS (SELECT ST_ExteriorRing(ST_Union(a.geom)) geom FROM tblh a JOIN tbli b ON ST_Intersects(a.geom, b.geom))
            SELECT ST_Intersection(a.geom, b.geom) geom FROM tblj a JOIN tbla b ON ST_Intersects(a.geom, b.geom);
$BODY$
LANGUAGE SQL;

Run

SELECT ST_CentraLAxisLongestLineFromPolygons_VD(geom, 100) AS geom FROM <source_table>

Stack Exchange Network

Selecting the longest connected lines from PostGIS ST_ApproximateMedialAxis

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
postgis
sql
topology
dissolve
connectivity-analysis
or ask your own question.

Linked

Hot Network Questions

Selecting the longest connected lines from PostGIS ST_ApproximateMedialAxis

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged postgissqltopologydissolveconnectivity-analysis or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
postgis
sql
topology
dissolve
connectivity-analysis
or ask your own question.