Methods: polygonize and bounding_polygon are incomplete for non-contiguous 2d topologies

[Huite]

Polygonize collects all possible polygons:

        collection = shapely.polygonize(shapely.linestrings(coordinates[edges]))
        polygon = max(collection.geoms, key=lambda x: _bbox_area(x.bounds))

This will work as desired with holes and such, but not if there e.g. two non-contiguous part of the mesh.

The correct check here would be to find the largest bounding boxes that are not contained by any other. Arguably the return object should be a MultiPolygon.