Python shapely.geometry.MultiPolygon() Examples

def AttemptCollapse(name, zone):
  overlaps = False
  for i in range(0,len(zone)):
    for j in range(i+1, len(zone)):
      if zone[i].overlaps(zone[j]):
        overlaps = True
        break
    if overlaps:
      break

  if overlaps:
    print 'Found overlapping zone on %s!' % name
    mp = SMultiPolygon(zone)
    collapsed = cascaded_union(mp)
    if type(collapsed) == SPolygon:
      return [collapsed]
    else:
      print 'Got multipolygon len %d' % len(collapsed.geoms)
      return collapsed.geoms
  else:
    print 'Non-overlapping zone on %s (%d)' % (name, len(zone))
    return zone 

def surround_with_holes(geometry, hole_spacing, hole_radius, padding, max_distance):
    """
    Surrounds the given geometry with holes, which are arranged in a square lattice around the structure.
    This can be used for generating vortex traps like presented in https://doi.org/10.1103/PhysRevApplied.11.064053

    :param geometry: The geometry around which the holes are generated
    :param hole_spacing: Spacing between the holes
    :param hole_radius: Radius of the holes
    :param padding: Padding around the geometry
    :param max_distance: Maximum distance of a hole from the geometry
    :return: Shapely object, which describes the holes
    """
    geometry = geometric_union(geometry if isinstance(geometry, (tuple, list)) else (geometry,))
    buffer_around_waveguide = geometry.buffer(max_distance)
    area_for_holes = prep(buffer_around_waveguide.difference(geometry.buffer(hole_radius + padding)))
    area = buffer_around_waveguide.bounds
    points = (Point(x, y) for x in np.arange(area[0], area[2], hole_spacing) for y in
              np.arange(area[1], area[3], hole_spacing))
    return MultiPolygon([point.buffer(hole_radius) for point in points if area_for_holes.contains(point)]) 

def _force_geometry_2d(geometry):
    """ Convert a geometry to 2d
    """
    if geometry["type"] in ["Polygon", "MultiLineString"]:
        geometry["coordinates"] = [
            _force_linestring_2d(l) for l in geometry["coordinates"]
        ]
    elif geometry["type"] in ["LineString", "MultiPoint"]:
        geometry["coordinates"] = _force_linestring_2d(geometry["coordinates"])
    elif geometry["type"] == "Point":
        geometry["coordinates"] = geometry["coordinates"][:2]
    elif geometry["type"] == "MultiPolygon":
        geometry["coordinates"] = [
            [_force_linestring_2d(l) for l in g] for g in geometry["coordinates"]
        ]

    return geometry 

def _line_iterator(obj):
    if isinstance(obj, (sg.LineString)):
        yield obj
    elif isinstance(obj, (sg.MultiLineString)):
        for obj2 in obj.geoms:
            yield obj2
    elif isinstance(obj, (sg.Polygon)):
        yield sg.LineString(obj.exterior)
        for obj2 in obj.interiors:
            yield sg.LineString(obj2)
    elif isinstance(obj, (sg.MultiPolygon)):
        for obj2 in obj.geoms:
            yield sg.LineString(obj2.exterior)
            for obj3 in obj2.interiors:
                yield sg.LineString(obj3)
    else:
        try:
            tup = tuple(obj)
        except TypeError:
            raise TypeError('Could not use type %s' % type(obj))
        else:
            for obj2 in tup:
                for line in _line_iterator(obj2):
                    yield line 

def PolyWithoutSmallHoles(poly, min_hole_area_km2=0.5):
  """Returns input |shapely| geometry with small internal holes removed.

  Args:
    poly: A shapely polygon geometry (Polygon or MultiPolygon) defined in WGS84
       or in NAD83 (lon, lat) coordinates (degrees).
  min_hole_area_km2: the minimum area for holes to be kept (in km2).

  Returns:
    (a |shapely| geometry) The input geometry with all small holes removed.
  """
  if isinstance(poly, sgeo.Polygon):
    interiors = [interior
                 for interior in poly.interiors
                 if GeometryArea(sgeo.Polygon(interior)) > min_hole_area_km2]
    if len(interiors) == len(poly.interiors):
      return poly
    else:
      return sgeo.Polygon(poly.exterior, interiors)
  elif isinstance(poly, sgeo.MultiPolygon):
    return sgeo.MultiPolygon([PolyWithoutSmallHoles(p, min_hole_area_km2)
                              for p in poly])
  else:
    raise ValueError('Input geometry is not a shapely Polygon or MultiPolygon') 

