4

Polygon imageI want to create a centerline from a polygon and I have tried using both ST_StraightSkeleton and ST_ApproximateMedialAxis functions in PostGIS.

It gives an error with a message

straight skeleton of Polygon with touching interior rings is not implemented

and in other instances

CGAL failed to create straightSkeleton

I have tried alternative answers on this platform generating similar errors.

Here is the code

''' 
SELECT
    ogc_fid, 
    subtype, 
    subtypecod,
    ST_ApproximateMedialAxis(geom) as geom
FROM 
    sometable
'''

Is there a proper way to go about this issue?

Sample data <https://www.filehosting.org/file/details/896986/Polygon_5.zip>
14
  • 1
    ST_LongestLine() ? Sep 28 '20 at 23:26
  • Database error: function st_longestline(geometry) does not exist HINT: No function matches the given name and argument types. You might need to add explicit type casts. Thanks but the function did not work.
    – Krismars
    Sep 29 '20 at 0:29
  • 1
    This is what I'm talking about postgis.net/docs/ST_LongestLine.html if you don't have it available that's another issue... Sep 29 '20 at 1:42
  • 3
    You need to clarify what kind of polygons we are talking about, it would be nice if you added a picture in your question, polygons can be simple type of square, rectangle, etc. and complex... Sep 29 '20 at 5:58
  • 2
    try merging the "touching interior rings"
    – Ian Turton
    Sep 29 '20 at 7:24
3

As you can see here, the problem is that CGAL cannot create a skeleton when you have an interior ring that touch the exterior ring (there is a hole in a geom that touches the exterior boundary in one point). You should prepare your geometry first, to go around this issue.

Exemple of geometry like that:

SELECT ST_GeomFromText('POLYGON((0 0, 1 1, 0 2, -1 1, 0 0 , 0.5 0.6, 0 1, -0.5 0.6, 0 0, 0 0))');

exemple of geom

Here is what you can do:

  • If you don't care about the hole, you can try st_convexhull (or st_concavehull) to get only the hull, or even ST_Buffer. You will get only the outside geom (the outside square on the exemple above)
  • If you care, start by doing st_makevalid, and then you can use ST_Buffer to modify the geometry to make sure to remove the place where ther is the problem, either by growing the bridge (buffer positive) or by removing it (buffer negative). You can apply the 2 in a row to make sure that you didn't modified the other points positions

Negative then positive:

st_buffer(st_buffer(st_makevalid(geom), -0.01), 0.01)

negative then positive

Positive then negative:

st_buffer(st_buffer(st_makevalid(geom), 0.01), -0.01)

positive then negative

Then you can make your skeleton (exemple with the last one):

enter image description here

1
  • 1
    I have tried your suggestions and the error message "CGAL failed to create straightSkeleton" still comes up. Thanks
    – Krismars
    Sep 30 '20 at 16:59
1

So, if the street were approximately equal to the width as shown in the figure below, enter image description here

it would be possible to use this customizable geo tool

(for my example, the width of the street is 40 meters),

Configure and run the script:

WITH 
  tbla AS (SELECT ST_ExteriorRing(geom) geom FROM (SELECT ((ST_Dump(geom)).geom) geom FROM street) foo),
  tblb AS (SELECT ((ST_Dump(geom)).geom) geom FROM street),
  tblc AS (SELECT ((ST_DumpPoints(geom)).geom) geom FROM street),
  tbld AS (SELECT ST_Centroid(ST_UnaryUnion(Unnest(ST_ClusterWithin(geom, 60)))) geom FROM tblc),
  tble AS (SELECT ST_ExteriorRing((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tbld),
  tblf 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 tble) AS geom)
  SELECT ST_Intersection (a.geom, b.geom) geom FROM tblf a JOIN tblb b ON ST_Within(a.geom, b.geom)

Result

enter image description here

Delete possible artifacts,

and that's all...

Edit

For your terrible data set, I would have cleaned it up and run a double process, for example:

WITH 
  tbla AS (SELECT ST_ExteriorRing(geom) geom FROM (SELECT ((ST_Dump(geom)).geom) geom FROM "Polygon_5") foo),
  tblb AS (SELECT ((ST_Dump(geom)).geom) geom FROM "Polygon_5"),
  tblc AS (SELECT ((ST_DumpPoints(geom)).geom) geom FROM "Polygon_5"),
  tbld AS (SELECT ST_Centroid(ST_UnaryUnion(Unnest(ST_ClusterWithin(geom, 10)))) geom FROM tblc),
  tble AS (SELECT ST_ExteriorRing((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tbld),
  tblf 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 tble) AS geom),
  tblg AS (SELECT ST_Intersection (a.geom, b.geom) geom FROM tblf a JOIN tblb b ON ST_Within(a.geom, b.geom)),
  tblh AS (SELECT ((ST_DumpPoints(geom)).geom) geom FROM tblg),
  tbli AS (SELECT ST_ExteriorRing((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblh),
  tblj 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 tbli) AS geom)
  (SELECT ST_Intersection (a.geom, b.geom) geom FROM tblj a JOIN tblb b ON ST_Within(a.geom, b.geom))

Removed the artifacts and a bit of manual refinement...

Good luck in knowing...🙂

2
  • There is another option, to play with an internal one-way buffer... Oct 1 '20 at 6:28
  • There is another option, polygons in a line and ST_VoronoiDiagramsFromLines... Apr 5 at 19:41

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.