17

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 ?

obchardon
  • 1,724
  • 2
  • 13
  • 21
  • 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 09: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 09: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 09: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 09:37
  • @JohnBarça is that not the same ? if it's not shared so it "erase" the shared part ? no ? – obchardon Mar 31 '16 at 09:45
  • Without time to think thoroughly I feel that WHERE does not restrict the query as you suppose but query starts to make cartesian product. – user30184 Mar 31 '16 at 09:46
  • I think you might need to back up and understand why st_difference(t1.geom, t2.geom) from t1,t2 where st_intersects(t1.geom,t2.geom); does not give the expected results, as @user30184 has suggested. – John Powell Mar 31 '16 at 09:47
  • I understand the difference between st_instersects and st_overlaps, so as @user30184 said my query just don't do what i want. So i will edit my question. But does exist a function that is similar to the erase function in arcgis. A function that just delete the common part of a geometry ? – obchardon Mar 31 '16 at 09:56
  • Difference is right function. Run this in SQL window and you'll see select ST_AsText(ST_Difference(ST_GeomFromText('POLYGON (( 207 373, 207 719, 526 719, 526 373, 207 373 ))'),ST_GeomFromText('POLYGON (( 354 470, 354 628, 787 628, 787 470, 354 470 ))'))) Result is POLYGON (( 207 373, 207 719, 526 719, 526 628, 354 628, 354 470, 526 470, 526 373, 207 373 )) that is, part of geom1 that is not covered by geom2. Your problem is probably in what you get as geom1/geom2 pair when you process the whole layer. – user30184 Mar 31 '16 at 10:06
  • I suppose what you really want is something like "Take first feature from layer one as fid1. Select all features from layer 2 which overlap this feature. Union the features found from layer 2 and call union as eraser. Compute difference (fid1 - eraser). Take second feature from layer 1 and continue the same cycle." I am sure that this is possible to do with SQL but the query will be longer than the one that you wrote. – user30184 Mar 31 '16 at 10:19
  • @user30184 Ok great that's exactly what i need, i will try to implement something. I come back after that. – obchardon Mar 31 '16 at 10:23
  • Crude way would be to union the whole layer 2 into one big geometry and use that as an eraser for each feature of layer 1. However, such geometry may become huge. Still worth having a try. – user30184 Mar 31 '16 at 10:28
  • Something like: SELECT ST_Difference(ST_Union(t1.geom), ST_Union(t2.geom)) FROM t1, t2 though you might need to use ST_Dump to get polygons back or ST_CollectionExtract. For a small dataset, this should be fine, performance wise. I ran some tests and ST_Overlaps took about 1.5x longer than ST_Intersects, but this will depend so much on polygon complexity, index use, and any other functions being used. – John Powell Mar 31 '16 at 12:07
  • 1
    Would be nice to show EXPLAIN results of theses 2 queries.... – WKT Mar 31 '16 at 12:22
  • @azimo i edited my question instead ! – obchardon Mar 31 '16 at 13:59
  • Explain is of no use whatsoever, as ST_Overlaps and ST_Intersects both use the same condition, namely bounding box intersection, ie, &&, to determine whether to proceed with an actual intersects/overlaps check. This is why the number of rows is 13085 in both cases, but ST_Overlaps will take longer, as it has more work to do. How much more, is dependent on the complexity of the polygons and how many overlap. I would suggest you run the query that myself and user30184 have suggested, ie, the difference of two unions,see if that produces the correct results, and then we can work on optimizing it. – John Powell Mar 31 '16 at 17:28

1 Answers1

21

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));
Matt
  • 1,672
  • 16
  • 24
obchardon
  • 1,724
  • 2
  • 13
  • 21
  • 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 02 '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