TransWikia.com

Multipolygon supposed to be contained by another but returned False with Shapely

Geographic Information Systems Asked by Pateheo on December 16, 2020

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

3 Answers

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

Answered by Kadir Şahbaz on December 16, 2020

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

Answered by gene on December 16, 2020

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.

Answered by dr_jts on December 16, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP