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):

enter image description here

But what I am after is this:

enter image description here

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, 2018 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.
    – mikitk
    Feb 19, 2018 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, 2018 at 21:45
  • yes, that's exactly right.
    – mikitk
    Feb 19, 2018 at 21:47

3 Answers 3


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.
    – jpmc26
    Feb 20, 2018 at 3:22
  • 1
    @jpmc26, what would be a good practice for the case you describe? Initial reduction of the set?
    – mikitk
    Feb 20, 2018 at 11:17
  • 1
    @mikitk I was working on an alternative answer, but I wasn't able to finish it yet. I'll post it later.
    – jpmc26
    Feb 20, 2018 at 11:37
  • 1
    @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
    – geozelot
    Feb 21, 2018 at 0:09
  • 1
    For newer versions of postgis, I had to wrap the builtin st_intersection into a separate function because of the new gridSize parameter. So I replaced the sfunc parameter with the wrapped function. Sep 21, 2021 at 17:12

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.

Example data with total intersection

Here's a solution, with some example data:

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

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

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 $$
   i integer;
   tmpGeom geometry;
    tmpGeom := geom[1];
    FOR i IN 1..array_length(geom,1) LOOP
      tmpGeom:= ST_Intersection(tmpGeom,geom[i]);
    return tmpGeom;
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.
    – ziggy
    Feb 20, 2018 at 15:11
  • tmpGeom := geom[1]; why geom[1]? shouldnt it just be an empty array?
    – ziggy
    Feb 20, 2018 at 15:11
  • FOR i IN 1..array_length(geom,1) LOOP 1..array_length(geom,1) LOOP don't understand this syntax
    – ziggy
    Feb 20, 2018 at 15:12
  • 1
    defining geom[1] 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, 2018 at 19:12
  • 1
    @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.
    – geozelot
    Feb 20, 2018 at 23:52

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.