title: Shapely intersection with multiple MultiPolygons where a percentage of Multipoloygons intersect?
body: <p>I'm using Shapely's polygons for human-generated data. Multiple humans were asked to draw polygons around certain features in images. For each image, we thus have n <code>MultiPolygon</code>'s, where n equals the number of participants per image.</p>
  <p>I can plot each of these <code>Multipolygon</code>'s.</p>
  <pre><code>fig, ax = plt.subplots()
  
  for ii, multi_poly in enumerate(multi_polys):
      for poly in multi_poly.geoms:
          x,y = poly.exterior.xy
          plt.plot(x,y, c=colors[ii])
  </code></pre>
  <p><a href="https://i.stack.imgur.com/lYmsM.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/lYmsM.png" alt="enter image description here" /></a></p>
  <p>We can see that at some locations, the <code>Multipolygon</code>'s overlap, while at others there is no overlap.</p>
  <p>I wish to get the overlap, or intersections, of these polygons.
  This should be trivial, as I can do something like:</p>
  <pre><code>intersection = multi_a.intersection(multi_b) \
              .intersection(multi_c) \
              .intersection(multi_d) \
              .inters... 
  </code></pre>
  <p>I can plot this intersection on the previous plot and we see:</p>
  <p><a href="https://i.stack.imgur.com/JTMgZ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/JTMgZ.png" alt="enter image description here" /></a></p>
  <p>This looks pretty good. However, this method only returns the regions where all <code>Multipoloygon</code>'s overlap.  Is there a way to get the intersection where 75% of the polygons overlap? Or where 50% overlap?</p>
  <p>A code example: The following dummy data gives this figure:</p>
  <pre><code>P1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
  P2 = Polygon([(2.5, 2), (3, 2), (3, 3), (2.5, 3)])
  multi_a = MultiPolygon([P1, P2])
  
  P1 = Polygon([(-1, -1), (-1, 2),(2, 2), (2, -1)])
  P2 = Polygon([(3, 3), (4, 2), (4, 4), (3, 4)])    
  multi_b = MultiPolygon([P1,P2])  
  
  P1 = Polygon([(-2, 4), (2.2, 4),(2.2, -2), (-2, -2)])
  P2 = Polygon([(-1.5, 3), (-1.1, 3), (-1.5, -1), (-1.1, -1)])    
  multi_c = MultiPolygon([P1,P2])
  
  P1 = Polygon([(2.5, -1), (3.2, -1),(3.2, -2), (-2, -2)])
  P2 = Polygon([(3, 0), (4, 0), (3, 1), (4, 1)])
  multi_d = MultiPolygon([P1,P2])
  </code></pre>
  <p><a href="https://i.stack.imgur.com/1ATYX.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/1ATYX.png" alt="enter image description here" /></a></p>
  <p>On these four multipolygons, the intersection method would returns no intersection as there is not a single spot that all four multipolygons occupy. However, the blue square, accentuated with yellow marker, is occupied by the blue, orange and green polygon. Thus 75% of the multipolygons overlap at this location.</p>
  <p>Is there a way (preferably using Shapely) to get the locations where 75% of polygons overlap?</p>
  <hr />
  <p>The accepted answer appeared to break down in certain conditions.
  I found the problem: the boundaries of the contain shape sometimes overlap the inter polygons. I can replace <code>if geom.contains(polygon):</code> with <code>if geom.contains(polygon.buffer(-1)):</code> and it works exactly as intended.</p>
  <p><a href="https://i.stack.imgur.com/SbkjA.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/SbkjA.png" alt="enter image description here" /></a></p>