Using ST_Difference to remove overlapping features?

Question

I am trying to use ST_Difference to create a set of polygons (processing.trimmedparcelsnew) that do not contain any of the area covered by another set of polygons (test.single_geometry_1) using PostGis 2.1 (and Postgres SQL 9.3). Here is my query:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference
FROM 
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig;

But the resulting polygons have not been trimmed, instead they seem to have been split where they intersect with the other layer. I have tried just running the select without putting the result in a table and everything else I can think of, but I can't seem to get this function to work.

I have attached a picture of the result

---

After comments, I have tried adding a WHERE clause. I want the parcels that have no intersections, and the intersecting areas of the other parcels removed (the layer test.single_geometry represents contamination I want removed from my parcels). I tried an intersect but of course I actually want the non intersections so I am now trying a disjoint. I have also tried adding the orig to my table but the documentation for ST_Difference (http://postgis.net/docs/ST_Difference.html) does say it returns the exact geometry I need (a geometry that represents that part of geometry A that does not intersect with geometry B), so I am confused as to why I would want the original polygon in my table instead. Anyway, here is my modified code:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference, orig.geom AS geom
FROM 
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig
WHERE ST_Disjoint(orig.geom, cont.geom);

---

Following from dbaston's answer I have now tried:

CREATE TABLE processing.parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

The result of this is just a copy of test.multi_geometry_1. Though now the splitting is no longer occurring.

I tried the earlier version, but again just get a copy of test.multi_geometry_1:

CREATE TABLE processing.parcels_trimmed_no_coalesce AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

I am starting to wonder if there is something else I am doing wrong? The proceeding statement is:

DROP TABLE IF EXISTS processing.parcels_trimmed_no_coalesce;

And I am running the queries from the PostgreSQL SQL query window and Openjump.

The statement I use to see the table is:

SELECT * FROM processing.parcels_trimmed_no_coalesce;

---

In the interest of simplification I have now reduced this query down to just:

SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.geometriestocutagainst b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.geometriestocut a;

This still results in just the original polygons (test.geometriestocut) when the desired result is the original trimmed against test.geometriestocutagainst.

Answer 1

A self-join allows you to operate on the relationship between pairs of two features. But I don't think you're interested in pairs: for each feature, you want to operate on the relationship between that feature and all other features in your dataset. You can accomplish this with a subquery expression:

CREATE TABLE parcels_trimmed AS
SELECT id, ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                FROM parcels b
                                WHERE ST_Intersects(a.geom, b.geom)
                                  AND a.id != b.id))
FROM parcels a;

You might see something odd in the results, though. Parcels that don't have any overlaps are being dropped entirely! That's because the ST_Union aggregate on an empty recordset is going to be NULL, and ST_Difference(geom, NULL) is NULL. To smooth this over, you need to wrap your ST_Difference call in a COALESCE:

CREATE TABLE parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM parcels b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM parcels a;

This means that if the result of ST_Difference is NULL, the coalesced expression will evaluate to the original geometry.

The above query will remove overlapping areas from your domain entirely. If you instead want to pick a winner, you could do a.id < b.id, or some other criterion, instead of a.id != b.id.

Answer 2

I had the same problem as you. I don't know if you already found a solution to your problem, but I modified the accepted answer above and I got what I wanted.

CREATE TABLE parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Collect(b.geom) 
                                         FROM parcels b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         )), a.geom)
FROM parcels a;

Answer 3

I use ST_DifferenceAgg() from the PostGIS Addons. You have to merge the two tables together, have a unique identifier and an index on the geometry column. Here is a short example:

WITH overlappingtable AS (
  SELECT 1 id, ST_GeomFromText('POLYGON((0 1, 3 2, 3 0, 0 1), (1.5 1.333, 2 1.333, 2 0.666, 1.5 0.666, 1.5 1.333))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((1 1, 3.8 2, 4 0, 1 1))')
  UNION ALL
  SELECT 3 id, ST_GeomFromText('POLYGON((2 1, 4.6 2, 5 0, 2 1))')
  UNION ALL
  SELECT 4 id, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))')
  UNION ALL
  SELECT 5 id, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))')
)
SELECT a.id, ST_DifferenceAgg(a.geom, b.geom) geom
FROM overlappingtable a,
     overlappingtable b
WHERE a.id = b.id OR -- Make sure to pass at least once the polygon with itself
      ((ST_Contains(a.geom, b.geom) OR -- Select all the containing, contained and overlapping polygons
        ST_Contains(b.geom, a.geom) OR
        ST_Overlaps(a.geom, b.geom)) AND
       (ST_Area(a.geom) < ST_Area(b.geom) OR -- Make sure bigger polygons are removed from smaller ones
        (ST_Area(a.geom) = ST_Area(b.geom) AND -- If areas are equal, arbitrarily remove one from the other but in a determined order so it's not done twice.
         a.id < b.id)))
GROUP BY a.id
HAVING ST_Area(ST_DifferenceAgg(a.geom, b.geom)) > 0 AND NOT ST_IsEmpty(ST_DifferenceAgg(a.geom, b.geom));

This will merge the overlapping parts with the biggest overlapping polygon. If you want to keep the overlapping part separated look at the ST_splitAgg() example.