Using ST_Difference

Question

I have 2 tables Table A with this fields (379 rows in this table)

The geom of table A contains this (QGIS only used for better visualization)

The table B contains this fields (97 rows in this table)

The Field cobertura5 contains a buffer of 20 km

I made the difference operation on QGIS first and the result is this (298 rows on QGIS)

In PostGIS, I've tried this sentence looking for the same result

Create Table prueba AS
Select ST_Difference(tableA.geom,tableB.cobertura5) From tableA, tableB

And get this on QGIS (35162 rows) (Exactly the same as tableA)

I get stuck there with no ideas of what to do maybe a WHERE or maybe include another operation on the sentences?

The multiple tables in a FROM clause style query fell out of favor nigh on 30 years ago, when JOIN was added to the SQL standard. One of the main issues with the FROM list was forgetting to add a WHERE clause element to enforce a relationship between the tables, resulting in polynomial expansion. Since not all features overlap all other features, you're going to have nearly a hundred copies of every row in A in the result set. I suggest finding a good SQL tutorial, where you can work on the JOIN clause before integrating spatial types and functions into the mix. — Vince, Apr 17 '21 at 00:11

score 4 · Accepted Answer · answered Apr 19 '21 at 20:49

4

This recent post has a good solution to this problem. For your data the code looks like:

SELECT ST_Multi(COALESCE(
         ST_Difference(a.geom, blade.geom),
         a.geom
       )) AS geom
FROM   tableA AS a
CROSS JOIN LATERAL (
  SELECT ST_Union(b.cobertura5) AS geom
  FROM   tableB AS b
  WHERE  ST_Intersects(a.geom, b.cobertura5) 
) AS blade;

answered Apr 19 '21 at 20:49

dr_jts

5,038
18
15

Thanks, that works really well and give me the expecting results – Jeff_Barahona Apr 20 '21 at 03:41
This will work as far as PostGIS is able to union the complete coverage into a single polygon. That's not the case of you have hundreds of thousands of polygons. – Pierre Racine Apr 20 '21 at 14:17
The query only unions the polygons which intersect each polygon from the A table. That should work fine for the dataset shown. – dr_jts Apr 20 '21 at 18:21
@dr_jts With this dataset, no doubts. We need ST_DifferenceAgg() in geos! – Pierre Racine Apr 21 '21 at 13:28
1

@PierreRacine handling a massive amount of data in reasonable memory is going to require a new approach to overlay. Hopefully that can get worked on sometime. – dr_jts Apr 22 '21 at 18:19
@dr_jts The simple ST_DifferenceAgg() works very fine on millions of polygons. I'm just saying it would be nice to have a C++ version of it. – Pierre Racine Apr 22 '21 at 18:37
How is the performance? It seems like it must run millions of Difference operations? – dr_jts Apr 22 '21 at 19:19
For each tableA polygon, it only ST_Difference() the overlapping polygons from tableB. So about 2 or 3 tableB polygons per tableA polygons. Pretty fast and robust actually. Il would be faster if the aggregate state functions would be in C++. Other than simply doing the difference between two layers like here, we use it to remove overlaps on the about 40 millions forest stand polygons covering the Canadian boreal forest. No memory overflow like ST_Union() would do. Have a look at the code! – Pierre Racine Apr 23 '21 at 14:10
But if there are only a few intersecting B polys per A poly, then isn't the above query plenty fast too? And it uses the spatial index. – dr_jts Apr 23 '21 at 20:03
Because of the LATERAL JOIN hum? I guess you are right. Wonder if that would work to remove overlaps. – Pierre Racine Apr 23 '21 at 20:20

score 2 · Answer 2 · answered Apr 19 '21 at 15:19

In PostGIS each polygon from table B is removed from each polygon from table A individually and the resulting difference is not cumulative. i.e. polygon B1 and B2 of table B are removed from two different instances of polygon A1 of table A producing two polygons instead of one.

The trick is to accumulate all the differences from table B for each polygon of table A. Cumulate here means aggregate. Unfortunately there is still no ST_Difference() aggregate function in PostGIS. But one has been developed in the PostGIS Addons, a set of advanced plpgsql functions for PostGIS..

So download the Addons and execute them in your database. This will enable those functions.

Then you can try this query:

SELECT ST_DifferenceAgg(tableA.geom, tableB.geom) geom
FROM tableA, tableB
WHERE ST_Intersects(tableA.geom, tableB.geom)
GROUP BY tableA.geom;

Works well, it was easy with the Addons. Thanks – Jeff_Barahona Apr 20 '21 at 03:43 — Jeff_Barahona, Apr 20 '21 at 03:43

Using ST_Difference

2 Answers2

Linked