title: Substracting inner rings from a list of polygons
body: <p>I have a large list of polygons (>10^6) most of which are non-intersecting but some of these polygons are hole of another polygon (~10^3 cases). Here a image to explain the problem, the smaller polygon is a hole in the larger polygon but both are independent polygon in the list of polygons. <a href="https://i.stack.imgur.com/VrExH.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/VrExH.png" alt="A large polygon with smaller polygon on the inside."></a></p>
  
  <p>Now I would like to efficiently determine which polygons are holes and substract the holes i.e. subtract the smaller polygons which lie completely inside in another polygon and return a list of "cleaned" polygons. A pair of hole and parent polygon should be transformed like this (so basically hole subtracted from the parent):
  <a href="https://i.stack.imgur.com/zMjZu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/zMjZu.png" alt="enter image description here"></a></p>
  
  <p>There are plenty of similar questions on Stackoverflow and gis.stackexchange.com but I haven't found one that actually solves this problems. Here are some related questions: 
  1. <a href="https://gis.stackexchange.com/questions/5405/using-shapely-translating-between-polygons-and-multipolygons">https://gis.stackexchange.com/questions/5405/using-shapely-translating-between-polygons-and-multipolygons</a>
  2. <a href="https://gis.stackexchange.com/questions/319546/converting-list-of-polygons-to-multipolygon-using-shapely">https://gis.stackexchange.com/questions/319546/converting-list-of-polygons-to-multipolygon-using-shapely</a></p>
  
  <p>Here is a sample code. </p>
  
  <pre><code>from shapely.geometry import Point
  from shapely.geometry import MultiPolygon
  from shapely.ops import unary_union
  import numpy as np
  
  #Generate a list of polygons, where some are holes in others; 
  def generateRandomPolygons(polygonCount = 100, areaDimension = 1000, holeProbability = 0.5):
      pl = []
      radiusLarge = 2 #In the real dataset the size of polygons can vary
      radiusSmall = 1 #Size of holes can also vary
  
      for i in range(polygonCount):
          x, y = np.random.randint(0,areaDimension,(2))
          rn1 = np.random.random(1)
          pl.append(Point(x, y).buffer(radiusLarge))
          if rn1 &lt; holeProbability: #With a holeProbability add a hole in the large polygon that was just added to the list
              pl.append(Point(x, y).buffer(radiusSmall))
      return pl
  
  polygons = generateRandomPolygons()
  print(len(pl))
  </code></pre>
  
  <p>Output looks like this: <a href="https://i.stack.imgur.com/2J8r8.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/2J8r8.png" alt="enter image description here"></a></p>
  
  <p>Now how can I create a new list of polygons with the holes removed. Shapely provides functions to subtract one polygon from another (<a href="https://shapely.readthedocs.io/en/stable/manual.html#spatial-analysis-methods" rel="nofollow noreferrer">difference</a>) but is there a similar function for lists of polygons (maybe something like <a href="https://shapely.readthedocs.io/en/latest/manual.html#shapely.ops.unary_union" rel="nofollow noreferrer">unary_union</a> but where overlaps are removed)? Alternative how to efficiently determine which are holes and then subtract them from the larger polygons? </p>