Aggregate version of ST_Intersection

Question

I would like to use an aggregate version of the ST_Intersection function in PostGIS. In a table I have multiple polygons and would like to calculate the intersection of all of them. Ideally it would even be possible to group by one of the other columns. In short, I would like to use ST_Intersection as an aggregate function, or something similar to

SELECT groupid, ST_Intersection(geometry) FROM tablename GROUP BY groupid;

The result of this would be a table with a line for each groupid. The resulting geometry for each groupid should be the intersection of all polygons with this groupid.

Is there a way to do this in PostGIS?

Example: My table contains four records, each containing a polygon. They are shown in the left column in the picture below. In the picture assume that the squares are every time the same square (e.g. [-1,1]x[-1,1]). Each of the records has a groupid, shown in text and with a color. In the output I want a record for each groupid. The polygon in the output should be the intersection of all input polygons from that group. In the picture they are shown in the right column. In this pictures that means the record for groupid = 1 (green) should contain the intersection of all green input polygons and analogous for groupid = 2 (red). Here, the intersection of both green triangles is a smaller triangle and the intersection of the red rectangles is a smaller red square.

The example in code: The code for this example would look like this.

 CREATE TEMPORARY TABLE simpleexample
(
    id INT, 
    groupid INT,
    polygon GEOMETRY
);

INSERT INTO simpleexample VALUES
(1, 1, ST_GeomFromText('POLYGON((-1 1, -1 -1, 1 -1, -1 1))')), -- lowerleft half of square [-1,1]x[-1,1]
(2, 1, ST_GeomFromText('POLYGON((1 1, -1 -1, 1 -1, 1 1))')), -- lowerright half of square [-1,1]x[-1,1]
(3, 2, ST_GeomFromText('POLYGON((-1 0, 1 0, 1 -1, -1 -1, -1 0))')), -- lower half of square [-1,1]x[-1,1]
(4, 2, ST_GeomFromText('POLYGON((0 1, 0 -1, 1 -1, 1 1, 0 1))')); -- left half of square [-1,1]x[-1,1]

-- making the intersections manually works fine
SELECT ST_AsText(ST_Intersection(ST_GeomFromText('POLYGON((-1 1, -1 -1, 1 -1, -1 1))'), ST_GeomFromText('POLYGON((1 1, -1 -1, 1 -1, 1 1))')));
-- correctly returns "POLYGON((1 -1,-1 -1,0 0,1 -1))"
SELECT ST_AsText(ST_Intersection(ST_GeomFromText('POLYGON((-1 0, 1 0, 1 -1, -1 -1, -1 0))'), ST_GeomFromText('POLYGON((0 1, 0 -1, 1 -1, 1 1, 0 1))')));
-- correctly returns "POLYGON((0 0,1 0,1 -1,0 -1,0 0))"

-- doing is all in one go doesn't work, this is what I would like to do
SELECT groupid, ST_Intersection(polygon) FROM simpleexample GROUP BY groupid
-- gives error "function st_intersection(geometry) does not exist"
-- I want this to be the two previously calculated lines in one table.

Obviously, I would have to do this for many more polygons than is handy for copy-pasting.

its a little unclear of what you want to do here. if you are looking to AGGREGATE geometries than just ST_UNION is the geometric aggregate. think of it like a SUM or COUNT function but for geometries — ziggy, Jan 31 '18 at 13:25
do you want to check if two geomrties are intersecting? in this case use ST_intersects(geomA, geomB) which returns a bool value and can used to aggregate geomtries see: https://postgis.net/docs/ST_Intersects.html — LaughU, Jan 31 '18 at 13:36
@ziggy I see where the confusion comes from: I am not looking for a function that aggregates geometries, but for an aggregated function in the database sense (like your examples SUM and COUNT are, but also MEAN etc). Just like a SUM takes the sum of all values, I am looking for a way to calculate the intersection of all polygons. — Some wandering yeti, Jan 31 '18 at 13:55
@LaughU I don't want to know if they intersect, but the intersection of the two polygons itself. So the result should be a geometry, not boolean. That is what ST_Intersection does. But I don't want to do it for two polygons but for a few thousand. In fact, I have about 1000 groups of each 1000 polygons. For each group I want the intersection of all 1000 polygons in that group. And that is where copy-pasting is no longer an option. — Some wandering yeti, Jan 31 '18 at 13:59
so you want to intersect each row by the others?? this seems like a case for a self join — ziggy, Jan 31 '18 at 15:14
adding a picture of what your original simpleexample looks like would help a lot — ziggy, Jan 31 '18 at 15:22
I believe you need to do a recursive query, similar in spirit to this answer. You would have to modify it somewhat, but the basic idea is the same . — John Powell, Jan 31 '18 at 16:20

ziggy · Answer 1 · 2018-01-31T16:49:14.207

I am going to attempt an answer based on your comments:

First

drop table if exists simpleexample;
 CREATE table simpleexample
(
    id INT, 
    groupid INT,
    polygon GEOMETRY
);

INSERT INTO simpleexample VALUES
(1, 1, ST_GeomFromText('POLYGON((-1 1, -1 -1, 1 -1, -1 1))')), -- lowerleft half of square [-1,1]x[-1,1]
(2, 1, ST_GeomFromText('POLYGON((1 1, -1 -1, 1 -1, 1 1))')), -- lowerright half of square [-1,1]x[-1,1]
(3, 2, ST_GeomFromText('POLYGON((-1 0, 1 0, 1 -1, -1 -1, -1 0))')), -- lower half of square [-1,1]x[-1,1]
(4, 2, ST_GeomFromText('POLYGON((0 1, 0 -1, 1 -1, 1 1, 0 1))')); -- left half of square [-1,1]x[-1,1]

select * from simpleexample;

running this yields

NOW if you want to run an intersection on JUST the polygons that have the same groupid you need to run the table against itself like this

select a.id a_id,a.groupid a_groupid,b.id b_id,b.groupid b_groupid,st_intersection(a.polygon,b.polygon) geometry
    from simpleexample a cross join simpleexample b
    where a.groupid=b.groupid and a.id<>b.id

the where conditions are if the groupid's equal eachother and the ids do not equal eachother to unsure you are not intersecting the features on themselves. this query will duplicate the features because it is a cross join that retains all the ids

this is just to illustrate what the cross join does.

Finally if you just solely want the geometric intersection of the features within each group just pull out the id columns and only return the intersection and groupid like this

select st_union(geometry) geom,groupid 
    from(select a.groupid,st_intersection(a.polygon,b.polygon) geometry
            from simpleexample a cross join simpleexample b
            where a.groupid=b.groupid and a.id<>b.id
            ) t
group by groupid

UPDATE based on your comments

try this query, essentially it will assess if there are any duplicate geometric intersections (which there will be for each id intersection) and return only 1 record for that geometry and aggregate which id combination formed that geometry.

select distinct on(geometry) geometry, array_agg(distinct geometry) geoms,groupid,array_agg('['|| aid || ',' || bid || ']')
        from(select a.id aid,b.id bid, a.groupid,st_intersection(a.polygon,b.polygon) geometry
                from simpleexample a cross join simpleexample b
                where a.groupid=b.groupid and a.id<>b.id
                ) t
    group by geometry,groupid

Thank you very much. This does indeed what I want for the simple example (added a picture to my question as well). But if I understand it correctly the st_union at the end will return the union of all intersections between two features with the same groupid, while I would like to keep making intersections. Example: I add a fifth record with the same feature as the fourth. The intersection of id 3 with 4 or 5 would give the square that is now returned as intersection between 3 and 4. But the intersection between 4 and 5 would return 4 (or 5, it's the same). And making the union of those... — Some wandering yeti, Jan 31 '18 at 16:03
... would give feature 4 (or 5, still the same), while I would like the result to be that small square we had as intersection of 3 and 4. — Some wandering yeti, Jan 31 '18 at 16:04
Basically, I would like to replace the st_union in your solution by st_intersection. But unfortunately st_intersection can't be used as an aggregate function the way st_union can. — Some wandering yeti, Jan 31 '18 at 16:07
okay now I understand better, you may need some sort of recursive option here. my first select query might be the route you are trying to go but it will give duplicates for each id intersection pair...like 1 intersects 2 will be the same as 2 intersects 1...is keeping their original id's important? — ziggy, Jan 31 '18 at 16:22

dr_jts · Answer 2 · 2021-03-16T18:10:44.230

2

This answer presents a user-defined aggregate function to perform the intersection of a set of geometries. Here's the (slightly modified) code and an example:

CREATE AGGREGATE ST_IntersectionAgg (
  basetype = geometry,
  stype = geometry,
  sfunc = ST_Intersection
);

WITH data(geom) AS (VALUES
( ST_Buffer(ST_Point(0,0), 0.75) ),
( ST_Buffer(ST_Point(1,0), 0.75) ),
( ST_Buffer(ST_Point(0.5,1), 0.75) )
)
SELECT ST_IntersectionAgg(geom) FROM data;

edited Mar 16 '21 at 18:10

answered Mar 16 '21 at 17:05

dr_jts

5,038
18
15

1

Link-only answers, even links to GIS SE, are subject to removal if they don't provide an actual answer to the question within the Answer body. What you have here is more appropriate as a comment than as an Answer. – Vince Mar 16 '21 at 17:23

score 1 · Accepted Answer · answered Feb 09 '18 at 11:18

Based on the input of ziggy and John Powell, I managed to write a recursive query that does what I wanted to obtain.

--  add an index number that counts the polygons within each groupid
WITH RECURSIVE indexed AS (SELECT ROW_NUMBER() OVER(PARTITION BY groupid) idx, groupid, polygon FROM simpleexample),

-- for each groupid, find the highest value of idx
lastindex AS (SELECT groupid, MAX(idx) maxidx FROM indexed GROUP BY groupid),

-- recursively make all the intersections 
t AS (  
    SELECT idx, groupid, polygon FROM indexed WHERE idx = 1     
    UNION ALL   
    SELECT i.idx, i.groupid, ST_Intersection(i.polygon, t.polygon) polygon 
    FROM indexed i JOIN t ON i.idx = t.idx + 1 AND i.groupid = t.groupid 
)

SELECT
  t.groupid,
  polygon
FROM  t 
JOIN lastindex -- join to filter out the highest idx
  ON t.groupid = lastindex.groupid AND t.idx = lastindex.maxidx;

Aggregate version of ST_Intersection

3 Answers3