ST_Polygonize creating single polygon

Question

I have a situation that is very similar to the issue described here. However, the accepted answer does not work for me.

I am essentially attempting to create polygons from contour lines generated using ST_Contour.

Here is my query:

WITH inputs AS (
    SELECT 10 AS resample_factor
),
bbox AS (
    SELECT
        ST_Transform(
            ST_MakeEnvelope(
                -93.1585693359375,
                46.71915170604123,
                -92.24121093750001,
                47.28109232255293,
                4326
            ),
            3857
        ) AS geom
),
-- clip raster to AOI
clip_rast AS (
    SELECT
        ST_Union(
            ST_Clip(rast, 1, geom, TRUE)
        ) as rast
    FROM _l_d30_cc_2060_standage_orig_cog
    JOIN bbox ON ST_Intersects(rast, geom)
),
-- Resample raster 
resampled_rast AS (
    SELECT ST_Resample(
        rast, 
        (ST_PixelWidth(rast) * resample_factor),
        (ST_PixelHeight(rast) * resample_factor), 
        algorithm => 'Cubic'
    ) AS rast
    FROM clip_rast, inputs
),
-- Generate contours
contours AS (
    SELECT (ST_Contour(
        rast, 
        bandnumber => 1, 
        level_interval => 10
    )).*
    FROM resampled_rast
),
rast_outline AS (
    SELECT ST_ExteriorRing(ST_ConcaveHull(
    (ST_DumpAsPolygons(ST_Reclass(rast, 1,
    '(0-1000):1',
    '16BUI', -1))).geom, 0)) AS geom
--  SELECT ST_Polygon(rast) AS geom
    FROM resampled_rast
),
-- Create lines from combination of raster boundary and contours
all_lines AS (
    SELECT
        value,
        (ST_Dump(
            ST_Node(
                ST_Collect(
                    ro.geom,                
                    c.geom
                )
            )
        )).geom AS geom
    FROM contours AS c, rast_outline AS ro
),
-- Convert lines to polygons
all_polys AS (
    SELECT
        value,
        (ST_Dump(ST_Polygonize(geom))).geom AS geom
    FROM all_lines
    GROUP BY value
)
SELECT value, geom
FROM all_polys;

When I select all_polys, I see this in pgAdmin's geometry viewer (a single polygon):

However, when I select all_lines I see this:

This is the result I am after (if these were polygons instead of linestrings).

I tried running a negative buffer on the exterior rast_outline to ensure that the contour lines overlapped with it, but that didn't help.

Would using ST_Buffer for all_lines do it? – user30184 Oct 27 '23 at 05:52 — user30184, Oct 27 '23 at 05:52

score 3 · Answer 1 · answered Oct 27 '23 at 20:29

I'm going to put down my eventual solution here just in case it helps anyone in the future.

I was after polygon contours, however polygonize => true was not working for me. I was getting all zeroes for values. However, the line output was working fine. So I needed to combine the line values with the polygon geometry and this was my query, in the end:

WITH inputs AS (
    SELECT 10 AS resample_factor
),
bbox AS (
    SELECT
        ST_Transform(
            ST_MakeEnvelope(
                -93.6585693359375,
                47.21915170604123,
                -92.74121093750001,
                47.78109232255293,
                4326
            ),
            3857
        ) AS geom
),
-- clip raster to AOI
clip_rast AS (
    SELECT
        ST_Union(
            ST_Clip(rast, 1, geom, TRUE)
        ) as rast
    FROM _l_d30_cc_2060_standage_orig_cog
    JOIN bbox ON ST_Intersects(rast, geom)
),
-- Resample raster 
resampled_rast AS (
    SELECT ST_Resample(
        rast, 
        (ST_PixelWidth(rast) * resample_factor),
        (ST_PixelHeight(rast) * resample_factor), 
        algorithm => 'Cubic'
    ) AS rast
    FROM clip_rast, inputs
),
-- Generate contours
contours AS (
    SELECT (ST_Contour(
        rast, 
        bandnumber => 1, 
        level_interval => 10,
        level_base => 10
    )).*
    FROM resampled_rast
),
-- Generate contours as polygons
poly_contys AS (
    SELECT (ST_Contour(
        rast, 
        bandnumber => 1, 
        level_interval => 10,
        level_base => 10,
        polygonize => true
    )).*
    FROM resampled_rast
),
-- join values to polygons
poly_contys_with_vals AS (
    SELECT c.id, c.value, pc.geom
    FROM contours AS c, poly_contys AS pc
    WHERE ST_CoveredBy(c.geom, pc.geom)
)
SELECT DISTINCT ON (c.geom) c.value, c.geom, row_number() OVER () AS id
FROM poly_contys_with_vals AS c

This gives the result I want:

Any suggestion for how to speed it up or other approaches are welcome.

jbalk · Accepted Answer · 2023-10-27T18:40:22.350

2

This part of your query is causing every contour geom to be joined with the raster boundary:

ST_Collect(
  ro.geom,                
  c.geom
)

Change that to a table union:

all_lines AS (
    SELECT
        value,
        (ST_Dump(
            ST_Node(
                    geom               
            )
        )).geom AS geom
    FROM contours
union
select 9999 value, geom from rast_outline
)

Now st_polygonize should work as expected. You can also use st_makepolygon, but you have to make sure the linestring is closed:

select value, st_setsrid(st_makepolygon(geom),3857) as geom from 
(select value, 
  case when not st_isclosed(geom) then st_addpoint(geom,st_startpoint(geom))
  else geom end geom
  from lines
) foo;

edited Oct 27 '23 at 18:40

answered Oct 27 '23 at 17:10

jbalk

7,505
1
17
39

Thanks for the help, @jbalk. This answer is correct is the sense in that it solves the problem of ST_Polygonize producing a single polygon. It's not the complete answer for my situation because every contour line that isn't closed is being combined with the exterior (or background) polygon of the exterior ring, and those values are being replaced by the 9999 assigned to it. But I suppose this next issue deserves its own separate question. – rumski20 Oct 27 '23 at 19:36

ST_Polygonize creating single polygon

2 Answers2