How to transform isolines to isopolygons with postgis?

Question

I have a table a postgis table of isolines that's defined like this:

CREATE TABLE myisolines
(
  gid serial NOT NULL,
  isotime timestamp without timezone,
  val numeric(10,4),
  geom geometry(LineString,4326)
);

Visually this linestring objects look like this:

I know the spatial extent of my data, so I can add a Bbox, so the LineStrings can be sort of closed.

I want to create a table of isopolygons myisopolygons from myisolines table, with polygons, that won't overlap but create a continuous surface and have a column val with the lowest val of isolines, from which the polygon was formed. I understand it can be formed from self-closed isoline (island), or isoline closed with bbox, in that case the val should be taken from that particular isoline. Visually it should look like this:

I thought that I can create topology somehow and then transform faces to polygons, but I don't understand how to do it properly. How can this be done?

Another option would be to recursively use a difference function between bbox and each created polygon, but I guess that's not the right way to do it, and definitely not fast at all.

Answer 1

Here's a solution using ST_Polygonize. It generates a polygon for each boundary, and provides the min and max elevation covered by the polygon. The algorithm cannot distinguish between a peak and a depression and will return the same elevation for both the min and max in these cases.

WITH closed_contours AS (
    SELECT 
      ST_Union(geom) AS geom 
    FROM 
      (SELECT geom FROM contours 
       UNION ALL 
       SELECT ST_SetSRID(ST_Boundary(ST_Expand(ST_Extent(geom), -1e-10)), 4326) 
       FROM contours) 
sq)

SELECT
  poly_id, 
  min(polys.geom) AS geom, 
  min(elevation)  AS min_elev, 
  max(elevation)  AS max_elev
FROM
  (SELECT row_number() OVER () AS poly_id, geom FROM
      (SELECT 
         (ST_Dump(ST_Polygonize(geom))).geom
       FROM closed_contours) dump
  ) polys
INNER JOIN contours ON ST_Intersects(polys.geom, contours.geom)
GROUP BY poly_id;

The WITH clause of the query "closes" any open contours by unioning them with the slightly-contracted extent of the existing contours. (The extent is contracted to wash over any round-off errors resulting from the use of ST_Extent, which produces a single-precision box, with ST_Polygonize, which requires perfectly closed and noded inputs in doulbe-precision). If your contours are already closed (i.e., you're working with an island), then this step can be omitted.

Answer 2

I am not very exerienced, but I would try function geometry ST_MakePolygon(geometry outerlinestring, geometry[] interiorlinestrings);

Answer 3

Using your bbox and iterating over each contour line, you can use ST_ConcaveHull to convert each region into a polygon.