Python shapely.ops.unary_union() Examples

def _construct_centerline(self):
        vertices, ridges = self._get_voronoi_vertices_and_ridges()
        linestrings = []
        for ridge in ridges:
            if self._ridge_is_finite(ridge):
                starting_point = self._create_point_with_restored_coordinates(
                    x=vertices[ridge[0]][0], y=vertices[ridge[0]][1]
                )
                ending_point = self._create_point_with_restored_coordinates(
                    x=vertices[ridge[1]][0], y=vertices[ridge[1]][1]
                )
                linestring = LineString((starting_point, ending_point))

                if self._linestring_is_within_input_geometry(linestring):
                    linestrings.append(linestring)

        if len(linestrings) < 2:
            raise exceptions.TooFewRidgesError

        return unary_union(linestrings) 

def geom(self, union=False):
        """
        Converts the Rectangles to a shapely geometry.

        Parameters
        ----------
        union: boolean (default=False)
            Whether to compute a union between the geometries

        Returns
        -------
        A shapely geometry
        """
        boxes = [box(*g) for g in self.array([0, 1, 2, 3])]
        nboxes = len(boxes)
        if not nboxes:
            geom = GeometryCollection()
        elif nboxes == 1:
            geom = boxes[0]
        else:
            geom = MultiPolygon(boxes)
        return unary_union(geom) if union else geom 

def geom(self, union=False):
        """
        Converts the Contours to a shapely geometry.

        Parameters
        ----------
        union: boolean (default=False)
            Whether to compute a union between the geometries

        Returns
        -------
        A shapely geometry
        """
        geoms = expand_geoms([g['geometry'] for g in path_to_geom_dicts(self)])
        ngeoms = len(geoms)
        if not ngeoms:
            geom = GeometryCollection()
        elif ngeoms == 1:
            geom = geoms[0]
        else:
            geom = MultiLineString(geoms)
        return unary_union(geom) if union else geom 

def geom(self, union=False):
        """
        Converts the Path to a shapely geometry.

        Parameters
        ----------
        union: boolean (default=False)
            Whether to compute a union between the geometries

        Returns
        -------
        A shapely geometry
        """
        geoms = expand_geoms([g['geometry'] for g in path_to_geom_dicts(self)])
        ngeoms = len(geoms)
        if not ngeoms:
            geom = GeometryCollection()
        elif ngeoms == 1:
            geom = geoms[0]
        else:
            geom = MultiLineString(geoms)
        return unary_union(geom) if union else geom 

def geom(self, union=False):
        """
        Converts the Points to a shapely geometry.

        Parameters
        ----------
        union: boolean (default=False)
            Whether to compute a union between the geometries

        Returns
        -------
        A shapely geometry
        """
        points = [Point(x, y) for (x, y) in self.array([0, 1])]
        npoints = len(points)
        if not npoints:
            geom = GeometryCollection()
        elif len(points) == 1:
            geom = points[0]
        else:
            geom = MultiPoint(points)
        return unary_union(geom) if union else geom 

def _create_border(self, geometry: HybridGeometry, width, append=None):
        altitude = (np.vstack(chain(*(mesh.tolist() for mesh in geometry.faces)))[:, :, 2].max()+1)/1000
        geometry = self.buffered_bbox.intersection(geometry.geom)

        lines = tuple(chain(*(
            ((geom.exterior, *geom.interiors) if isinstance(geom, Polygon) else (geom,))
            for geom in getattr(geometry, 'geoms', (geometry,))
        )))

        if not lines:
            return np.empty((0, 3, 3+len(append)))

        lines = unary_union(lines).buffer(width, cap_style=CAP_STYLE.flat, join_style=JOIN_STYLE.mitre)

        vertices, faces = triangulate_polygon(lines)
        triangles = np.dstack((vertices[faces], np.full((faces.size, 1), fill_value=altitude).reshape((-1, 3, 1))))

        return self._append_to_vertices(triangles.astype(np.float32), append) 

