4

I have a dataset with two layers of polygons:

Layer A: Areas of interest from the corine land cover inventory dataset

Layer B: Road network resulting from buffered linestrings that were aggregated by st_union (in combination with st_dump to expand the multipolygon)

Below are layer A and B (excerpt from dataset): Layer A and B

I want to subtract layer B from layer A by using st_difference. My SQL query looks like this:

CREATE TABLE diff AS
SELECT st_difference(clc.geom, roads.geom) AS geom
FROM clc, roads
WHERE st_intersects(clc.geom, roads.geom);

When looking at the result I was very confused at first: Difference of A - B

There was no difference in the upper part! While troubleshooting I changed the layer opacity and I noticed, that the difference was below a polygon with the same shape but without the difference: Difference of A - B with transparency

Somehow the output of st_difference results in multiple overlapping polygons and some of them are the result of the difference and some are not.

Layer A and B do not have multiple overlapping polygons themselves.

1
6

You have done a full cross table join, so you have intersected every road against every land cover polygon in turn, so you have (N-1) polygons that don't intersect for ever polygon that does (where N is the number of roads).

To get the answer you expect you need to union the land cover polygons and then subtract the roads from that polygon (or union the roads and subtract that from each polygon in turn). So something like:

SELECT st_difference(clc.geom, roads.geom) AS geom
FROM clc, (SELECT st_union(geom) from roads) as roads;
1
  • You likely want to add an index driven (spatial) filter condition into the correlated sub-query.
    – geozelot
    Oct 24 at 11:43

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.