# Looping through table to get single intersection from N>2 geometries using PostGIS?

It's really clear how to apply ST_Intersection when you have only two geometries. But I am struggling with the situation where you have N>2 shapes and you want to compute a single one that represents intersection of all of them as a single geometry.

With the help of ST_Intersection: Intersection of all geometries in a table I can get intersections for all pairs of geometries in a single table. But it is not exactly what I am trying to achieve though.

Let me try to illustrate the problem. What I have is this (simplified to only three geometries): But what I am after is this: I have over 5000 shapes in the table, so obviously I can't specify them explicitly.

I imagine a loop could be a solution, something like:

``````result = union(geometries)
for geom in geometries:
result = intersection(result, geom)
``````

How would you achieve this with PostgreSQL/PostGIS?

• do all of your geometries in one table and in one field? and do you want to get just one geometry (the result of the intersection of all of them)?
– Moh
Feb 19 '18 at 21:31
• Records are in one table, one polygon per field (there are multiple polygons actually, but they don't overlap). I want only one resulting geometry, yes. Feb 19 '18 at 21:41
• so all of your polygons in one (column)?. and you would like as your drawing shows to get just one polygon (intersection result of the of all of the intersected polygons). Sorry I repeat my question just to understand totally
– Moh
Feb 19 '18 at 21:45
• yes, that's exactly right. Feb 19 '18 at 21:47

This is a great application for a user-defined aggregate function, and I'm a bit wondering why this particular aggregate doesn't already exist in PostGIS. At its core, an aggregate function needs to do nothing more than iterate over a set of rows, maintaining a state (of type `stype`), and repeatedly calling a function (`sfunc`) that transforms a state and a row into a new state.

If this sounds like the `reduce` operator in languages like JavaScript, you're right - it's the exact same thing.

Aggregates can be complex to define (see docs), and it doesn't help that there are multiple versions of the syntax. But this is a pretty simple case:

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

Here's an example of this aggregate in action, producing the inside of a Venn diagram:

``````CREATE table test (geom geometry);
INSERT INTO test SELECT ST_Buffer(ST_MakePoint(0,0), 0.75);
INSERT INTO test SELECT ST_Buffer(ST_MakePoint(1,0), 0.75);
INSERT INTO test SELECT ST_Buffer(ST_MakePoint(0.5,1), 0.75);

SELECT full_intersection(geom) FROM test;
``````
• It doesn't exist because this would be very inefficient for a large number of geometries, except in the case that few intersect and you rapidly end up with the empty geometry. Feb 20 '18 at 3:22
• @jpmc26, what would be a good practice for the case you describe? Initial reduction of the set? Feb 20 '18 at 11:17
• @mikitk I was working on an alternative answer, but I wasn't able to finish it yet. I'll post it later. Feb 20 '18 at 11:37
• Good explanation of aggregates. I don't think it would be hugely inefficient, if you first fed it the dumped result of ST_Union. Feb 20 '18 at 14:26
• @dbaston is there no constraint to prohibit stype becoming empty/NULL during iteration? a fallback check would be nice...one could possibly implement a custom function/wrapper for ST_Intersection, but that seems somewhat clumsy...+1 anyways Feb 21 '18 at 0:09

You can do the following using a loop by converting your geom to an array:

``````create or replace function ST_IntersectionArray(geom geometry[]) returns geometry as \$\$
declare
i integer;
tmpGeom geometry;
begin
tmpGeom := geom;
FOR i IN 1..array_length(geom,1) LOOP
tmpGeom:= ST_Intersection(tmpGeom,geom[i]);
END LOOP;
return tmpGeom;
end;
\$\$
LANGUAGE plpgsql;
``````

And the result (table=your table name)

``````SELECT ST_IntersectionArray(ARRAY(select geom from table));
``````

to check! insert your result in a new table (result)

``````Create table result as
SELECT ST_IntersectionArray(ARRAY(select geom from table));

select * from table result
``````

Tested (the result in QGIS) working.

Reference (ST_Intersection of multiple polygons in one table)

PS: I tried the answer using aggregate, it gives me the same with my answer, actually perfect to use aggregate. but sometimes using arrays has an advantage if someone wants to work on a just specific number of polygons

• is there anyway you can explain the code from the begin to end statement. Feb 20 '18 at 15:11
• tmpGeom := geom; why geom? shouldnt it just be an empty array? Feb 20 '18 at 15:11
• FOR i IN 1..array_length(geom,1) LOOP 1..array_length(geom,1) LOOP don't understand this syntax Feb 20 '18 at 15:12
• defining geom to access just one geom each time in the loop iteration. The second line is to loop over all geom and add one geom as row in the array
– Moh
Feb 20 '18 at 19:12
• @ziggy that's assigning the first array element to tmpGeom, you just need one to start with for `ST_Intersection` - in the loop, assign the result of `ST_Intersection` between the initial geom and the next to tmpGeom - repeat. you could start with 2 in the `FOR` loop here btw. Feb 20 '18 at 23:52

The aggregate function answer is perhaps the most elegant way to do this. But another option is to use a recursive CTE. This may provide more flexibility in some situations. It's also a nice use of the standard recursive CTE "running total" pattern. Here's a solution, with some example data:

``````WITH RECURSIVE
data(geom) AS (VALUES
('POLYGON ((50 150, 100 200, 250 50, 200 0, 50 150))'),
('POLYGON ((200 200, 250 150, 100 0, 50 50, 200 200))'),
('POLYGON ((120 200, 180 200, 180 0, 120 0, 120 200))'),
('POLYGON ((50 130, 50 70, 250 70, 250 130, 50 130))')
),
seq AS
(
SELECT geom::geometry, ROW_NUMBER() OVER () AS i
FROM data
),
intgeoms AS (
SELECT i, geom, geom AS intgeom
FROM seq
WHERE i = 1
UNION ALL
SELECT seq.i, seq.geom, ST_Intersection(intgeom, seq.geom) AS intgeom
FROM intgeoms INNER JOIN seq
ON seq.i = intgeoms.i + 1
)
SELECT i, intgeom
FROM intgeoms
ORDER BY i DESC LIMIT 1;
``````

The way it works is:

• the `seq` query computes a sequential index for each geometry in the set to be intersected
• the `intgeoms` query interates over each geometry in the sequence, keeping a "running total" of the intersections