def getBridgesPoly(o):
    if not hasattr(o, 'bridgespolyorig'):
        bridgecollectionname = o.bridges_collection_name
        bridgecollection = bpy.data.collections[bridgecollectionname]
        shapes = []
        bpy.ops.object.select_all(action='DESELECT')

        for ob in bridgecollection.objects:
            if ob.type == 'CURVE':
                ob.select_set(state=True)
        bpy.context.view_layer.objects.active = ob
        bpy.ops.object.duplicate();
        bpy.ops.object.join()
        ob = bpy.context.active_object
        shapes.extend(curveToShapely(ob, o.use_bridge_modifiers))
        ob.select_set(state=True)
        bpy.ops.object.delete(use_global=False)
        bridgespoly = sops.unary_union(shapes)

        # buffer the poly, so the bridges are not actually milled...
        o.bridgespolyorig = bridgespoly.buffer(distance=o.cutter_diameter / 2.0)
        o.bridgespoly_boundary = o.bridgespolyorig.boundary
        o.bridgespoly_boundary_prep = prepared.prep(o.bridgespolyorig.boundary)
        o.bridgespoly = prepared.prep(o.bridgespolyorig) 

def geom(self, union=False):
        """
        Converts the Path to a shapely geometry.

        Parameters
        ----------
        union: boolean (default=False)
            Whether to compute a union between the geometries

        Returns
        -------
        A shapely geometry
        """
        geoms = expand_geoms([g['geometry'] for g in polygons_to_geom_dicts(self)])
        ngeoms = len(geoms)
        if not ngeoms:
            geom = GeometryCollection()
        elif ngeoms == 1:
            geom = geoms[0]
        else:
            geom = MultiPolygon(geoms)
        return unary_union(geom) if union else geom 

def _union_in_blocks(contours, block_size):
    """
    Generator which yields a valid shape for each block_size multiple of
    input contours. This merges together the contours for each block before
    yielding them.
    """

    n_contours = len(contours)
    for i in range(0, n_contours, block_size):
        j = min(i + block_size, n_contours)

        inners = []
        for c in contours[i:j]:
            p = _contour_to_poly(c)
            if p.type == 'Polygon':
                inners.append(p)
            elif p.type == 'MultiPolygon':
                inners.extend(p.geoms)
        holes = unary_union(inners)
        assert holes.is_valid

        yield holes 

def _GetAllExclusionZones():
  """Read all exclusion zones."""
  global _exclusion_zones_gbs
  global _exclusion_zones_p90
  if _exclusion_zones_gbs is None:
    kml_file = os.path.join(CONFIG.GetNtiaDir(), EXCLUSION_ZONE_FILE)
    zones = _ReadKmlZones(kml_file, data_fields=['freqRangeMhz'])
    gbs_zones = []
    p90_zones = []
    for name, zone in zones.items():
      freq_range = _SplitFreqRange(zone.freqRangeMhz)
      if (3550, 3650) in freq_range:
        gbs_zones.append(zone.geometry)
      elif (3650, 3700) in freq_range:
        p90_zones.append(zone.geometry)
      else:
        raise ValueError('Zone %s: unsupported freq range %r',
                         name, freq_range)
    _exclusion_zones_gbs = ops.unary_union(gbs_zones)
    _exclusion_zones_p90 = ops.unary_union(p90_zones) 

def test_snap_bounds_to_zoom():
    bounds = (-180, -90, -60, -30)
    for pixelbuffer in [0, 5, 10]:
        for metatiling in [1, 2, 4]:
            pyramid = BufferedTilePyramid(
                "geodetic", pixelbuffer=pixelbuffer, metatiling=metatiling
            )
            for zoom in range(3, 5):
                snapped_bounds = mapchete.config.snap_bounds(
                    bounds=bounds,
                    pyramid=pyramid,
                    zoom=zoom
                )
                control_bounds = unary_union([
                    t.bbox for t in pyramid.tiles_from_bounds(bounds, zoom)
                ]).bounds
                assert snapped_bounds == control_bounds 

def minimal_set(polys):
    """Remove overlaps from a set of polygons.

    :param polys: List of polygons.
    :returns: List of polygons with no overlaps.
    """
    normal = polys[0].normal_vector
    as_2d = [p.project_to_2D() for p in polys]
    as_shapely = [Polygon(p) for p in as_2d]
    lines = [p.boundary for p in as_shapely]
    borders = unary_union(lines)
    shapes = [Polygon2D(p.boundary.coords) for p in polygonize(borders)]
    as_3d = [p.project_to_3D(polys[0]) for p in shapes]
    if not almostequal(as_3d[0].normal_vector, normal):
        as_3d = [p.invert_orientation() for p in as_3d]
    return [p for p in as_3d if p.area > 0] 

