10

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.

Connected Lines

  • 3
    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 Sep 7 '18 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 Sep 8 '18 at 5:25
  • 3
    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 Sep 8 '18 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 Sep 10 '18 at 4:36
  • 2
    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 Sep 10 '18 at 9:22
3

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

Look at the result and share it with others :-). Well, we all landed and I hope you will tell others about this flight and the axial lines of the Cordillera mountain system we saw :-).

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

With you there was a free, cheerful and brave pilot of high aerobatics!

Original solutions.

The script is called: ST_LongestAxisPolygonsFromVoronoiDiagrams.

Warning: Be careful while creating such engines and in your flights! In fact, I had to deal with a lot of technical problems, and it was possible to write more than one scientific paper or dissertation during the circumvention of which. This geoinstrument was created on the basis of my answer: https://gis.stackexchange.com/a/318733/120129, and its modifications, allowing to create such geoinstruments, for example: ST_ShortestAxisPolygonsFromVoronoiDiagrams and others.

(!Not to use for simple shapes like circle, triangle, rectangle!)

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

I suppose I've managed to cheer someone up!

Have a nice day, everybody.

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;

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
)
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, 888) 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, 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_MakeLine(ST_Centroid(geom)) geom FROM tbll
$BODY$
LANGUAGE SQL;

SELECT ST_LongestAxisPolygonsFromVoronoiDiagrams_gs(geom) 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;

Сreation...

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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