def _GeoJsonToShapelyGeometry(geometry):
  """Returns a |shapely| geometry from a GeoJSON geometry.

  Args:
    geometry: A dict or string representing a GeoJSON geometry.

  Raises:
    ValueError: If invalid GeoJSON geometry is passed.
  """
  if isinstance(geometry, basestring):
    geometry = json.loads(geometry)
  if not isinstance(geometry, dict) or 'type' not in geometry:
    raise ValueError('Invalid GeoJSON geometry.')

  if 'geometries' in geometry:
    return sgeo.GeometryCollection([_GeoJsonToShapelyGeometry(g)
                                    for g in geometry['geometries']])
  geometry = sgeo.shape(geometry)
  if isinstance(geometry, sgeo.Polygon) or isinstance(geometry, sgeo.MultiPolygon):
    geometry = geometry.buffer(0)
  return geometry 

def test_shrinks_polygon(self):
    with open(os.path.join(TEST_DIR, 'test_shrink.json'), 'r') as fd:
      ppa = json.load(fd)
    geometry = ppa['features'][0]['geometry']
    self.assertTrue(geometry['type'], 'Polygon')
    spoly = utils.ToShapely(geometry)
    mpoly = sgeo.MultiPolygon([sgeo.Point(0,0).buffer(1),
                               sgeo.Point(2,0).buffer(0.1)])
    poly1 = utils.ShrinkAndCleanPolygon(geometry, 1e-2)
    poly2 = utils.ShrinkAndCleanPolygon(spoly, 1e-2)
    self.assertTrue(poly2.area < spoly.area)
    with self.assertRaises(ValueError):
      poly = utils.ShrinkAndCleanPolygon(mpoly, 1e-2)
    self.assertEqual(poly1['type'], 'Polygon')
    self.assertTrue(isinstance(poly2, sgeo.Polygon))
    spoly1 = utils.ToShapely(poly1)
    self.assertEqual(poly2.difference(spoly1).area, 0)
    self.assertEqual(spoly1.difference(poly2).area, 0) 

def test_remove_small_holes(self):
    # Test all degenerate shapes (points, lines) have area zero
    # Use a small sw=quare centered around 60degree latitude, so radius is half
    # on the longitude
    square = sgeo.Polygon([(4, 59), (6, 59), (6, 61), (4, 61)])
    hole = sgeo.Polygon([(5, 60), (5.01, 60.01), (5, 60.02)])
    square_with_hole = sgeo.Polygon([(4, 59), (6, 59), (6, 61), (4, 61)],
                                    [[(5, 60), (5.01, 60.01), (5, 60.02)]])
    square_cleaned = utils.PolyWithoutSmallHoles(square_with_hole,
                                                 min_hole_area_km2= 0.62)
    multipoly = sgeo.MultiPolygon([square_with_hole, hole])
    multipoly_cleaned = utils.PolyWithoutSmallHoles(multipoly,
                                                    min_hole_area_km2= 0.62)
    area = utils.GeometryArea(square)
    area_with_hole = utils.GeometryArea(square_with_hole)
    self.assertAlmostEqual(area - area_with_hole, 0.617, 3)
    self.assertEqual(square, square_cleaned)
    self.assertEqual(multipoly_cleaned, sgeo.MultiPolygon([square, hole])) 

