PostGIS functions for generating centerline from polygon

Question

I 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>

Answer 1

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))');

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)

Positive then negative:

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

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

Answer 2

So, if the street were approximately equal to the width as shown in the figure below,

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

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