9

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

enter image description here


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.

  • You did not specify a WHERE clause, so you may have polynomial expansion in the resulting table. How many rows are in trimmedparcelsnew ? – Vince Jul 23 '15 at 16:10
  • If you only want the difference where they intersect, you could try adding WHERE ST_Intersects(orig.geom, cont.geom). Otherwise, the difference of two polygons that do not intersect, is the original polygon. – John Powell Jul 23 '15 at 16:11
  • There are 24 rows in trimmed parcel new, I want the difference even when they do not intersect, so would I correct that I need to use orig.geom in the table rather than difference? – Mart Jul 24 '15 at 8:18
  • A disjoint test should produce polynomial expansion -- each feature in cont to appear once for every feature in orig that doesn't overlap, and the difference will never change an input geometry – Vince Jul 24 '15 at 12:50
  • Ok thanks for the clarification, but I am still not sure then why the original code is not working. If ST_Difference(orig.geom, cont.geom) returns the geometries in a that do not intersect with b, than why does the table contain the split geometry's rather than the geometries in a that do not intersect b. – Mart Jul 24 '15 at 13:09
14

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.

  • Thanks for the response, unfortunately I am having trouble getting this to work for me, instead just end up with the original polygon (a). I will edit my question with more info. Thanks again. – Mart Jul 27 '15 at 13:35
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;
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.