Using ST_Difference to remove overlapping features?
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
.
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;
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.