Wrong result of unary_union: islands in polygons being unioned get "lost"

[theroggy]

When union-ing two polygon that are both ~donut shaped, the result doesn't have a hole in it anymore.

Originally reported here: https://github.com/shapely/shapely/issues/1633

Visual (left: green + blue: overlapping donut-like shapes / right: result of union):

Script to reproduce:

import shapely
import shapely.wkt
from shapely import plotting
import matplotlib.pyplot as plt

# TODO shapely bug
poly1 = shapely.wkt.loads(
    "POLYGON ((-70.1162076595306 -1.534802639761859, -69.11618088109869 0.1972327069961337, -58.117783208841075 -6.152921874253107, -82.91750004113058 -49.107291593330054, -92.53033522436974 -36.35721051634815, -71.84825846863409 -0.5348026424198639, -70.1162076595306 -1.534802639761859), (-70.1162076595306 -1.534802639761859, -90.1362204660169 -36.210482111666536, -83.13912274760473 -45.49115337178261, -60.849818553267646 -6.88494590501217, -70.1162076595306 -1.534802639761859))")
poly2 = shapely.wkt.loads(
    "POLYGON ((-54.71170009216383 -7.253412859063702, -80.4528216636855 -51.8383434227556, -82.83439152620278 -50.463343426410354, -85.03022855188276 -52.1188795198414, -95.62396528443203 -38.067769761534805, -93.42812825875205 -36.41223366810376, -95.80969812126934 -35.037233671758514, -73.50445135763194 3.59658713514124, -71.12288149511465 2.2215871387959965, -69.76459062209818 4.612726177890158, -69.7335565179279 4.595097198664126, -71.09184739094437 2.2039581595699644, -69.71681057060047 4.585506761362204, -55.718233134337225 -3.4968642609262157, -57.09326995468111 -5.878412862718456, -54.71170009216383 -7.253412859063702), (-60.84981855326765 -6.884945905012167, -70.1162076595306 -1.5348026397618595, -90.1362204660169 -36.210482111666536, -83.13912274760473 -45.4911533717826, -60.84981855326765 -6.884945905012167))")
print(f"{poly1.is_valid=}")
print(f"{poly2.is_valid=}")

u = shapely.unary_union([poly1, poly2])

if u.area <= (poly1.area + poly2.area) + (eps := 1):
    print(f"Unary union OK: area ({u.area}) <= {(poly1.area + poly2.area)=}")
else:
    print(f"Unary union not OK: area ({u.area}) > {(poly1.area + poly2.area)=}")

union = poly1.union(poly2)
if union.area <= (poly1.area + poly2.area) + (eps := 1):
    print(f"Union OK: area ({union.area}) <= {(poly1.area + poly2.area)=}")
else:
    print(f"Union not OK: area ({union.area}) > {(poly1.area + poly2.area)=}")

fig, ax = plt.subplots(ncols=2, sharey=True)

plotting.plot_polygon(ax=ax[0], polygon=poly1, color="blue", alpha=0.5)
plotting.plot_polygon(ax=ax[0], polygon=poly2, color="green", alpha=0.5)
plotting.plot_polygon(ax=ax[1], polygon=u, color="red", alpha=0.5)
plt.show()

[pramsey]

This one replicates in JTS OverlayNGRobust, @dr-jts