def _union_in_blocks(contours, block_size):
    """
    Generator which yields a valid shape for each block_size multiple of
    input contours. This merges together the contours for each block before
    yielding them.
    """

    n_contours = len(contours)
    for i in range(0, n_contours, block_size):
        j = min(i + block_size, n_contours)

        inners = []
        for c in contours[i:j]:
            p = _contour_to_poly(c)
            if p.type == 'Polygon':
                inners.append(p)
            elif p.type == 'MultiPolygon':
                inners.extend(p.geoms)
        holes = unary_union(inners)
        assert holes.is_valid

        yield holes 

def ToShapely(geometry):
  """Returns a |shapely| geometry from a generic geometry or a GeoJSON.

  The original geometry structure is maintained, for example GeometryCollections
  and MultiPolygons are kept as is. To dissolve them, apply the
  `shapely.ops.unary_union()` routine on the output.

  Args:
    geometry: A generic geometry or a GeoJSON geometry (dict or str). A generic
      geometry is any object implementing the __geo_interface__ supported by all
      major geometry libraries (shapely, ..)
  """
  if isinstance(geometry, sgeo.base.BaseGeometry):
    return geometry
  try:
    return _GeoJsonToShapelyGeometry(geometry.__geo_interface__)
  except AttributeError:
    return _GeoJsonToShapelyGeometry(geometry) 

def silhoueteOffset(context, offset):
    bpy.context.scene.cursor.location = (0, 0, 0)
    ob = bpy.context.active_object
    if ob.type == 'CURVE' or ob.type == 'FONT':
        silhs = curveToShapely(ob)
    else:
        silhs = getObjectSilhouete('OBJECTS', [ob])

    polys = []
    mp = shapely.ops.unary_union(silhs)
    mp = mp.buffer(offset, resolution=64)
    shapelyToCurve('offset curve', mp, ob.location.z)

    return {'FINISHED'} 

def absorb_by_shared_perimeter(sources, targets, relative_threshold=None):
    if len(sources) == 0:
        return targets

    if len(targets) == 0:
        raise IndexError("targets must be nonempty")

    inters = intersections(sources, targets, area_cutoff=None).buffer(0)
    assignment = assign_to_max(inters.length)

    if relative_threshold is not None:
        under_threshold = (
            sources.area / assignment.map(targets.area)
        ) < relative_threshold
        assignment = assignment[under_threshold]

    sources_to_absorb = GeoSeries(
        sources.groupby(assignment).apply(unary_union), crs=sources.crs,
    )

    result = targets.union(sources_to_absorb)

    # The .union call only returns the targets who had a corresponding
    # source to absorb. Now we fill in all of the unchanged targets.
    result = result.reindex(targets.index)
    did_not_absorb = result.isna() | result.is_empty
    result.loc[did_not_absorb] = targets[did_not_absorb]

    return result 

def holes_of_union(geometries):
    """Returns any holes in the union of the given geometries."""
    geometries = get_geometries(geometries)
    if not all(
        isinstance(geometry, (Polygon, MultiPolygon)) for geometry in geometries
    ):
        raise TypeError("all geometries must be Polygons or MultiPolygons")

    union = unary_union(geometries)
    series = holes(union)
    series.crs = geometries.crs
    return series 

def hybrid_union(geoms):
    if not geoms:
        return HybridGeometry(GeometryCollection(), ())
    if len(geoms) == 1:
        return geoms[0]
    add_faces = {}
    for other in geoms:
        for crop_id, faces in other.add_faces.items():
            add_faces[crop_id] = add_faces.get(crop_id, ()) + faces
    return HybridGeometry(geom=unary_union(tuple(geom.geom for geom in geoms)),
                          faces=tuple(chain(*(geom.faces for geom in geoms))),
                          add_faces=add_faces,
                          crop_ids=reduce(operator.or_, (other.crop_ids for other in geoms), set())) 

def clip_altitudes(self, new_geometry, new_altitude=None):
        # register new geometry with an altitude
        # a geometry with no altitude will reset the altitude information of its area as if nothing was ever there
        if self.last_altitude is not None and self.last_altitude > new_altitude:
            raise ValueError('Altitudes have to be ascending.')

        if new_altitude in self.altitudes:
            self.altitudes[new_altitude] = unary_union([self.altitudes[new_altitude], new_geometry])
        else:
            self.altitudes[new_altitude] = new_geometry 