def test_cover_geometry_empty_geoms(tiler):
    """Empty geometries should return empty iterators."""
    assert not cover_geometry(tiler, geometry.Point(), 0) == True
    assert not cover_geometry(tiler, geometry.Point(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiPoint(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiPoint(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.LineString(), 0) == True
    assert not cover_geometry(tiler, geometry.LineString(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiLineString(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiLineString(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.Polygon(), 0) == True
    assert not cover_geometry(tiler, geometry.Polygon(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiPolygon(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiPolygon(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.GeometryCollection(), 0) == True
    assert not cover_geometry(tiler, geometry.GeometryCollection(), [0, 1]) == True 

def main():
    mm = 25.4
    w = 8.5 * mm
    h = 11 * mm
    p = 0.5 * mm
    states = MultiPolygon(load_shapes(STATES))
    # counties = MultiPolygon(load_shapes(COUNTIES))
    # counties = counties.simplify(100)
    g = xy.GCode.from_geometry(states)
    g = g.rotate(90)
    g = g.scale_to_fit(w, h, p).move(w / 2, p, 0.5, 0)
    # im = g.render(0, 0, w, h, 96 / mm)
    # im.write_to_png('usa.png')
    # g.save('usa.nc')
    device = xy.Device()
    device.home()
    device.gcode(g) 

def main():
    paths = []
    for path in PATHS:
        path = xy.parse_svg_path(path)
        path.append(path[0])
        paths.extend(path)
    polygons = [geometry.Polygon(x) for x in paths]
    lines = geometry.MultiPolygon(polygons)
    for i in range(4):
        n = 3 - i
        o = i * 10
        for j in range(-n, n + 1):
            paths += convert(lines.buffer(o + j * 0.667))
    drawing = xy.Drawing(paths).scale(1, -1).rotate_and_scale_to_fit(315, 380, step=90)
    im = drawing.render()
    im.write_to_png('frog.png')
    # xy.draw(drawing) 

def __init__(self, shape_list, crs, bbox_size):
        """
        :param shape_list: A list of geometrical shapes describing the area of interest
        :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon)
        :param crs: Coordinate reference system of the shapes in `shape_list`
        :type crs: CRS
        :param bbox_size: Physical size in metres of generated bounding boxes. Could be a float or tuple of floats
        :type bbox_size: int or (int, int) or float or (float, float)
        """
        super().__init__(shape_list, crs)

        self.bbox_size = self._parse_split_parameters(bbox_size, allow_float=True)

        self.shape_geometry = Geometry(self.area_shape, self.crs).transform(CRS.WGS84)

        self.utm_grid = self._get_utm_polygons()

        self._make_split() 

def test_results(self):
        for test_case in self.test_cases:
            with self.subTest(msg='Test case {}'.format(test_case.name)):
                iterator = test_case.request

                features = list(iterator)
                dates = iterator.get_dates()
                geometries = iterator.get_geometries()
                tiles = iterator.get_tiles()

                for result_list, expected_type in [(features, dict), (dates, datetime.datetime),
                                                   (geometries, MultiPolygon), (tiles, tuple)]:
                    self.assertEqual(len(result_list), test_case.result_len)
                    for result in result_list:
                        self.assertTrue(isinstance(result, expected_type),
                                        msg='Expected type {}, got type {}'.format(expected_type, type(result))) 

def GetPart90ExclusionZones():
  """Returns all Part90 federal exclusion zones as a |shapely.MultiPolygon|."""
  _GetAllExclusionZones()
  return _exclusion_zones_p90 

def project_geometry(geometry, crs, to_latlong=False):
    """
    Project a shapely Polygon or MultiPolygon from WGS84 to UTM, or vice-versa

    Parameters
    ----------
    geometry : shapely Polygon or MultiPolygon
        the geometry to project
    crs : int
        the starting coordinate reference system of the passed-in geometry
    to_latlong : bool, optional
        if True, project from crs to WGS84, if False, project
        from crs to local UTM zone

    Returns
    -------
    geometry_proj, crs : tuple (projected Shapely geometry, crs of the
    projected geometry)
    """
    gdf = gpd.GeoDataFrame()
    gdf.crs = crs
    gdf.name = 'geometry to project'
    gdf['geometry'] = None
    gdf.loc[0, 'geometry'] = geometry
    gdf_proj = project_gdf(gdf, to_latlong=to_latlong)
    geometry_proj = gdf_proj['geometry'].iloc[0]
    return geometry_proj, gdf_proj.crs 

def test_togeojson_with_multi_polygon(self):
    loop1 = [[-95.0, 40.0], [-95.5, 40.5], [-95.5, 40.0], [-95.0, 40.0]]
    loop2 = [[-95.2, 40.2], [-95.3, 40.2], [-95.3, 40.3], [-95.2, 40.2]]
    mpoly = sgeo.MultiPolygon([sgeo.Polygon(loop1), sgeo.Polygon(loop2)])

    json_poly = json.loads(utils.ToGeoJson(mpoly))
    self.assertIn('type', json_poly)
    self.assertEqual(json_poly['type'], 'MultiPolygon')
    self.assertIn('coordinates', json_poly)
    self.assertTrue(json_poly['coordinates'] == [[loop1], [list(reversed(loop2))]]) 

def _force_geometry_ccw(geometry):
    if geometry["type"] == "Polygon":
        return _force_polygon_ccw(geometry)
    elif geometry["type"] == "MultiPolygon":
        oriented_polygons = [
            _force_polygon_ccw({"type": "Polygon", "coordinates": g})
            for g in geometry["coordinates"]
        ]
        geometry["coordinates"] = [g["coordinates"] for g in oriented_polygons]
        return geometry
    else:
        return geometry 

def test_complex_polygons(self):
    square = sgeo.Polygon([(4, 59), (6, 59), (6, 61), (4, 61)])
    hole = sgeo.Polygon([(5, 60), (5.1, 60.1), (5, 60.2)])
    square_with_hole = sgeo.Polygon(square.exterior, [hole.exterior])
    multipoly = sgeo.MultiPolygon([square_with_hole, hole])
    self.assertAlmostEqual(utils.GeometryArea(square_with_hole),
                           utils.GeometryArea(square) - utils.GeometryArea(hole), 7)
    self.assertAlmostEqual(utils.GeometryArea(multipoly), 24699.71, 2) 

def GetUsBorder():
  """Gets the US border as a |shapely.MultiPolygon|.

  This is a composite US border for simulation purposes only.
  """
  global _border_zone
  if _border_zone is None:
    kml_file = os.path.join(CONFIG.GetFccDir(), USBORDER_FILE)
    zones = _ReadKmlZones(kml_file)
    _border_zone = ops.unary_union(zones.values())
  return _border_zone 

def multipolygon():
    polygon_1 = geometry.Polygon([[0, 0], [0, 4], [4, 4], [4, 0]])
    polygon_2 = geometry.Polygon([[5, 5], [5, 9], [9, 9], [9, 5]])

    return geometry.MultiPolygon([polygon_1, polygon_2]) 

def _extract_polygons_from_input_geometry(self):
        if isinstance(self._input_geometry, MultiPolygon):
            return (polygon for polygon in self._input_geometry)
        else:
            return (self._input_geometry,) 

def __init__(self):
        # Units (in or mm)
        self.units = Geometry.defaults["init_units"]
        
        # Final geometry: MultiPolygon or list (of geometry constructs)
        self.solid_geometry = None

        # Attributes to be included in serialization
        self.ser_attrs = ['units', 'solid_geometry']

        # Flattened geometry (list of paths only)
        self.flat_geometry = [] 

def input_geometry_is_valid(self):
        """Input geometry is of a :py:class:`shapely.geometry.Polygon`
        or a :py:class:`shapely.geometry.MultiPolygon`.

        :return: geometry is valid
        :rtype: bool
        """
        if isinstance(self._input_geometry, Polygon) or isinstance(
            self._input_geometry, MultiPolygon
        ):
            return True
        else:
            return False 

def _check_result(self, tesselation, orig_gdf, unique_id):
        """
        Check whether result of tessellation matches buildings and contains only Polygons.
        """
        # check against input layer
        ids_original = list(orig_gdf[unique_id])
        ids_generated = list(tesselation[unique_id])
        if len(ids_original) != len(ids_generated):
            import warnings

            self.collapsed = set(ids_original).difference(ids_generated)
            warnings.warn(
                "Tessellation does not fully match buildings. {len} element(s) collapsed "
                "during generation - unique_id: {i}".format(
                    len=len(self.collapsed), i=self.collapsed
                )
            )

        # check MultiPolygons - usually caused by error in input geometry
        self.multipolygons = tesselation[tesselation.geometry.type == "MultiPolygon"][
            unique_id
        ]
        if len(self.multipolygons) > 0:
            import warnings

            warnings.warn(
                "Tessellation contains MultiPolygon elements. Initial objects should be edited. "
                "unique_id of affected elements: {}".format(list(self.multipolygons))
            ) 

def _point_array(self, objects, unique_id):
        """
        Returns lists of points and ids based on geometry and unique_id.
        """
        points = []
        ids = []
        for idx, row in tqdm(objects.iterrows(), total=objects.shape[0]):
            if row["geometry"].type in ["Polygon", "MultiPolygon"]:
                poly_ext = row["geometry"].boundary
            else:
                poly_ext = row["geometry"]
            if poly_ext is not None:
                if poly_ext.type == "MultiLineString":
                    for line in poly_ext:
                        point_coords = line.coords
                        row_array = np.array(point_coords[:-1]).tolist()
                        for i, a in enumerate(row_array):
                            points.append(row_array[i])
                            ids.append(row[unique_id])
                elif poly_ext.type == "LineString":
                    point_coords = poly_ext.coords
                    row_array = np.array(point_coords[:-1]).tolist()
                    for i, a in enumerate(row_array):
                        points.append(row_array[i])
                        ids.append(row[unique_id])
                else:
                    raise Exception("Boundary type is {}".format(poly_ext.type))
        return points, ids 

def buffered_limit(gdf, buffer=100):
    """
    Define limit for :class:`momepy.Tessellation` as a buffer around buildings.

    See :cite:`fleischmann2020` for details.

    Parameters
    ----------
    gdf : GeoDataFrame
        GeoDataFrame containing building footprints
    buffer : float
        buffer around buildings limiting the extend of tessellation

    Returns
    -------
    MultiPolygon
        MultiPolygon or Polygon defining the study area

    Examples
    --------
    >>> limit = mm.buffered_limit(buildings_df)
    >>> type(limit)
    shapely.geometry.polygon.Polygon

    """
    study_area = gdf.copy()
    study_area["geometry"] = study_area.buffer(buffer)
    built_up = study_area.geometry.unary_union
    return built_up 

def _get_flat_poly_iter(cls, poly):
        if (isinstance(poly, shgeo.MultiPolygon) or
                isinstance(poly, shgeo.MultiLineString) or
                isinstance(poly, shgeo.GeometryCollection)):
            yield from poly
        else:
            yield poly 

def test_quadrat_cut_geometry(simple_polygon):
    multipolygon = load.quadrat_cut_geometry(geometry=simple_polygon,
                                             quadrat_width=0.5,
                                             min_num=3,
                                             buffer_amount=1e-9)

    assert isinstance(multipolygon, geometry.MultiPolygon)
    assert len(multipolygon) == 4 

def isolation_geometry(self, offset):
        """
        Creates contours around geometry at a given
        offset distance.

        :param offset: Offset distance.
        :type offset: float
        :return: The buffered geometry.
        :rtype: Shapely.MultiPolygon or Shapely.Polygon
        """
        return self.solid_geometry.buffer(offset) 

def _fix_geometry(geometry):
    shapely_geometry = shape(geometry)
    if not shapely_geometry.is_valid:
        buffered = shapely_geometry.buffer(0.0)
        # this will fix some invalid geometries, including bow-tie geometries, but for others it will return an empty geometry
        if buffered and (
            (type(buffered) == Polygon and buffered.exterior)
            or (type(buffered) == MultiPolygon and len(buffered.geoms) > 0)
        ):
            return mapping(buffered)
    return geometry