14

I'm trying to use ST_difference and ST_intersection with two layers of multipolygons.

I want something like:

enter image description here

With the input on the left and the result on the right.

MY DATA:

t1.geom contains 1000 differents polygons.

t2.geom contains 1000 differents polygons.

THE QUERIES:

I try to create (in a naive way) two queries to perform these operations

st_intersection(t1.geom, t2.geom)
from t1,t2
where st_intersects(t1.geom,t2.geom);

Unfortunately it doesn't work, or takes forever to run.

I would like a query similar to the tools erase and intersect from ArcGIS.

QUESTION:

How to perform those operations efficiently ?

16
  • You should change only one thing at a time. How long does it take to run st_difference(t1.geom, t2.geom) from t1,t2 where st_intersects(t1.geom,t2.geom);
    – user30184
    Mar 31, 2016 at 9:17
  • @user30184 I have already try this option, but the result doesn't make sense, that don't give me t1.geom - t2.geom, in fact i don't understand the result, i just want to erase t2.geom that overlaps t1.geom
    – obchardon
    Mar 31, 2016 at 9:29
  • I agree with what user30184 says, but would add that ST_Overlaps is a more complete test than ST_Intersects, as it is false in the case of containment. ST_Intersects will return true as soon as any two line segment from geom1 and geom2 intersect, whereas ST_Overlaps may have to test for containment, which depending on the complexity of the polygons, may take longer. Mar 31, 2016 at 9:31
  • st_difference(t1.geom, t2.geom) will give you the part of t1.geom that is not shared with t2.geom. It will not erase t2.geom, if I have understood you correctly. Mar 31, 2016 at 9:37
  • 1
    Would be nice to show EXPLAIN results of theses 2 queries....
    – WKT
    Mar 31, 2016 at 12:22

1 Answer 1

19

I finally found an answer by myself.

For st_difference():

with temp as 
(
  select   b.gid, st_union(a.geom) as geom
  from     t1 b join t2 a on st_intersects(a.geom, b.geom)
  group by b.gid
) 
select st_difference(b.geom,coalesce(t.geom, 'GEOMETRYCOLLECTION EMPTY'::geometry)) as newgeom
from t1 b left join temp t on b.gid = t.gid

It works like a charm and runs fast. The query is similar for st_intersection().

In addition, some of my polygons weren't valid so I corrected the geometry with

update t1
set geom = ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3));
update t2
set geom = ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3));
2
  • in the coalesce function what is the purpose of 'GEOMETRYCOLLECTION EMPTY'::geometry? I am assuming you have it there in case the t.geom is null, but how could you turn a 'GEOMETRYCOLLECTION EMPTY'::geometry into a real geometry?
    – ziggy
    Aug 2, 2017 at 13:36
  • It appears that ST_Difference(geom, 'GEOMETRYCOLLECTION EMPTY'::geometry) is just the object itself, i..e, geom. I think the purpose of the COALESCE statement is therefore to make sure that the query returns a result even in cases when a and b don't intersect.
    – John
    Sep 22, 2018 at 21:12

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.