Python shapely.ops.cascaded_union() Examples

def get_point_in_polygon_score(self, lane_seq: List[int],
                                   xy_seq: np.ndarray, city_name: str,
                                   avm: ArgoverseMap) -> int:
        """Get the number of coordinates that lie insde the lane seq polygon.

        Args:
            lane_seq: Sequence of lane ids
            xy_seq: Trajectory coordinates
            city_name: City name (PITT/MIA)
            avm: Argoverse map_api instance
        Returns:
            point_in_polygon_score: Number of coordinates in the trajectory that lie within the lane sequence

        """
        lane_seq_polygon = cascaded_union([
            Polygon(avm.get_lane_segment_polygon(lane, city_name)).buffer(0)
            for lane in lane_seq
        ])
        point_in_polygon_score = 0
        for xy in xy_seq:
            point_in_polygon_score += lane_seq_polygon.contains(Point(xy))
        return point_in_polygon_score 

def export_svg(self, scale_factor=0.00):
        """
        Exports the Gemoetry Object as a SVG Element

        :return: SVG Element
        """
        # Make sure we see a Shapely Geometry class and not a list
        geom = cascaded_union(self.flatten())

        # scale_factor is a multiplication factor for the SVG stroke-width used within shapely's svg export

        # If 0 or less which is invalid then default to 0.05
        # This value appears to work for zooming, and getting the output svg line width
        # to match that viewed on screen with FlatCam
        if scale_factor <= 0:
            scale_factor = 0.05

        # Convert to a SVG
        svg_elem = geom.svg(scale_factor=scale_factor)
        return svg_elem 

def subtract_polygon(self, points):
        """
        Subtract polygon from the given object. This only operates on the paths in the original geometry, i.e. it converts polygons into paths.

        :param points: The vertices of the polygon.
        :return: none
        """
        if self.solid_geometry is None:
            self.solid_geometry = []

        #pathonly should be allways True, otherwise polygons are not subtracted
        flat_geometry = self.flatten(pathonly=True)
        log.debug("%d paths" % len(flat_geometry))
        polygon=Polygon(points)
        toolgeo=cascaded_union(polygon)
        diffs=[]
        for target in flat_geometry:
            if type(target) == LineString or type(target) == LinearRing:
                diffs.append(target.difference(toolgeo))
            else:
                log.warning("Not implemented.")
        self.solid_geometry=cascaded_union(diffs) 

def buffer(self, buf_distance):
        selected = self.get_selected()

        if len(selected) == 0:
            self.app.inform.emit("[warning] Nothing selected for buffering.")
            return

        if not isinstance(buf_distance, float):
            self.app.inform.emit("[warning] Invalid distance for buffering.")
            return

        pre_buffer = cascaded_union([t.geo for t in selected])
        results = pre_buffer.buffer(buf_distance)
        self.add_shape(DrawToolShape(results))

        self.replot() 

def union(self):
        """
        Makes union of selected polygons. Original polygons
        are deleted.

        :return: None.
        """

        results = cascaded_union([t.geo for t in self.get_selected()])

        # Delete originals.
        for_deletion = [s for s in self.get_selected()]
        for shape in for_deletion:
            self.delete_shape(shape)

        # Selected geometry is now gone!
        self.selected = []

        self.add_shape(DrawToolShape(results))

        self.replot() 

def cutpath(self):
        selected = self.get_selected()
        tools = selected[1:]
        toolgeo = cascaded_union([shp.geo for shp in tools])

        target = selected[0]
        if type(target.geo) == Polygon:
            for ring in poly2rings(target.geo):
                self.add_shape(DrawToolShape(ring.difference(toolgeo)))
            self.delete_shape(target)
        elif type(target.geo) == LineString or type(target.geo) == LinearRing:
            self.add_shape(DrawToolShape(target.geo.difference(toolgeo)))
            self.delete_shape(target)
        else:
            self.app.log.warning("Not implemented.")

        self.replot() 

def alpha_shape(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    """
    if len(points) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull

    coords = np.array([point.coords[0] for point in points])
    tri = Delaunay(coords)
    triangles = coords[tri.vertices]
    a = ((triangles[:,0,0] - triangles[:,1,0]) ** 2 + (triangles[:,0,1] - triangles[:,1,1]) ** 2) ** 0.5
    b = ((triangles[:,1,0] - triangles[:,2,0]) ** 2 + (triangles[:,1,1] - triangles[:,2,1]) ** 2) ** 0.5
    c = ((triangles[:,2,0] - triangles[:,0,0]) ** 2 + (triangles[:,2,1] - triangles[:,0,1]) ** 2) ** 0.5
    s = ( a + b + c ) / 2.0
    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
    circums = a * b * c / (4.0 * areas)
    filtered = triangles[circums < (1.0 / alpha)]
    edge1 = filtered[:,(0,1)]
    edge2 = filtered[:,(1,2)]
    edge3 = filtered[:,(2,0)]
    edge_points = np.unique(np.concatenate((edge1,edge2,edge3)), axis = 0).tolist()
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return cascaded_union(triangles), edge_points 

def _get_merged_polygon(self, cluster):
        """Merge polygons using shapely (Internal).

        Given a single cluster from _get_merge_clusters_from_df(), This creates
        and merges polygons into a single cascaded union. It first dilates the
        polygons by buffer_size pixels to make them overlap, merges them,
        then erodes back by buffer_size to get the merged polygon.

        """
        buffer_size = self.merge_thresh + 3
        nest_polygons = []
        for nestinfo in cluster:
            nest = dict(self.edge_contours[nestinfo['roiname']].loc[
                nestinfo['nid'], :])
            roitop = self.roiinfos[nestinfo['roiname']]['top']
            roileft = self.roiinfos[nestinfo['roiname']]['left']
            coords = _parse_annot_coords(
                nest, x_offset=roileft, y_offset=roitop)
            nest_polygons.append(Polygon(coords).buffer(buffer_size))
        merged_polygon = cascaded_union(nest_polygons).buffer(-buffer_size)
        return merged_polygon

    # %% ===================================================================== 

def cascaded_union_with_alt(polyalt: AirspaceList) -> AirspaceList:
    altitudes = set(alt for _, *low_up in polyalt for alt in low_up)
    slices = sorted(altitudes)
    if len(slices) == 1 and slices[0] is None:
        simple_union = cascaded_union([p for p, *_ in polyalt])
        return [ExtrudedPolygon(simple_union, float("-inf"), float("inf"))]
    results: List[ExtrudedPolygon] = []
    for low, up in zip(slices, slices[1:]):
        matched_poly = [
            p
            for (p, low_, up_) in polyalt
            if low_ <= low <= up_ and low_ <= up <= up_
        ]
        new_poly = ExtrudedPolygon(cascaded_union(matched_poly), low, up)
        if len(results) > 0 and new_poly.polygon.equals(results[-1].polygon):
            merged = ExtrudedPolygon(new_poly.polygon, results[-1].lower, up)
            results[-1] = merged
        else:
            results.append(new_poly)
    return results


# -- Methods below are placed here because of possible circular imports -- 

def get_names_rings(area):
    names = []
    all_rings =[]

    for idx, a in enumerate(area):
        rings = [LinearRing(p) for p in a['polygons']]
        if len(rings) > 1:
            u = cascaded_union(rings)
        else:
            u = rings[0]
        all_rings.append(u.envelope)
        if a.get('name'):
            names.append(a.get('name'))

    names = sorted(names)
    
    return names, all_rings 

def _area_at_zoom(self, zoom):
        if zoom not in self._cache_area_at_zoom:
            # use union of all input items and, if available, intersect with
            # init_bounds
            if "input" in self._params_at_zoom[zoom]:
                input_union = cascaded_union([
                    self.input[get_hash(v)].bbox(self.process_pyramid.crs)
                    for k, v in _flatten_tree(self._params_at_zoom[zoom]["input"])
                    if v is not None
                ])
                self._cache_area_at_zoom[zoom] = input_union.intersection(
                    box(*self.init_bounds)
                ) if self.init_bounds else input_union
            # if no input items are available, just use init_bounds
            else:
                self._cache_area_at_zoom[zoom] = box(*self.init_bounds)
        return self._cache_area_at_zoom[zoom] 

def offset(self, dist, side):
        if dist == 0:
            return _SidedPolygon(self.poly, self.level)
        if dist < 0:
            side = not side
            dist = abs(dist)
        if (side == c.OUTSIDE and self.isFeature) or (side == c.INSIDE and not self.isFeature):
            return _SidedPolygon(self.poly.buffer(dist), self.level)
        try:
            buffPoly = self.poly.exterior.buffer(dist)
            if len(buffPoly.interiors) > 1:
                inPoly = cascaded_union([Polygon(i) for i in buffPoly.interiors])            
            else:
                inPoly = Polygon(buffPoly.interiors[0])
            return _SidedPolygon(inPoly, self.level)
        except Exception:
            return None 

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 LinesToPolyHoles(lines):
    # get all of the polys, both outer and holes
    # if there are holes, you'll get them double:
    # once as themselves and again as holes in a boundary
    polys = list(polygonize(lines))

    # merge to get just the outer
    bounds =  cascaded_union(polys)

    if not iterable(bounds):
        bounds = [bounds]

    retval = []
    for bound in bounds:
        for poly in polys:
            if Polygon(poly.exterior).almost_equals(bound):
                retval.append(poly)

    return retval 

def geohashes_to_polygon(geohashes):
    """
    :param geohashes: array-like. List of geohashes to form resulting polygon.
    :return: shapely geometry. Resulting Polygon after combining geohashes.
    """
    return cascaded_union([geohash_to_polygon(g) for g in geohashes]) 

def test_voronoi_geopandas_with_plot():
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

    # focus on South America, convert to World Mercator (unit: meters)
    south_am = world[world.continent == 'South America'].to_crs(epsg=3395)
    cities = cities.to_crs(south_am.crs)  # convert city coordinates to same CRS!

    # create the bounding shape as union of all South American countries' shapes
    south_am_shape = cascaded_union(south_am.geometry)
    south_am_cities = cities[cities.geometry.within(south_am_shape)]  # reduce to cities in South America

    # convert the pandas Series of Point objects to NumPy array of coordinates
    coords = points_to_coords(south_am_cities.geometry)

    # calculate the regions
    poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, south_am_shape)

    assert isinstance(poly_shapes, list)
    assert 0 < len(poly_shapes) <= len(coords)
    assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes])

    assert np.array_equal(points_to_coords(pts), coords)

    assert isinstance(poly_to_pt_assignments, list)
    assert len(poly_to_pt_assignments) == len(poly_shapes)
    assert all([isinstance(assign, list) for assign in poly_to_pt_assignments])
    assert all([len(assign) == 1 for assign in poly_to_pt_assignments])   # in this case there is a 1:1 correspondance

    fig, ax = subplot_for_map()

    plot_voronoi_polys_with_points_in_area(ax, south_am_shape, poly_shapes, pts, poly_to_pt_assignments)

    return fig 

def ReadShapeFile(ShapeFile):

    """
    Open shapefile and create Polygon dictionary
    returns dictionary of shapely Polygons

    MDH

    """

    #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 

def create_geometry(self):
        # TODO: This takes forever. Too much data?
        self.solid_geometry = cascaded_union([geo['geom'] for geo in self.gcode_parsed]) 

def polygon_shapes_from_voronoi_lines(poly_lines, geo_shape=None, shapes_from_diff_with_min_area=None):
    """
    Form shapely Polygons objects from a list of shapely LineString objects in `poly_lines` by using
    [`polygonize`](http://toblerity.org/shapely/manual.html#shapely.ops.polygonize). If `geo_shape` is not None, then
    the intersection between any generated polygon and `geo_shape` is taken in case they overlap (i.e. the Voronoi
    regions at the border are "cut" to the `geo_shape` polygon that represents the geographic area holding the
    Voronoi regions). Setting `shapes_from_diff_with_min_area` fixes rare errors where the Voronoi shapes do not fully
    cover `geo_shape`. Set this to a small number that indicates the minimum valid area of "fill up" Voronoi region
    shapes.
    Returns a list of shapely Polygons objects.
    """

    # generate shapely Polygon objects from the LineStrings of the Voronoi shapes in `poly_lines`
    poly_shapes = []
    for p in polygonize(poly_lines):
        if geo_shape is not None and not geo_shape.contains(p):    # if `geo_shape` contains polygon `p`,
            p = p.intersection(geo_shape)                          # intersect it with `geo_shape` (i.e. "cut" it)

        if not p.is_empty:
            poly_shapes.append(p)

    if geo_shape is not None and shapes_from_diff_with_min_area is not None:
        # fix rare cases where the generated polygons of the Voronoi regions don't fully cover `geo_shape`
        vor_polys_union = cascaded_union(poly_shapes)   # union of Voronoi regions
        diff = np.array(geo_shape.difference(vor_polys_union), dtype=object)    # "gaps"
        diff_areas = np.array([p.area for p in diff])    # areas of "gaps"
        # use only those "gaps" bigger than `shapes_from_diff_with_min_area` because very tiny areas are generated
        # at the borders due to floating point errors
        poly_shapes.extend(diff[diff_areas >= shapes_from_diff_with_min_area])

    return poly_shapes 

def merge(boxes):
        polygons = [box.to_polygon() for box in boxes]
        union = cascaded_union(polygons)
        union_aabb = union.envelope
        return BBox.from_axis_aligned_bbox(AxisAlignedBBox(*union_aabb.bounds)) 

def bounds(self):
        return cascaded_union(tuple(item.geometry.buffer(0)
                                    for item in chain(self.altitudeareas.all(), self.buildings.all()))).bounds 

def union_extracts(extracts):
    unioned = cascaded_union([ make_polygon(p) for p in extracts ])

    if unioned.geom_type == 'Polygon':
        return [ polygon_to_extract(unioned) ]
    else:
        return [ polygon_to_extract(geom) for geom in unioned ] 

def _merge_polygons(self, poly_list):
        if self.buffer_size > 0:
            poly_list = [j.buffer(self.buffer_size) for j in poly_list]
            merged_polys = cascaded_union(poly_list).buffer(-self.buffer_size)
        else:
            merged_polys = cascaded_union(poly_list)
        return merged_polys

    # %% ===================================================================== 

def fix_multipolygon(mpoly):
    # try to fix it by buffering by 1 as close by contours are most likely to be of same object
    poly = []
    for i, item in enumerate(mpoly):
        p = item.buffer(1)   #maybe 0
        if p.is_valid:
            poly.append(p)
    poly = cascaded_union(poly).buffer(0)
    if not poly.geom_type == 'Polygon':
        poly = poly.convex_hull
    return poly 

def fix_poly(polys):
    poly = []
    for p in polys:
        p = p.buffer(0)
        if p.is_valid:
            poly.append(p)
    # may return empty poly check downstream
    poly = cascaded_union(poly).buffer(0)
    return poly


#polygon need to be shifted as its part of larger tile 

def area_at_zoom(self, zoom=None):
        """
        Return process bounding box for zoom level.

        Parameters
        ----------
        zoom : int or None
            if None, the union of all zoom level areas is returned

        Returns
        -------
        process area : shapely geometry
        """
        if not self._init_inputs:
            return box(*self.init_bounds)
        if zoom is None:
            if not self._cache_full_process_area:
                logger.debug("calculate process area ...")
                self._cache_full_process_area = cascaded_union([
                    self._area_at_zoom(z) for z in self.init_zoom_levels]
                ).buffer(0)
            return self._cache_full_process_area
        else:
            if zoom not in self.init_zoom_levels:
                raise ValueError("zoom level not available with current configuration")
            return self._area_at_zoom(zoom) 

def compute_bounds(all_rings):  
    if len(all_rings) > 1:
        combined = cascaded_union(all_rings)
    else:
        combined = all_rings[0]

    combined = combined.envelope
    bounds_not_wrapped = copy.deepcopy(combined.bounds)
    bounds=(bounds_not_wrapped[0], convert_to_bounded_lon(bounds_not_wrapped[1]), bounds_not_wrapped[2], convert_to_bounded_lon(bounds_not_wrapped[3]))
    return bounds 

def _get_level_geometries(level):
        buildings = level.buildings.all()
        buildings_geom = cascaded_union([building.geometry for building in buildings])
        spaces = {space.pk: space for space in level.spaces.all()}
        holes_geom = []
        for space in spaces.values():
            if space.outside:
                space.geometry = space.geometry.difference(buildings_geom)
            columns = [column.geometry for column in space.columns.all()]
            if columns:
                columns_geom = cascaded_union([column.geometry for column in space.columns.all()])
                space.geometry = space.geometry.difference(columns_geom)
            holes = [hole.geometry for hole in space.holes.all()]
            if holes:
                space_holes_geom = cascaded_union(holes)
                holes_geom.append(space_holes_geom.intersection(space.geometry))
                space.geometry = space.geometry.difference(space_holes_geom)

        for building in buildings:
            building.original_geometry = building.geometry

        if holes_geom:
            holes_geom = cascaded_union(holes_geom)
            holes_geom_prep = prepared.prep(holes_geom)
            for obj in buildings:
                if holes_geom_prep.intersects(obj.geometry):
                    obj.geometry = obj.geometry.difference(holes_geom)

        results = []
        results.extend(buildings)
        for door in level.doors.all():
            results.append(door)

        results.extend(spaces.values())
        return results 

def set_space_outside(apps, schema_editor):
    Section = apps.get_model('mapdata', 'Section')
    for section in Section.objects.all():
        building_geometries = cascaded_union(tuple(building.geometry for building in section.buildings.all()))
        for space in section.spaces.all():
            if space.geometry.intersection(building_geometries).area / space.geometry.area < 0.5:
                space.outside = True
                space.save() 

def union(self):
        """
        Runs a cascaded union on the list of objects in
        solid_geometry.

        :return: None
        """
        self.solid_geometry = [cascaded_union(self.solid_geometry)]