def save(self, last_update, new_update):
        self.finalize()

        for level_id, geometries in self._geometries_by_level.items():
            geometries = unary_union(geometries)
            if geometries.is_empty:
                continue
            history = MapHistory.open_level(level_id, mode='base', default_update=last_update)
            history.add_geometry(geometries.buffer(1), new_update)
            history.save()
        self.reset() 

def _get_unary_union(self, level_id):
        union = self._unary_unions.get(level_id)
        if union is None:
            union = unary_union(self._geometries_by_level[level_id])
            self._unary_unions[level_id] = union
        return union 

def copdem_zone(lon_ex, lat_ex):
    """Returns a list of Copernicus DEM tarfile and tilename tuples
    """

    # path to the lookup shapefiles
    gdf = gpd.read_file(get_demo_file('RGI60_COPDEM_lookup.shp'))

    # intersect with lat lon extents
    p = _extent_to_polygon(lon_ex, lat_ex, to_crs=gdf.crs)
    gdf = gdf.loc[gdf.intersects(p)]

    # COPDEM is global, if we miss any tile it is worth an error
    if (len(gdf) == 0) or (not unary_union(gdf.geometry).contains(p)):
        raise InvalidDEMError('Could not find all necessary Copernicus DEM '
                              'tiles. This should not happen in a global DEM. '
                              'Check the RGI-CopernicusDEM lookup shapefile '
                              'for this particular glacier!')

    flist = []
    for _, g in gdf.iterrows():
        cpp = g['CPP File']
        eop = g['Eop Id']
        eop = eop.split(':')[-2]
        assert 'Copernicus' in eop

        flist.append((cpp, eop))

    return flist 

def getObjectOutline(radius, o, Offset):  # FIXME: make this one operation independent
    # circle detail, optimize, optimize thresold.

    polygons = getOperationSilhouete(o)

    i = 0
    # print('offseting polygons')

    if Offset:
        offset = 1
    else:
        offset = -1

    outlines = []
    i = 0
    # print(polygons, polygons.type)
    for p1 in polygons:  # sort by size before this???
        print(p1.type, len(polygons))
        i += 1
        if radius > 0:
            p1 = p1.buffer(radius * offset, resolution=o.circle_detail)
        outlines.append(p1)

    print(outlines)
    if o.dont_merge:
        outline = sgeometry.MultiPolygon(outlines)
    # for ci in range(0,len(p)):
    #	outline.addContour(p[ci],p.isHole(ci))
    else:
        # print(p)
        outline = shapely.ops.unary_union(outlines)
    # outline = sgeometry.MultiPolygon([outline])
    # shapelyToCurve('oboutline',outline,0)
    return outline 

def getAmbient(o):
    if o.update_ambient_tag:
        if o.ambient_cutter_restrict:  # cutter stays in ambient & limit curve
            m = o.cutter_diameter / 2
        else:
            m = 0

        if o.ambient_behaviour == 'AROUND':
            r = o.ambient_radius - m
            o.ambient = getObjectOutline(r, o, True)  # in this method we need ambient from silhouete
        else:
            o.ambient = spolygon.Polygon(((o.min.x + m, o.min.y + m), (o.min.x + m, o.max.y - m),
                                          (o.max.x - m, o.max.y - m), (o.max.x - m, o.min.y + m)))

        if o.use_limit_curve:
            if o.limit_curve != '':
                limit_curve = bpy.data.objects[o.limit_curve]
                # polys=curveToPolys(limit_curve)
                polys = curveToShapely(limit_curve)
                o.limit_poly = shapely.ops.unary_union(polys)
                # for p in polys:
                #	o.limit_poly+=p
                if o.ambient_cutter_restrict:
                    o.limit_poly = o.limit_poly.buffer(o.cutter_diameter / 2, resolution=o.circle_detail)
            o.ambient = o.ambient.intersection(o.limit_poly)
    o.update_ambient_tag = False 

def get_bounds(feature_collection):
    """Get a bounding box for a FeatureCollection"""
    shape_lst = [shape(f['geometry']) for f in feature_collection['features']]
    return unary_union(shape_lst).bounds 

