1

So,

I wrote a recursive query that works well, maybe only for me :-),

[https://gis.stackexchange.com/a/333562/120129]

WITH 
    tbla AS (SELECT (ST_DumpPoints(geom)).geom FROM adm_polygons),
    tblb AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tbla),
    tblc AS (SELECT ST_Intersection(a.geom, b.geom) geom FROM tblb a JOIN poly_extent b ON ST_Intersects(a.geom, b.geom)),
    tbld AS (SELECT ST_Difference(a.geom, b.geom) geom FROM tblc a JOIN adm_polygons b ON ST_Intersects(a.geom, b.geom) AND ST_Overlaps(a.geom, b.geom))
    SELECT ST_Union(geom) geom FROM tbld;

See the result in figure,

enter image description here

Then I decided to wrap it in an SQL function:

CREATE OR REPLACE FUNCTION ST_CarvesPolygons(
    geom1 geometry, 
    geom2 geometry)
RETURNS geometry
AS $BODY$
WITH 
    tbla AS (SELECT (ST_DumpPoints(geom1)).geom),
    tblb AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tbla),
    tblc AS (SELECT ST_Intersection(geom, geom2) geom FROM tblb WHERE ST_Intersects(geom, geom2)),
    tbld AS (SELECT ST_Difference(geom, geom1) geom FROM tblc WHERE ST_Intersects(geom, geom1) AND ST_Overlaps(geom, geom1))
    SELECT ST_Union(geom) geom FROM tbld
$BODY$
LANGUAGE SQL;

SELECT ST_CarvesPolygons(a.geom, b.geom) geom FROM adm_polygons a, poly_extent b 

but the result also changed, see figure,

enter image description here

What's wrong with the function syntax?

Maybe it's a bug?

Tested by: PostgreSQL 11.1, PostGIS 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1...

I have to deal with it, because there are other recursive queries waiting for me in the queue that ask me to wrap them in the SQL function :-)....

So why the results of work with WITH... and the results of the work of the corresponding SQL function in the body of which is added a query with WITH... gives different results!

0

3 Answers 3

3

Update:

To answer the actual issue with your query within a function.

First off: it's important to note here that, since you don't specify an extent_to geometry in ST_VoronoiPolygons, by default it extents the voronoi diagram to an envelope of a 50% enlarged bbox of the input geometry!

Now, in the original query from the linked post you create ST_VoronoiPolygons for a collection of all geometries in the given table, then apply your geometric operations; since a collection of 50% enlarged bboxes (see note above) covers the large polygon you want to cut out from, your geometric operations work as expected.

However, when calling your function, only a single geometry is processed by ST_VoronoiPolygons, and the enlarged bbox it uses as extent doesn't cover the large polygon to cut from. Your ST_Intersection then cuts out the small area of voronoi polygons from the large one, and from that intersection then the adm polygons. And it does that for every single adm geometry you pass into the function (since you pass them one by one)!

Hope that clears things up; if you want to make this work, you need to have the function work with a collection of input geometries. Then, however, the following approaches are better suited!

Original answer:


A few things in general:

  • technically, the query you link is not recursive
  • your query currently applies whatever logic of your function to a pair of two geometries from each table (cross joined), and returns a result for every possible pairing (as mentioned by @dr_jts)
  • a simple solution to that linked question could be

    SELECT ST_Difference(ST_Envelope(geom), geom) AS geom
    FROM   (
        SELECT ST_Collect(geom) AS geom
        FROM   <polygons>
    ) q;
    

    if the given Polygons are not overlapping, or

    SELECT ST_Difference(ST_Envelope(geom), geom) AS geom
    FROM   (
      SELECT ST_Union(geom) AS geom
      FROM   <polygons>
    ) q;
    

    for when they do overlap.

    If comparing to a generic polygon, use

    SELECT ST_Difference(b.geom, a.geom) AS geom
    FROM   (
        SELECT ST_Collect(geom) AS geom   -- ST_Union(geom) AS geom
        FROM   <polygons_used_as_blade>
    ) AS a,
    (
        SELECT ST_Envelope(ST_Collect(geom)) AS geom
        FROM   <polygon_used_to_cut_from>
    ) AS b
    ;
    

However:

If you intend to carve out all overlapping polygons for each extent, a generic and way more performant query solution for two large tables would be to collect only those polygons (blade) that actually overlap with each polygon that needs to be cut (input), e.g.

