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
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
In your case, a.touches(A)
returns False
because documentation says:
touches
returnsTrue
if the objects have at least one point in common and their interiors do not intersect with any part of the other.
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:
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:
a:
a.symmetric_difference(A):
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
.
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
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP