9

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:

enter image description here

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

enter image description here

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:

enter image description here

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.

  • ST_Split or ST_BuildArea are good candidates for your problem – nickves Apr 26 '16 at 0:11
3

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.

|improve this answer|||||
0

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

|improve this answer|||||
  • That doesn't really answer the question fully. – John Powell Feb 25 '15 at 7:58
0

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

|improve this answer|||||

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.