# Turn a MultiPolygon with Nested Shells into a Polygon with Holes

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],
}

``````

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
``````

# `outer`

• Nice! This is a smart approach. unfortunately the data does not reliably feature the "outer" as the first item. Wish it did, seems like it should be a part of the geojson standard. Commented Oct 12, 2021 at 20:14
• Maybe sort by area? Largest == outer? `multi = MultiPolygon(sorted(multi, reverse=True, key=lambda p: p.area))` Commented Oct 12, 2021 at 20:53
• Thats a perfect solution if we can assume there is only one "shell" - in our case we technically can't do that :( Commented Oct 13, 2021 at 17:49

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.