6

There are two multi-polygons a and A where a is supposed to be contained inside A. Using shapely, the following difference operator

b = a.difference(A)
print(b)

gives

POLYGON EMPTY

However testing A.contains(a) returns False; a.touches(A) also gives False.

How could this be possible with Shapely?

The wkt descriptions for A and a are hosted in the following Gist link (due to characters limit) https://gist.github.com/vodp/169b12f507a9730c0dd2ac86ca2589db

Edits

So what I am trying to test is assessing the equivalence between

a.difference(A) -> False

and

A.contains(a) -> True

But A.contains(a) returns False instead. If I were buffering a bit, for example

A.contains(a.buffer(-1))

then I get True. Notice that both

aa = a.buffer(-1)
aa.difference(A) -> EMPTY

and

a.difference(A) -> EMPTY

Does that mean for contains to work, none of the point of a lying in the boundary of A ?

Any chance this is because of rounding errors? I try to get the minimal buffering value and the smallest one I could get before .contains() return False is that:

A.contains(a.buffer(-0.000000001)) # 1e-9
2
  • did you check your geometries are both valid?
    – Ian Turton
    Commented Dec 6, 2020 at 12:10
  • Both are valid.
    – Pateheo
    Commented Dec 6, 2020 at 14:08

3 Answers 3

6

In your case, a.touches(A) returns False because documentation says:

touches returns True if the objects have at least one point in common and their interiors do not intersect with any part of the other.

enter image description here

Also A.contains(a) returns False because a lot of points of a lie in the extrior of A. Documentation says:

contains returns True if no points of other lie in the exterior of the object and at least one point of the interior of other lies in the interior of object.

a.difference(A) returns empty geometry because all interior of a is in A. But A.difference(a) gives:

enter image description here

1
  • Thank you. I was wrong when testing with touches indeed. What I do not understand is, if a.difference(A) returns empty geometry because all interior of a is in A then why A.contains(a) returns False? object.contains(other): Returns True if no points of other lie in the exterior of the object and at least one point of the interior of other lies in the interior of object.
    – Pateheo
    Commented Dec 6, 2020 at 13:51
5

The solution is a.symmetric_difference(A)

Reread The Shapely User Manual

object.symmetric_difference(other)
Returns a representation of the points in this object not in the other geometric object, and the points in the other not in this geometric object.

and not

object.difference(other)
Returns a representation of the points making up this geometric object that do not make up the other object.

nor object.contains(other) or object.touches(other) as Kadir Şahbaz says

A:

enter image description here

a:

enter image description here

a.symmetric_difference(A):

enter image description here

5
  • Then is it true if I say a is a subset of A thus A.contains(a) should return True ?
    – Pateheo
    Commented Dec 6, 2020 at 13:55
  • No, you must understand the binary predicates of shapely, because A.intersects(a) is true and therefore A cannot contains a
    – gene
    Commented Dec 6, 2020 at 14:13
  • By Set theory, intersects is True if contains is True. They are not exclusive aren't they? That is true for Shapely too.
    – Pateheo
    Commented Dec 6, 2020 at 14:18
  • No because intersects returns True if the boundary(exterior) or interior of the object intersect in any way with those of the other and contains returns True if no points of other lie in the exterior of the object.
    – gene
    Commented Dec 6, 2020 at 14:36
  • My argument is not contradicting yours. What I mean is: if A.contains(a) -> True then A.intersects(a) -> True as well; the inverse is not always true. And if A.difference(a) -> EMPTY and A != EMPTY and a != EMPTY then a must be contained inside A which means A.contains(a) must return True.
    – Pateheo
    Commented Dec 6, 2020 at 14:41
2

The inconsistency between a.difference(A) = POLYGON EMPTY and A.contains(a) = FALSE is due to numerical precision issues. The polygon a contains many vertices which are probably intended to lie exactly on the line segments of A. But due to the limitations of finite precision floating-point, a point almost never lies exactly on a line segment.

As visible in the image below, some of these vertices lie outside the polygon A.

Vertices not lying exactly on line segments

The contains predicate computation (in GEOS) is evaluated with high precision, and so detects that there are vertices of a which lie outside A. The reason that difference is not affected is that it contains logic to improve robustness which ends up snapping these vertices to the line.

It's hoped that a future version of JTS/GEOS will provide spatial predicates that can use a tolerance value to handle situations like this.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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