Turn a MultiPolygon with Nested Shells into a Polygon with Holes

Question

I am reading GeoJSON into shapely objects and need help,

I have a massive MultiPolygon with several nested Polygons in it. I believe this is why it is not considered a valid shape?

as_shape.is_valid
-> False

explain_validity(as_shape)
-> 'Nested shells[-118.752368927 36.63641068]'

Are the inner polygons actually the problem?

Right now we attempt to fix all invalid poly/multipolygons we get using the first buffer ring, and throw an error if that does not make the shape valid

        first_buffer = as_shape.buffer(0)
        if not first_buffer.is_valid:
            raise Exception('Invalid geo json contains invalid geometries')

        return mapping(first_buffer)

In Many cases this works, but in the case of a multipolygon with several complex holes, the MultiPolygon gets oversimplified (holes go away and we get the outer shell).

What's the right way to make all inner polygons holes while keeping the outer MultiPolygon in tact?

My attempts so far seem overly complicated --

        all_ps = list(as_shape)
        get = {
            'shells': [ g for g in all_ps if g and not len([ g.buffer(0).within(i) for i in all_ps if i!=g ]) > 2],
            'holes': [ g for g in all_ps if g and len([ g.buffer(0).within(i) for i in all_ps if i!=g ]) > 2],
        }
    

Answer 1

You could clip out the nested shells using outer.difference(inner)

from shapely.geometry import Point, MultiPolygon
from shapely.validation import explain_validity


multi = MultiPolygon([Point(0, 0).buffer(5), Point(0, 0).buffer(2)])

multi = MultiPolygon(sorted(multi, reverse=True, key=lambda p: p.area))  # ensure largest is 1st

outer = multi[0]  # 1st part of the multipoly
inner = multi[1:]  # handle overlapping inner polys

for poly in inner:
    outer = outer.difference(poly)

print(explain_validity(multi)) 
print(explain_validity(outer))    

>>> Nested shells[2 0]
>>> Valid Geometry

multi

outer

Answer 2

I am going to answer with what worked for us, but it is also built on a key assumption about your underlying data, so I am not going to accept my own answer since YMMV. If there is some wildly clever way to solve this universally, perhaps using something I don't know about geojson, please share!

---

What we ended up determining is that we only got this Nested Shells issue on MultiPolygons that had 3 levels of Polygons on top of each other. What we determined was that the "holes" were being represented by the polygons that were "one level" down (covered only by the outer most geometries). So we ended up with something like this

    all_ps: MultiPolygon = shp.geoms

    solid = MultiPolygon(
        [
            poly for poly in all_ps 
            if len([p for p in all_ps if p.covers(poly)]) % 2 != 0
        ]
    )
    with open('fire.json', 'w') as fp:
        fp.write(json.dumps(mapping(fire)))

    holes = MultiPolygon(
        [
            poly for poly in all_ps if len([p for p in all_ps if p.covers(poly)]) % 2 == 0 
        ]
    )

    return MultiPolygon([
        Polygon(
            shell=solid.exterior.coords,
            holes=[nf.exterior.coords for nf in holes if fire.covers(nf)]
        ) for fire in holes
    ])

This works for us. It MIGHT have been good to distinguish exterior solids to better capture actual Polygon holes, but we found that this works fine.