7

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);

where explain return:

"Nested Loop  (cost=0.15..49833.44 rows=149759 width=1796)"
"  ->  Seq Scan on affectation_com_1 t1  (cost=0.00..1447.85 rows=13085 width=923)"
"  ->  Index Scan using affectation_com_geom_idx on affectation_com t2  (cost=0.15..0.83 rows=1 width=873)"
"        Index Cond: (t1.geom && geom)"
"        Filter: _st_intersects(t1.geom, geom)"

AND

st_difference(t1.geom, t2.geom)
from t1,t2
where st_overlaps(t1.geom,t2.geom);

where explain return:

"Nested Loop  (cost=0.15..12768.08 rows=149759 width=1796)"
"  ->  Seq Scan on oca1032s_affectation_com_1 t1  (cost=0.00..1447.85 rows=13085 width=923)"
"  ->  Index Scan using oca1032s_affectation_com_geom_idx on oca1032s_affectation_com t2  (cost=0.15..0.83 rows=1 width=873)"
"        Index Cond: (t1.geom && geom)"
"        Filter: _st_overlaps(t1.geom, 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 ?

  • 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 '16 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 '16 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. – John Powell Mar 31 '16 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. – John Powell Mar 31 '16 at 9:37
  • 1
    Would be nice to show EXPLAIN results of theses 2 queries.... – WKT Mar 31 '16 at 12:22
10

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));
  • 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 '17 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 '18 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.