SELECT ST_Difference(a.geom, b.geom) AS geom
FROM   input AS a
CROSS JOIN LATERAL (
  SELECT ST_Collect(geom) AS geom   -- ST_Union(geom) AS geom
  FROM   blade
  WHERE  ST_Intersects(a.geom, geom)
) AS b
;

or simply (note that this may be slower due to the GROUP BY)

SELECT ST_Difference(a.geom, ST_Collect(b.geom)) AS geom   -- ST_Difference(a.geom, ST_Union(b.geom)) AS geom
FROM   input AS a
JOIN   blade AS b
  ON   ST_Intersects(a.geom, b.geom)
GROUP BY
       a.geom
;

Conclusion:

A function would only make sense if it had full access to the blade table, to collect and cut the input geometry with all relevant blade geometries; this will likely add significant overhead due to function calling and individual table lookups.

8
  • Thank you very much for the detailed explanation of the different approaches, I will take some time to reflect on the information provided, however, at this point in time I do not understand why when I run the WITH query... with the same set of geodata, I get one result, and when I turn it around with an SQL function I get another result? Is it a function design violation? or PostgreSQL does not take the reference to the source table correctly? Or something else? Feb 27, 2020 at 14:07
  • 1
    @CyrilM to be honest, I don't even know why your results are wrong; you may want to edit that info onto the question. But, if you intend to have an adm area cut out from all extents that overlaps with it, the issue is A) that you are working with many extents (compared to only one in the linked question) that each have to be checked against many admin areas. And B) you are then calling your function to take only two (!) individual geometries per call!
    – geozelot
    Feb 27, 2020 at 14:17
  • @CyrilM having a quick look at your initial query in that other post you linked, the tabled query in the WITH block joins against the full table of 'blade' polygons, but here, you only ever have a single one to check against!
    – geozelot
    Feb 27, 2020 at 14:21
  • Absolutely right, I was not satisfied with the behavior of the ST_Difference() function, because it created separate cuts for each subdomain in terms of their number, and I decided to split all areas with the ST_VoronoiPolygons function and then make paired cuts... Feb 27, 2020 at 14:23
  • 1
    Cyril, as a general rule, you should put all the relevant info in the same question. It is hard to go back and forth between questions trying to figure out what has changed. It seems like you two queries are not functionally the same, so you should not really expect the same results. Feb 27, 2020 at 19:05
2

The result you are seeing is because the query in the function is running on each single geometry that is passed in. In the original query the tbla was computed as the entire set of points from all the input geometries.

Perhaps you can call the function in a query, passing in the ST_Collect of the input table, like this:

SELECT ST_CarvesPolygons(a.geom, b.geom) geom 
 FROM (SELECT ST_Collect(geom) AS geom FROM adm_polygons) a, poly_extent b;

UPDATE: In fact this doesn't work, because the function query also needs to work against a table of the input polygons, not the collection of them. More thought required...

8
  • Edited my answer to show suggested SQL.
    – dr_jts
    Feb 28, 2020 at 16:45
  • Ah right... I overlooked that your original query is using the individual adm_polygon geometries in some places. This is a bit trickier then. Perhaps you can extract the individual elements again inside the function with an ST_Dump?
    – dr_jts
    Feb 28, 2020 at 16:59
  • I'm not following your last comment. The error makes perfect sense - it's now a matter of reworking the queries to avoid passing an entire GEOMETRYCOLLECTION to the predicates.
    – dr_jts
    Feb 28, 2020 at 17:12
  • The problem is that the query is doing some operations on individual geometries, and some on an entire table. Functions are straightforward to write for individual objecs, but more complex fo tables. In fact, you have to use dynamic SQL to work with table arguments.
    – dr_jts
    Feb 28, 2020 at 17:16
  • 1
    Well yes, calling the function with the union will work, because the union is a MultiPolygon not a GeometryCollection. But as you say, that eliminates the performance benefit.
    – dr_jts
    Feb 28, 2020 at 18:51
0

Another thing you could do is to create a function to compute the tblc result (perhaps called ST_ClippedVoronoi). Note that the neatest way to do that is to make it a set-returning function.

Then you can implement the rest of the query as a regular SQL query using the resultants output from that function. This will make a simpler query, and will avoid having to use the adm_polygons table as an argument.

0

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.