6

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 '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. – mikitk 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. – mikitk Feb 19 '18 at 21:47
8

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 '18 at 3:22
  • 1
    @jpmc26, what would be a good practice for the case you describe? Initial reduction of the set? – mikitk Feb 20 '18 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 '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. – John Powell Feb 20 '18 at 14:26
  • 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 – ThingumaBob Feb 21 '18 at 0:09
4

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

  • is there anyway you can explain the code from the begin to end statement. – ziggy Feb 20 '18 at 15:11
  • tmpGeom := geom[1]; why geom[1]? shouldnt it just be an empty array? – ziggy 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 – ziggy Feb 20 '18 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 '18 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. – ThingumaBob Feb 20 '18 at 23:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.