I have OpenStreetMap data, which represents larger roads (e.g. dual carriageways, aka divided highways) as parallel lines.

The question is answered for ArcGIS or Skeletron here, for example using ESRI's Merge Divided Roads or Collapse Dual Lines to Centerline tools.

Does anyone have suggestions for how to do this within PostGIS? I am working on a server without ArcGIS, and would like to avoid the overhead of transferring the data back and forth across machines.

Update The linked question isn't really answered, and is looking for a solution that is not PostGIS. This should be kept open.


2 Answers 2


Maybe it's not a very elegant solution and it's based on the method of answering question:https://gis.stackexchange.com/a/295348/120129,

nevertheless, it's one of the variants of solving your question, set up a geoinstrument (for me it's a road with the name "road_border" 11 m wide, type - line (MultiLineString)),

run a simple geoinstrument :-),

      tbla AS (SELECT id, ((ST_DumpPoints(geom)).geom) geom FROM road_border),
      tblb AS (SELECT (ST_Buffer((ST_Dump(geom)).geom, 0.0001)) geom FROM tbla),
      tblc AS (SELECT ST_Centroid(ST_UnaryUnion(Unnest(ST_ClusterWithin(geom, 0.0001)))) geom FROM tblb),
      tbld AS (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblc),
      tble AS (SELECT (ST_Boundary(geom)) geom FROM ST_Dump((SELECT ST_Polygonize(geom) geom FROM (SELECT ST_Union(geom) geom FROM
      (SELECT ST_ExteriorRing(geom) geom FROM tbld) AS lines) AS foo))),
      tblf AS (SELECT ST_MakeLine(p1, p2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) p1,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom))) p2 FROM tble) AS geom),
      tblg AS (SELECT ST_Buffer((ST_Dump(ST_Union(ST_Buffer(geom, 0.0005)))).geom, -0.0005) geom FROM road_border)
      SELECT ST_Intersection (a.geom, b.geom) geom FROM tblf a JOIN tblg b ON ST_Within (a.geom, b.geom);

and see the result.

This script is called - ST_RoadAxisFromDelaunayTriangulation...

One more way.

The input is a simple road in the form of parallel lines as shown in Figure 1.

enter image description here

Figure 1.

Create GeoSQL function:

    geom GEOMETRY,
    n integer)
      geodata AS (SELECT (ST_Dump(geom)).geom),
      left_line AS (SELECT (ST_Dump(geom)).geom geom FROM geodata LIMIT 1),
      right_line AS (SELECT (a.geom) geom FROM geodata a WHERE NOT EXISTS (SELECT 1 FROM left_line b WHERE ST_Intersects(a.geom, b.geom))),
      intervals AS (SELECT generate_series (0, n) steps),
      lineinterpnt AS (SELECT steps AS stp, ST_LineInterpolatePoint(geom, steps/(SELECT count(steps)::float-1 FROM intervals)) geom FROM left_line, intervals GROUP BY intervals.steps, geom ORDER BY intervals.steps),
      sleepers AS (SELECT DISTINCT stp, ST_ShortestLine(a.geom, b.geom) geom FROM lineinterpnt a CROSS JOIN LATERAL 
            (SELECT (geom) geom FROM right_line) b ORDER BY stp)
        SELECT ST_MakeLine(geom) geom FROM (SELECT stp, ST_LineInterpolatePoint(geom, 0.5) geom FROM sleepers ORDER BY stp) foo;

Run GeoSQL function:

SELECT ST_RoadAxis(ST_Union(geom), 300) geom FROM <parall_line_name_table>

The result should be a centerline of the road as shown in Figure 2.

enter image description here

Figure 2.

The joint image is shown in Figure 3. enter image description here

Figure 3.

This script is called - ST_RoadAxis.

(-: FOGS :-)...


If you don't want to roll your own function and can put up with a little bit of innacuracy, this might possibly work.

Running ST_ConcaveHull i.e. ST_ConcaveHull(ST_Collect(geom),0.xx) (99% seemed to work as a parameter for my very limited test), then (so long as you have CGAL and SFGCAL compiled or linked to your binary ST_ApproximateMedialAxis would appear to be one way of doing it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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