Shapely intersection with multiple MultiPolygons where a percentage of Multipoloygons intersect?

One approach is to split all geometries to get a flat list of non-intersecting regions on the XY-plane, and then see how many of the original geometries contain each region. Any region that is contained by at least some threshold number of the original geometries can be added to the result. This is easier to explain using a combination of code and illustrations.

To start with, we'll need to fix one issue. The example you put together has a couple invalid geometries, which will cause Shapely to throw an error when attempting to query spatial relationships (e.g., calling contains or intersects). You can check that using the is_valid property, and get a more detailed information by calling explain_validity:

from shapely.geometry import Polygon
from shapely.validation import explain_validity

P2 = Polygon([(-1.5, 3), (-1.1, 3), (-1.5, -1), (-1.1, -1)])

>>> P2.is_valid
False

>>> explain_validity(P2)
'Self-intersection[-1.3 1]'

Basically, it's not happy about shapes like these being expressed as a Polygon when they should be multi-polygons:

So to make your example valid, some of your multi-polygons will have 3 rather than 2 constituent polygons:

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.3, 1)])
P3 = Polygon([(-1.5, -1), (-1.3, 1), (-1.1, -1)])
multi_c = MultiPolygon([P1,P2,P3])

P1 = Polygon([(2.5, -1), (3.2, -1),(3.2, -2), (-2, -2)])
P2 = Polygon([(3, 0), (4, 0), (3.5, 0.5)])
P3 = Polygon([(3.5, 0.5), (3, 1), (4, 1)])
multi_d = MultiPolygon([P1,P2,P3])

Hopefully your real source data has valid geometries (or you have some way of converting them to become valid geometries - which incidentally is an upcoming feature in Shapely 1.8, implemented via make_valid, but it is not yet released), since otherwise the method described below will not work.

With that out of the way, the first step is to get a flat list of non-intersecting regions from your list of shapes. To do this, we start with the original list of intersecting shapes (note darker shading where multiple shapes overlap):

Convert them to lines using linemerge (in combination with unaryunion):

Then polygonize the result:

It might not be clear from the picture, but the idea is none of these geometries intersect - some of these polygons have holes in them (in cases where one shape previously contained another). So this represents the "flat list of non-intersecting regions on the XY-plane" that I was referring to in the beginning.

The code for the process so far looks like this:

from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import linemerge, unary_union, polygonize

# Original shape definitions here (snipped)...

shapes = [multi_a, multi_b, multi_c, multi_d]

lines = unary_union(linemerge([geom.exterior for shape in shapes for geom in shape.geoms]))
polygons = list(polygonize(lines))

Next, we check each resulting region in the polygons list, and check how many shapes from the original list it intersected with. If it's above the threshold (defined here as 0.75 * len(shapes), then we add it to the result:

threshold = 0.75 * len(shapes)

def overlaps(polygon, shape):
    for geom in shape.geoms:
        if geom.contains(polygon):
            return True
    return False

result = []

for polygon in polygons:
    containing_shapes = []
    for shape in shapes:
        if overlaps(polygon, shape):
            containing_shapes.append(shape)
    if len(containing_shapes) >= threshold:
        result.append(polygon)

If you're dealing with a large dataset, checking intersections in a nested loop like that can be quite slow (O(N^2)), so you can probably speed it up using an STRtree:

from shapely.strtree import STRtree

# (Previous code here to get the flattened list of polygons...)

tree = STRtree([geom for shape in shapes for geom in shape.geoms])
result = []

for polygon in polygons:
    matches = [o.wkt for o in tree.query(polygon) if o.contains(polygon)]
    if len(matches) >= threshold:
        result.append(polygon)