13

I'm working with Shapely polygons and I need a way to delete all smaller polygons contained within a bigger polygon. I tried using the .contains() method which Shapely provides, but the method doesn't return True if a given smaller polygon isn't completely inside of the bigger, "parent" polygon.

Basically, I want a method like .contains() but that it returns True in case the inner polygon shares boundaries with the outer polygon like on the example image below.

here.

Here are the polygons from the picture presented in wkt format:

The green one:

POLYGON Z ((14.4265764858233823 45.3396418051734784 0.0000000000000000, 14.4267228266679606 45.3395430970275015 0.0000000000000000, 14.4266753563381904 45.3394727193694536 0.0000000000000000, 14.4265290154936121 45.3395714275154376 0.0000000000000000, 14.4265764858233823 45.3396418051734784 0.0000000000000000))`

The red one:

POLYGON Z ((14.4265450394689161 45.3395951840357725 0.0000000000000000, 14.4265695507109317 45.3395786509942837 0.0000000000000000, 14.4265802185605700 45.3395944667317679 0.0000000000000000, 14.4265982245953417 45.3395823215079616 0.0000000000000000, 14.4265715327703994 45.3395427492501426 0.0000000000000000, 14.4265290154936121 45.3395714275154376 0.0000000000000000, 14.4265450394689161 45.3395951840357725 0.0000000000000000))

I also tried using the .intersects() method but it returns True for polygons outside of a given polygon which have some common boundaries, which I do not want.

I hope you understand what I need and I'm thankful if someone knows a solution to this.

kfilipcic
  • 130
  • 1
  • 7
  • `contains` should normally work. Create, for example, two simple square polygons one inside another that share a side to see that it works. In your case, there is probably a problem with precision. Could you provide a definition of both polygons, so that we could reproduce the issue? – Georgy Jul 06 '19 at 14:52
  • I added the wkt strings of polygons to the question. I assume that is enough info about the polygons? It could be a problem with precision, but I need this particular precision. Also, thank you for the reply and edits from you and @eyllanesc. This is my first question here so I'm not quite sure what I am doing. – kfilipcic Jul 06 '19 at 15:53
  • The example polygons in wkt format that you have provided are actually trimmed and don't show all the digits (you probably just did `polygon.wkt`, right?). Because if I load them and use `green.contains(red)` then it will show `True`. To get all the digits, could you please use [`wkt.dumps`](https://shapely.readthedocs.io/en/stable/manual.html#shapely.wkt.dumps) with `trim=False` argument and add them in the question? So, it would be: `from shapely import wkt; wkt.dumps(polygon, trim=False)` – Georgy Jul 06 '19 at 15:58
  • I added the output from wkt.dumps. Yes, I just used polygon.wkt. (I'm not so experienced with using shapely). Thanks for the information – kfilipcic Jul 06 '19 at 16:04

2 Answers2

12

Normally the contains method should work when testing if one polygon is inside another and they have common borders. For example, if you would take the following simple example, it would work as expected:

from shapely.geometry import Polygon

a = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
b = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
a.contains(b)
# True

enter image description here

But quite often what happens is that due to precision errors the inner polygon comes out just a tiny bit from the outer and the test fails.

Here, for example, I plotted your polygons and zoomed in on the upper left intersection point:

import matplotlib.pyplot as plt

plt.plot(*green.exterior.xy, c='g')
plt.plot(*red.exterior.xy, c='r')

You can see that lines don't lie on each other perfectly:

enter image description here

There are several ways to deal with this problem. The first one, for example, was proposed in How to deal with rounding errors in Shapely and in some Shapely issues on GitHub:

  1. Shrink the smaller polygon or extend the bigger one a bit:

    big.contains(small.buffer(-1e-14))
# True
big.buffer(1e-14).contains(small)
# True

  2. Check that area of the smaller polygon that lies outside of the bigger one is close to zero:

    small.difference(big).area < 1e-14
# True

  3. Check if distances of each vertex of a smaller polygon to the bigger polygon are close to zero:

    from shapely.geometry import Point

vertices = map(Point, small.exterior.coords)
distances = map(big.distance, vertices)
all(distance < 1e-14 for distance in distances) 
# True

Probably there are more ways to perform the test, but I think these will suffice.

Georgy
  • 12,464
  • 7
  • 65
  • 73
  • The first solution you provided does the trick. Thanks! The second solution does work, but I'm not sure if it would sometimes return `True` for polygons close to the bigger polygon. As for the third solution, it doesn't work in my case because it returns `True` for outside polygons close to the bigger one, which I do not want in my case. – kfilipcic Jul 06 '19 at 18:06
  • @kfilipcic review the [tour] – eyllanesc Jul 06 '19 at 20:52
  • this small value is also helpful when simplifying polygons .. thank you – roberto tomás Jan 04 '22 at 17:31
1

Not properly an answer, but just a minimal sample in which the contains also fails:

from shapely.geometry import Polygon, Point

a = Polygon([
    (786.50468384, 805.6206089),
    (1500.32201405, 393.44262295),
    (1918.1206089, 513.34894614),
    (1920., 1080.),
    (792.03651158, 1077.95349473),
    (786.50468384, 882.43559719)
])

b = Point((1704, 1080))

print(a.contains(b))  # False
print(a.intersects(b))  # False
print(a.touches(b))  # False

enter image description here

Taking a look at this answer, it seems the proper predicate to use in such cases is intersect. However, it also fails for the example above.

Even with buffer, I had to use a eps of 1 to get a True response.

print(a.buffer(1).contains(b))
print(a.buffer(1).intersects(b))
print(a.buffer(1).touches(b))
pceccon
  • 9,379
  • 26
  • 82
  • 158