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

Question

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?

Answer 1

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;

Answer 2

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

Answer 3

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[1];
    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