How to use the shapely.geometry.MultiPolygon function in shapely

def randomMultiPolygon():
    return shapely.geometry.MultiPolygon([randomPolygon() for i in xrange(R.randint(1, 8))])

def load_limited_to(limited_to):
    require_geom_support()
    srs = SRS(limited_to['srs'])
    geom = limited_to['geometry']

    if not hasattr(geom, 'type'): # not a Shapely geometry
        if isinstance(geom, (list, tuple)):
            geom = bbox_polygon(geom)
        else:
            polygons = load_polygon_lines(geom.split('\n'))
            if len(polygons) == 1:
                geom = polygons[0]
            else:
                geom = shapely.geometry.MultiPolygon(polygons)

    return GeomCoverage(geom, srs, clip=True)

# this function uses a default of 4 pixel connectivity for grouping pixels into features
    shapes = rasterio.features.shapes(mask.astype(np.int16), mask > 0)

    polygons = []
    for shape, val in shapes:
        s = shapely.geometry.shape(shape).exterior

        if use_buffer:
            s = shapely.geometry.polygon.Polygon(s.buffer(buffer_size))
        else:
            s = shapely.geometry.polygon.Polygon(s)

        if s.area > min_polygon_area:
            polygons.append(s)

    mp = shapely.geometry.MultiPolygon(polygons)

    if isinstance(mp, shapely.geometry.Polygon):
        df = pd.DataFrame({
            'area_size': [mp.area],
            'poly': [mp],
            'image_id': [image_id]
        })
    else:
        df = pd.DataFrame({
            'area_size': [p.area for p in mp],
            'poly': [p for p in mp],
            'image_id': [image_id] * len(mp)
        })

    df = df.sort_values(by='area_size', ascending=False)
    df.loc[:, 'wkt'] = df.poly.apply(lambda x: shapely.wkt.dumps(x, rounding_precision=0))

from __future__ import division

import sys
import warnings

import param
import numpy as np
import shapely.geometry as sgeom

from cartopy import crs as ccrs
from shapely.geometry import (MultiLineString, LineString, MultiPolygon,
                              Polygon, LinearRing, Point, MultiPoint)
from holoviews.core.util import basestring

geom_types = (MultiLineString, LineString, MultiPolygon, Polygon,
              LinearRing, Point, MultiPoint)
line_types = (MultiLineString, LineString)
poly_types = (MultiPolygon, Polygon, LinearRing)


def wrap_lons(lons, base, period):
    """
    Wrap longitude values into the range between base and base+period.
    """
    lons = lons.astype(np.float64)
    return ((lons - base + period * 2) % period) + base


def project_extents(extents, src_proj, dest_proj, tol=1e-6):
    x1, y1, x2, y2 = extents

def _combine_to_multipolygon(df, column, *names):

    all_poly = [df[df[column] == name].geometry.values[0] for name in names]

    combined_poly = geometry.MultiPolygon(all_poly)

    df.loc[df[column] == names[0], "geometry"] = gp.GeoSeries(combined_poly).values

    for name in names[1:]:
        df = df.loc[df[column] != name]

    return df

scaled_gdf = gpd.sjoin(scaled_gdf,finnish_regions_gdf,how='inner',op='intersects')

    ### Switch the polygon geometry back to the 'geometry' field and drop uselesss columns
    scaled_gdf['geometry'] = scaled_gdf['polygon_geometry']
    scaled_gdf.drop(['index_right','polygon_geometry'],axis=1, inplace=True)

    ### Encode the region name with the One-hot encoding (= pandas dummy encoding) 
    ### method so machine learning can understand it better
    encoded_gdf = pd.get_dummies(scaled_gdf['NAMEFIN'])

    ### Join scaled gdf and encoded gdf together
    scaled_and_encoded_gdf = scaled_gdf.join(encoded_gdf).drop('NAMEFIN',axis=1)

    ### The resulting dataframe has Polygon and Multipolygon geometries. 
    ### This upcasts the polygons to multipolygon format so all of them have the same
    scaled_and_encoded_gdf["geometry"] = [MultiPolygon([feature]) if type(feature) == Polygon else feature for feature in scaled_and_encoded_gdf["geometry"]]
    print("Dataframe size after adding region name: " + str(len(scaled_and_encoded_gdf.index))+ " zip codes with " + str(len(scaled_and_encoded_gdf.columns)) + " columns")

    return scaled_and_encoded_gdf

def mask2polygons_layer(mask, epsilon=1.0, min_area=10.0):
    # first, find contours with cv2: it's much faster than shapely
    contours, hierarchy = cv2.findContours(((mask == 1) * 255).astype(np.uint8), cv2.RETR_CCOMP, cv2.CHAIN_APPROX_TC89_KCOS)

    # create approximate contours to have reasonable submission size
    if epsilon != 0:
        approx_contours = simplify_contours(contours, epsilon)
    else:
        approx_contours = contours

    if not approx_contours:
        return MultiPolygon()

    all_polygons = find_child_parent(hierarchy, approx_contours, min_area)

    # approximating polygons might have created invalid ones, fix them
    all_polygons = MultiPolygon(all_polygons)

    all_polygons = fix_invalid_polygons(all_polygons)

    return all_polygons

def keys(self):
        return self._storage.keys()

    def pop(self, *args, **kwargs):
        return self._storage.pop(*args, **kwargs)

    def update(self, dictionary):
        self._storage.update(dictionary)

    def values(self):
        return self._storage.values()


class SpatialCollection(AbstractCollection):
    _multi_cast = {'Point': MultiPoint, 'Polygon': MultiPolygon}

    def __init__(self, meta=None, key=None, crs=None, headers=None, value_keys=None):
        super(SpatialCollection, self).__init__()

        self._crs = None

        self.meta = meta
        self.key = key
        self.crs = crs
        self.headers = headers
        self.value_keys = value_keys

        self.ugeom = OrderedDict()

    @property
    def crs(self):

def get_box_split_bounds_geometry(lon_min: float, lat_min: float,
                                  lon_max: float, lat_max: float) -> shapely.geometry.base.BaseGeometry:
    box_1, box_2 = get_box_split_bounds(lon_min, lat_min, lon_max, lat_max)
    if box_2 is not None:
        return shapely.geometry.MultiPolygon(polygons=[shapely.geometry.box(*box_1), shapely.geometry.box(*box_2)])
    else:
        return shapely.geometry.box(*box_1)

"""

    #open shapefile and read shapes
    Shapes = fiona.open(ShapeFile)

    # get the input coordinate system
    Input_CRS = Proj(Shapes.crs)

    # Create a dictionary of shapely polygons
    PolygonDict = {}

    # loop through shapes and add to dictionary
    for Feature in Shapes:
        if Feature['geometry']['type'] == 'MultiPolygon':
            Shape = MultiPolygon(shape(Feature['geometry']))
            Value = float(Feature['properties']['ID'])
        elif Feature['geometry']['type'] == 'Polygon':
            Shape = Polygon(shape(Feature['geometry']))
            Value = float(Feature['properties']['ID'])

        #check for multipolygons
        if Value in PolygonDict:
            Polygons = [Shape, PolygonDict[Value]]
            Shape = cascaded_union(Polygons)
        #update dictionary
        PolygonDict[Value] = Shape

    return PolygonDict, Input_CRS