def draw_v_pair(pair, strip, l):
    # Drow polligons between two masks from test
    if l > 0:
        roi_points = [(0, 0), (img_side_len, 0),
                      (img_side_len*l, img_side_len), (0, img_side_len)]
        roi_poly = Polygon(roi_points)
        top_tile = find_points(pair[0])
        bottom_tile = find_points(pair[1])
        top_tile_pos, top_tile_neg = top_tile
        bottom_tile_pos, bottom_tile_neg = bottom_tile
        v_shift = l * img_side_len

        square_points = [(0, 0), (img_side_len, 0), (img_side_len, v_shift), (0, v_shift)]
        square_poly = Polygon(square_points)
        lines = []
        for i in range(len(top_tile_pos[bottom])):
            line = LineString([(top_tile_pos[bottom][i], 0),
                               (bottom_tile_pos[top][i], v_shift),
                               (bottom_tile_neg[top][i], v_shift),
                               (top_tile_neg[bottom][i], 0)])
            lines.append(line)

        merged = linemerge([square_poly.boundary, *lines])
        borders = unary_union(merged)
        polygons = []
        for poly in polygonize(borders):
            polygons.append(poly)
        masks = [mask_for_polygons([polygons[i]], (v_shift, img_side_len))
                 for i in range(0, len(polygons), 2)]
        mask = (np.any(masks, axis=0)*255).astype(np.uint8)
        return mask
    return None 

def make_clusters_geojson(self):
        """Makes a feature collection of aggregated feature geometries """
        geojson = LastUpdatedOrderedDict()
        geojson['type'] = 'FeatureCollection'
        geojson['features'] = []
        for cluster_group, feat_geoms in self.raw_features.items():
            geoms = []
            for feat_geometry in feat_geoms:
                geom_raw = shape(feat_geometry)
                geom = geom_raw.buffer(self.buffer_tolerance)
                if geom.is_valid:
                    geoms.append(geom)
            if len(geoms) < 1:
                # No valid geometries for this cluster, so skip
                continue
            union_geom = unary_union(geoms)
            poly_union = union_geom.convex_hull
            poly_union_simple = poly_union.simplify(self.simplify_tolerance)
            feat = LastUpdatedOrderedDict()
            feat['type'] = 'Feature'
            feat['properties'] = LastUpdatedOrderedDict()
            feat['properties'][self.cluster_property] = cluster_group
            feat['properties'][self.cluster_num_features_prop] = len(geoms)
            feat['geometry'] = LastUpdatedOrderedDict()
            geometry_type = poly_union_simple.geom_type
            coordinates = [list(poly_union_simple.exterior.coords)]
            v_geojson = ValidateGeoJson()
            c_ok = v_geojson.validate_all_geometry_coordinates(geometry_type,
                                                               coordinates)
            if not c_ok:
                coordinates = v_geojson.fix_geometry_rings_dir(geometry_type, coordinates)
            feat['geometry']['type'] = geometry_type
            feat['geometry']['coordinates'] = coordinates
            centroid = self.get_feature_centroid(feat['geometry'])
            feat['properties']['latitude'] = centroid[1]
            feat['properties']['longitude'] = centroid[0]
            geojson['features'].append(feat)
        return geojson 

def compute_bb(self) -> List[float]:
        """Compute the bounding box of all of the parts in the geometry.

        Returns
        -------
        List of [min_x, max_x, min_y, max_y].

        """
        all_shapes = list(self.parts.values())
        bbox_vertices = unary_union(all_shapes).envelope.exterior.coords.xy
        min_x = min(bbox_vertices[0])
        max_x = max(bbox_vertices[0])
        min_y = min(bbox_vertices[1])
        max_y = max(bbox_vertices[1])
        return [min_x, max_x, min_y, max_y] 

def GetUrbanAreas(simplify_deg=1e-3):
  """Gets the US urban area as a |shapely.GeometryCollection|.

  Note: Client code should cache it as expensive to load (and not cached here).

  Args:
    simplify_deg: if defined, simplify the zone with given tolerance (degrees).
      Default is 1e-3 which corresponds roughly to 100m in continental US.
  """
  kml_file = os.path.join(CONFIG.GetNtiaDir(), URBAN_AREAS_FILE)
  zones = _ReadKmlZones(kml_file, root_id_zone='Document', simplify=simplify_deg)
  urban_areas = sgeo.GeometryCollection(zones.values())  # ops.unary_union(zones.values())
  return urban_areas 

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