Intersection/Addition of multiple tables in PostGIS

Question

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.

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.