Calculando el polígono final después de restar todos los polígonos intersecantes utilizando GeoJson – Python.

Tengo un Polígono llamado “Jumeirah Islands Clusters” como se muestra en la figura. Quiero restar los archivos GeoJson del lago 1 y del lago 2 del Polígono de las Islas. Intenté la siguiente solución, pero no está funcionando para crear un Polígono final.

Lo intenté dos veces para obtener el resultado deseado, pero la salida mostrada no es deseada ya que el lago 2 no está restado en el paso final.

from shapely.geometry import Polygon, mapping

with open("file_polygon1.geojson") as f:
    feature = json.load(f)

if feature["geometry"]["type"] == "Polygon":
    polygon_1 = Polygon([(coor[0], coor[1]) for coor in feature["geometry"]["coordinates"][0]])

with open("file_polygon2.geojson") as f:
    feature_2 = json.load(f)

if feature_2["geometry"]["type"] == "Polygon":
    polygon_2 = Polygon([(coor[0], coor[1]) for coor in feature_2["geometry"]["coordinates"][0]])

new_geometry = mapping(polygon_2.difference(polygon_1))

new_feature = dict(type="Feature", id="", properties=dict(Name=""), geometry=dict(type=new_geometry["type"], coordinates=new_geometry["coordinates"]))

outjson = dict(type="FeatureCollection", features=[new_feature])