I want to intersect multiple table together. For each intersection geometry, I want to sum the value of all the layers intersecting at that place.

My current query does it by step, getting the intersection between two tables first and getting the new value by summing for each new geometry, and then finding the intersection of these two tables with a third table, etc. Eventually I can have up to 10 tables to combine and sum the value.

My current query to do so for 3 tables is as follow:

WITH intersect_ab AS (
SELECT ST_INTERSECTION(a.geometry, b.geometry) as geometry, a.val + b.val as val
FROM mytable_a a
INNER JOIN mytable_b b ON ST_INTERSECTS(a.geometry, b.geometry)
     intersect_ab_c AS (
SELECT ST_INTERSECTION(ab.geometry, c.geometry) as geometry, ab.val + c.val as val
FROM intersect_ab ab
INNER JOIN mytable_c c ON ST_INTERSECTS(ab.geometry, c.geometry)
SELECT geometry, val
FROM intersect_ab_c;

Here is the query plan associated with this query:

              QUERY PLAN
 CTE Scan on intersect_ab_c  (cost=2359.00..2371.08 rows=604 width=40) (actual time=0.735..707.197 rows=2201 loops=1)
   CTE intersect_ab
     ->  Nested Loop  (cost=0.14..179.04 rows=4857 width=40) (actual time=0.611..156.574 rows=1251 loops=1)
           ->  Seq Scan on mytable_b b  (cost=0.00..1.32 rows=32 width=211) (actual time=0.003..0.010 rows=32 loops=1)
           ->  Index Scan using mytable_a_gix on mytable_a a  (cost=0.14..4.79 rows=1 width=217) (actual time=0.115..0.400 rows=39 loops=32)
                 Index Cond: (geometry && b.geometry)
                 Filter: _st_intersects(geometry, b.geometry)
                 Rows Removed by Filter: 1
   CTE intersect_ab_c
     ->  Nested Loop  (cost=0.14..2179.96 rows=604 width=40) (actual time=0.734..705.353 rows=2201 loops=1)
           ->  CTE Scan on intersect_ab ab  (cost=0.00..97.14 rows=4857 width=40) (actual time=0.612..157.457 rows=1251 loops=1)
           ->  Index Scan using mytable_c_gix on mytable_c c  (cost=0.14..0.42 rows=1 width=135) (actual time=0.022..0.032 rows=2 loops=1251)
                 Index Cond: (ab.geometry && geometry)
                 Filter: _st_intersects(ab.geometry, geometry)
                 Rows Removed by Filter: 0
 Planning Time: 12.183 ms
 Execution Time: 707.623 ms

I have a geometry index on all tables (mytable_a, mytable_b and mytable_c). mytable_a has 1027 polygons, mytable_b has 32 polygons and mytable_c has 373 polygons.

Is there a way to make this more efficient? It seems like there ought to be for a "simple"/"likely common" spatial join like this one, but I'm not sure what the best approach would be.

1 Answer 1


At a high level, it's most efficient to compute a coverage of all the polygons. Then you can query each resultant polygon in the coverage against all the parent tables to sum up the values of the parents for that polygon.

To compute a polygonal coverage, extract the boundaries of the polygon, union all the lines to node them, and then polygonize the result. This is outlined in Polygon Averaging in PostGIS at the Crunchy Data blog.

  • That's a very interesting idea, which should work indeed. I decided to forgo that route for now due to the large size of my vector tables. For reference, I created a series of rasters, which I then align together in GDAL and use numpy to make quick calculations, before reuploading to PostGIS. It's not elegant, but it was faster than what I was trying to do and gave me what I wanted. Commented Apr 24, 2020 at 5:46

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.