In order to perform the basic skeletonization, we'll calculate the Voronoi polygons on the nodes that make up the original stream polygon. By default, the edges of the Voronoi polygons find the line that demarcates the midpoint between points. We will leverage this tendency by treating our lines like points—adding extra points to the lines and then converting the lines to a point set. This approach, in combination with the Voronoi approach, will provide an initial estimate of the polygon's centerline.

We will add extra points to our input geometries using the ST_Segmentize function and then convert the geometries to points using ST_DumpPoints:

CREATE TABLE chp06.voronoi_points AS(
  SELECT (ST_DumpPoints(ST_Segmentize(the_geom, 5))).geom AS the_geom 
    FROM chp06.voronoi_hydro
  UNION ALL
  SELECT (ST_DumpPoints(ST_Extent(the_geom))).geom AS the_geom 
    FROM chp06.voronoi_hydro
)

The following screenshot shows our polygons as a set of points if we view it on a desktop GIS:

The set of points in the preceding screenshot is what we feed into our Voronoi calculation:

CREATE TABLE chp06.voronoi AS(
  SELECT (ST_Dump(
    ST_SetSRID(
      ST_VoronoiPolygons(points.the_geom),
      3734))).geom as the_geom
FROM (SELECT ST_Collect(ST_SetSRID(the_geom, 3734)) as the_geom FROM chp06.voronoi_points) as points);

The following screenshot shows a Voronoi diagram derived from our points:

If you look closely at the preceding screenshot, you will see the basic centerline displayed in our new data. Now we will take the first step toward extracting it. We should index our inputs and then intersect the Voronoi output with the original stream polygon in order to clean the data back to something reasonable. In the extraction process, we'll also extract the edges from the polygons and remove the edges along the original polygon in order to remove any excess lines before our routing step. This is implemented in the following script:

CREATE INDEX chp06_voronoi_geom_gist
ON chp06.voronoi
USING gist(the_geom);
    
DROP TABLE IF EXISTS voronoi_intersect;
    
CREATE TABLE chp06.voronoi_intersect AS WITH vintersect AS (
  SELECT ST_Intersection(ST_SetSRID(ST_MakeValid(a.the_geom), 3734), 
  ST_MakeValid(b.the_geom)) AS the_geom 
  FROM Chp06.voronoi a, chp06.voronoi_hydro b
  WHERE ST_Intersects(ST_SetSRID(a.the_geom, 3734), b.the_geom)
),
linework AS (
  SELECT chp02.polygon_to_line(the_geom) AS the_geom 
  FROM vintersect
),
polylines AS (
  SELECT ((ST_Dump(ST_Union(lw.the_geom))).geom)
    ::geometry(linestring, 3734) AS the_geom 
  FROM linework AS lw
),
externalbounds AS (
  SELECT chp02.polygon_to_line(the_geom) AS the_geom 
  FROM voronoi_hydro
)
    
SELECT (ST_Dump(ST_Union(p.the_geom))).geom 
  FROM polylines p, externalbounds b 
  WHERE NOT ST_DWithin(p.the_geom, b.the_geom, 5);

Now we have a second-level approximatio of the skeleton (shown in the following screenshot). It is messy, but it starts to highlight the centerline that we seek: