# Why the result of executing a geometric user function is different from the result of execution of the corresponding CTE query

So,

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

``````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,

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,

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!

## 3 Answers

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.

• 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
• @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 `extent`s 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! 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! 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
• 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

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...

• Edited my answer to show suggested SQL. 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`? 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. 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. Feb 28, 2020 at 17:16
• 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. Feb 28, 2020 at 18:51

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.