Union on 3 Polygons results in a MultiPolygon with 1 Polygon

[theroggy]

Originally reported in this shapely issue: shapely/shapely#2001

Not a very problematic problem, but a bit weird: the unary_union of 3 Polygons results in a MultiPolygon with just one Polygon in it, so you would expect the result to be returned as a normal Polygon.

Script to reproduce:

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

_P1 = shapely.Polygon(
    [
        (-282.61216357247446, 148.76943812918358),
        (-281.4241833883043, 148.17032767858294),
        (-282.9997507209374, 145.04612450683473),
        (-282.20597816646364, 144.6458169603605),
        (-280.6304108338305, 147.7700201321087),
        (-279.40671534912155, 147.15289809336522),
        (-273.50213311412216, 158.86113433586473),
        (-276.70758133747506, 160.4776743716831),
        (-282.61216357247446, 148.76943812918358),
    ]
)
_P2 = shapely.Polygon(
    [
        (-279.6598724549748, 154.62355625043335),
        (-276.45442423162183, 153.00701621461496),
        (-276.4311126524318, 152.96193418063572),
        (-279.63656087578477, 154.5784742164541),
        (-279.6598724549748, 154.62355625043335),
    ]
)
_P3 = shapely.Polygon(
    [
        (-282.58885199328444, 148.72435609520434),
        (-281.4008718091143, 148.1252456446037),
        (-282.9764391417474, 145.0010424728555),
        (-282.1826665872736, 144.60073492638122),
        (-280.6070992546405, 147.72493809812946),
        (-279.38340376993153, 147.10781605938598),
        (-273.47882153493214, 158.8160523018855),
        (-276.68426975828504, 160.43259233770385),
        (-282.58885199328444, 148.72435609520434),
    ]
)

for a, b in [[_P1, _P2], [_P1, _P3], [_P2, _P3]]:
    intersection = a.intersection(b)
    union = shapely.unary_union([a, b])
    assert isinstance(intersection, shapely.Polygon)
    assert isinstance(union, shapely.Polygon)

# This becomes a MultiPolygon with one Polygon in it
union = shapely.unary_union([_P1, _P2, _P3])
print(type(union))
print(len(union.geoms))
plotter.plot_polygon(union)
plt.show()

[dr-jts]

Here's a GEOS repro:

bin/geosop -a "GEOMETRYCOLLECTION (POLYGON ((-282.61216357247446 148.76943812918358, -281.4241833883043 148.17032767858294, -282.9997507209374 145.04612450683473, -282.20597816646364 144.6458169603605, -280.6304108338305 147.7700201321087, -279.40671534912155 147.15289809336522, -273.50213311412216 158.86113433586473, -276.70758133747506 160.4776743716831, -282.61216357247446 148.76943812918358)), POLYGON ((-279.6598724549748 154.62355625043335, -276.45442423162183 153.00701621461496, -276.4311126524318 152.96193418063572, -279.63656087578477 154.5784742164541, -279.6598724549748 154.62355625043335)), POLYGON ((-282.58885199328444 148.72435609520434, -281.4008718091143 148.1252456446037, -282.9764391417474 145.0010424728555, -282.1826665872736 144.60073492638122, -280.6070992546405 147.72493809812946, -279.38340376993153 147.10781605938598, -273.47882153493214 158.8160523018855, -276.68426975828504 160.43259233770385, -282.58885199328444 148.72435609520434)))" unaryUnion

[dr-jts]

Not technically a bug, but it would be nice to reduce output to "lowest terms" in general.

[pramsey]

Erm, maybe? Or maybe "match inputs" as an underlying expectatio? (maybe? I'm not wedded to this) should maybe unary union on multipolygon to return multipolygon?

[dr-jts]

Erm, maybe? Or maybe "match inputs" as an underlying expectatio? (maybe? I'm not wedded to this) should maybe unary union on multipolygon to return multipolygon?

The usual approach in JTS/GEOS is for output to be as simple as possible. My vote is to continue with that policy, unless there's a good reason otherwise.