Python ops.cascaded_union函数代码示例

    def parse_query(self, query):
        super(MapPlotter, self).parse_query(query)

        self.projection = query.get('projection')

        self.area = query.get('area')

        names = []
        centroids = []
        all_rings = []
        data = None
        for idx, a in enumerate(self.area):
            if isinstance(a, basestring):
                a = a.encode("utf-8")
                sp = a.split('/', 1)
                if data is None:
                    data = list_areas(sp[0], simplify=False)

                b = [x for x in data if x.get('key') == a]
                a = b[0]
                self.area[idx] = a
            else:
                p = np.array(a['polygons'])
                p[:, :, 1] = pyresample.utils.wrap_longitudes(p[:, :, 1])
                a['polygons'] = p.tolist()
                del p

            rings = [LinearRing(po) for po in a['polygons']]
            if len(rings) > 1:
                u = cascaded_union(rings)
            else:
                u = rings[0]

            all_rings.append(u)
            if a.get('name'):
                names.append(a.get('name'))
                centroids.append(u.centroid)

        nc = sorted(zip(names, centroids))
        self.names = [n for (n, c) in nc]
        self.centroids = [c for (n, c) in nc]
        data = None

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

        self.combined_area = combined
        combined = combined.envelope

        self.centroid = list(combined.centroid.coords)[0]
        self.bounds = combined.bounds

        self.show_bathymetry = bool(query.get('bathymetry'))
        self.show_area = bool(query.get('showarea'))

        self.quiver = query.get('quiver')

        self.contour = query.get('contour')

    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 processpoly(polylist):
    print u"Полигоны"
    for p in polylist:
        pol = Polygon(p)
        print pol
    print u"Объединение"
    print cascaded_union([Polygon(p) for p in polylist])
    return 1

def fig_contains(fig1, fig2) -> float:
    pgz1 = polygonize_figure(fig1)
    pgz2 = polygonize_figure(fig2)
    if len(pgz1) > 0 and len(pgz2) > 0:
        plist1 = [Polygon(p) for p in pgz1]
        plist2 = [Polygon(p) for p in pgz2]
        u_polygon1 = cascaded_union(plist1)
        u_polygon2 = cascaded_union(plist2)
        if u_polygon1.contains(u_polygon2):
            return True
    return False

def handle_special_countries(c):
    c['Iceland'].patch = ops.cascaded_union([c['Iceland'].patch, c['England'].patch,
                                             c['Ireland'].patch]).convex_hull
    c['Madagascar'].patch = ops.cascaded_union([c['Madagascar'].patch,
                                             c['Mozambique'].patch]).convex_hull
    c['Madagascar'].__class__ = c['Iceland'].__class__
    c['Malaysia'].patch = c['Malaysia'].patch.buffer(1)
    c['Indonesia'].patch = c['Indonesia'].patch.buffer(1)
    c['Indonesia'].__class__ = c['Iceland'].__class__
    c['Malaysia'].__class__ = c['Iceland'].__class__
    c['Hawaii'].continent = ''
    c['Falkland Is.'].continent = ''

def fig_overlap_area(fig1, fig2) -> float:
    pgzfg1 = polygonize_figure(fig1)
    pgzfg2 = polygonize_figure(fig2)
    if len(pgzfg1) > 0 and len(pgzfg2) > 0:
        plist1 = [Polygon(p) for p in pgzfg1]
        plist2 = [Polygon(p) for p in pgzfg2]
        u_polygon1 = cascaded_union(plist1)
        u_polygon2 = cascaded_union(plist2)
        area = u_polygon1.intersection(u_polygon2).area
    else:
        area = 0.0
    return area

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])

    coords = np.array(points)

    print(coords)

    tri = Delaunay(coords)
    edges = set()
    edge_points = []
    # loop over triangles:
    # ia, ib, ic = indices of corner points of the
    # triangle
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]
        # Lengths of sides of triangle
        a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
        b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
        c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
        # Semiperimeter of triangle
        s = (a + b + c)/2.0
        # Area of triangle by Heron's formula
        area = math.sqrt(s*(s-a)*(s-b)*(s-c))
        circum_r = a*b*c/(4.0*area)
        # Here's the radius filter.
        #print circum_r
        if circum_r < 1.0/alpha:
            add_edge(edges, edge_points, coords, ia, ib)
            add_edge(edges, edge_points, coords, ib, ic)
            add_edge(edges, edge_points, coords, ic, ia)
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return (cascaded_union(triangles), edge_points)

    print (cascaded_union(triangles), edge_points)

def ConvertPolygons2LatLong(PolygonDict, CRS):

    """
    Convert coordinates to latlong using pyproj
    https://gis.stackexchange.com/questions/127427/transforming-shapely-polygon-and-multipolygon-objects

    returns shapely Polygon/Multipolygon Dictionary

    MDH

    """

    # define the output coordinate system
    Output_CRS = Proj({'init': "epsg:4326"})

    NewPolygonDict = {}

    #loop through polygon dict and convert to latlong
    for Key, Poly in PolygonDict.iteritems():

        if Poly.geom_type == 'Polygon':
            # get x and y as arrays
            x,y = Poly.exterior.coords.xy
            lon,lat = transform(CRS, Output_CRS, x, y)
            Shape = Polygon(zip(lon,lat))

            #check for multipolygons
            if Key in NewPolygonDict:
                Polygons = [Shape, NewPolygonDict[Key]]
                Shape = cascaded_union(Polygons)

            #add converted polygon to new dictionary
            NewPolygonDict[Key] = Shape

        elif Poly.geom_type == 'MultiPolygon':
            for SinglePoly in Poly:
                # get x and y as arrays
                x,y = SinglePoly.exterior.coords.xy
                lon,lat = transform(CRS, Output_CRS, x, y)
                Shape = Polygon(zip(lon,lat))

                #check for multipolygons
                if Key in NewPolygonDict:
                    Polygons = [Shape, NewPolygonDict[Key]]
                    Shape = cascaded_union(Polygons)

                #add converted polygon to new dictionary
                NewPolygonDict[Key] = Shape

    return NewPolygonDict

        def display(self, feature_name = None, *args):
            # allows function to be used by multiple processes, first option (when a feature_name is passed)
            # is for viewing data, second option is for viewing created geometries
            if feature_name:
                geom = self.datasets[self.current_dataset].features[feature_name][0]
                if geom.geom_type != ('Polygon' or 'MultiPolygon'):
                    self.dialog_text.set('This geometry is invalid. Please use a different dataset')
                    pass
                geom = cascaded_union(geom) #to dissolve multipolygons
                geom_bdry = geom.boundary
                minx, miny, maxx, maxy = self.bg.bounds
                w, h = maxx - minx, maxy - miny
                fig = plt.figure(1, figsize = (5, 5), dpi = 180, frameon = False)
                ax = fig.add_subplot(111)
                ax.set_xlim(minx,maxx)
                ax.set_ylim(miny,maxy)
                for poly in self.bg:
                    bg_patch = PolygonPatch(poly, fc = 'lightsage', ec = 'k', alpha = 1)
                    ax.add_patch(bg_patch)

                if geom.geom_type == 'MultiPolygon':
                    for poly in geom:
                        patch = PolygonPatch(poly, fc= 'teal', ec='navy', alpha = 0.5)
                        ax.add_patch(patch)
                else:
                    patch = PolygonPatch(geom, fc='teal', ec='navy', alpha = 0.5)
                    ax.add_patch(patch)
                plt.show()
            else:
                geom = args[0]
                name = args[1]
                geom = cascaded_union(geom) #to dissolve multipolygons
                minx, miny, maxx, maxy = self.bg.bounds
                w, h = maxx - minx, maxy - miny
                fig = plt.figure(1, figsize = (5, 5), dpi = 180, frameon = False)
                ax = fig.add_subplot(111)
                ax.set_xlim(minx,maxx)
                ax.set_ylim(miny,maxy)
                for poly in self.bg:
                    bg_patch = PolygonPatch(poly, fc = 'lightsage', ec = 'k', alpha = 1)
                    ax.add_patch(bg_patch)
                if geom.geom_type == 'MultiPolygon':
                    for poly in geom:
                        patch = PolygonPatch(poly, fc= 'teal', ec='navy', alpha = 0.5)
                        ax.add_patch(patch)
                else:
                    patch = PolygonPatch(geom, fc='teal', ec='navy', alpha = 0.5)
                    ax.add_patch(patch)
                plt.title(name)
                plt.show()

def test_split():
    print "Split donut by interior point"
    donut = sg.Point(0, 0).buffer(2.0).difference(sg.Point(0, 0).buffer(1.0))
    center = donut.centroid
    res = split_horiz_by_point(donut, center)
    print "Parts (should be 2): %s" % len(res)
    dif = cascaded_union(res).symmetric_difference(donut).area
    print "Difference: %s" % dif

    print "Split donut by exterior point"
    side = sg.Point(4,4)
    res = split_horiz_by_point(donut, side)
    print "Parts (should be 1): %s" % len(res)
    dif = cascaded_union(res).symmetric_difference(donut).area
    print "Difference: %s" % dif

	def generate_paths(self, shapes, tool_d):
		tool_r = tool_d / 2.0
		toolpath = []
		for shp in shapes:
			toolpath.append( shp.buffer(tool_r) )

		result = cascaded_union(toolpath)
		
		if result.geom_type == 'Polygon':
			opt_obj = result.simplify(tolerance=0.01, preserve_topology=False)
			tmp = self.add_holders(opt_obj.exterior, self.holder_size + tool_d, self.min_len)
			if type(tmp) == list:
				return tmp
			else:
				return [tmp]
			
		elif result.geom_type == 'MultiPolygon':
			tmp = []
			for f in result:
				opt_obj = f.simplify(tolerance=0.01, preserve_topology=False)
				tmp2 = self.add_holders(opt_obj.exterior, self.holder_size + tool_d, self.min_len)
				if type(tmp2) == list:
					for item in tmp2:
						tmp.append(item)
				else:
					tmp.append( tmp2 )
			return tmp

def alpha_shape(points, alpha):
    def add_edge(points, i, j):
        if (i, j) in edges or (j, i) in edges:
            return
        edges.add((i, j))
        edge_points.append((points[i], points[j]))
    tri = DelaunayTri(points)
    edges = set()
    edge_points = []
    for i1, i2, i3 in tri.vertices:
        x1, y1 = tri.points[i1]
        x2, y2 = tri.points[i2]
        x3, y3 = tri.points[i3]
        a = hypot(x1 - x2, y1 - y2)
        b = hypot(x2 - x3, y2 - y3)
        c = hypot(x3 - x1, y3 - y1)
        s = (a + b + c) / 2.0
        area = sqrt(s * (s - a) * (s - b) * (s - c))
        radius = a * b * c / (4 * area)
        if radius < 1.0 / alpha:
            add_edge(tri.points, i1, i2)
            add_edge(tri.points, i2, i3)
            add_edge(tri.points, i3, i1)
    shape = cascaded_union(list(polygonize(MultiLineString(edge_points))))
    return shape

def isGalleryCovered(camerasState):
    circles_array=[]
    for (i, camera) in enumerate(camerasState):
        if(camerasState[i]):
            circles_array.append(Gallery.cameras[i].circle)
    circles_union=cascaded_union(circles_array)
    return (circles_union.intersection(Gallery.gallery_polygon)==Gallery.gallery_polygon)

 def dissolve (self, inFile, outFile):
     # create dictionary for storing the uniqueRefs
     uniqueRefs = {}
     with fiona.open(inFile, 'r', encoding='utf-8') as input:
         input_driver = input.driver
         input_crs = input.crs
         input_schema = {'geometry': 'MultiLineString','properties': {'ref'.encode("utf-8"): 'str:254'}}
         with fiona.open(outFile, 'w', driver=input_driver, crs=input_crs, schema=input_schema, encoding='utf-8') as output:
             for item in input:
                 # extract the key, if the 'ref' attribute is NOT called 'ref' 
                 # you can insert the different attribute name HERE (and only HERE).
                 key = item['properties']['ids_and_re']
                 geom = shape(item['geometry'])
                 # find all motorways within the New Zealand mainland
                 # and remove all letters per  
                 newZeaBox = [(17920614.01, -4033681.682),(20362002, -4054837.565),(20357771.35, -6073108.484),(17683668.157,-6068877.308)]
                 newZeaPoly = Polygon(newZeaBox)
                 if geom.within(newZeaPoly):
                     key = re.sub(r'\D',"", key)
                 if not geom.type.startswith('Multi'):
                     geom = [geom]
                 for g in geom:
                     if key in uniqueRefs:
                         uniqueRefs[key].append(g)
                     else:
                         uniqueRefs[key] = [g]
             for key in uniqueRefs:
                 # omit lines that have blank 'ref' tags
                 if key is not None and key != 'None':
                     dissolve_feat = cascaded_union(uniqueRefs[key])
                     output.write({'geometry':mapping(dissolve_feat), 'properties': {'ref': key}})

def run_contour_paths(pseudo, epsilon, evals):
    # get paths
    paths = pseudo.contour_paths(epsilon)

    # check if pseudospectrum is correct by matching the parts of it
    import shapely.geometry as geom
    from shapely.ops import cascaded_union
    # create circles
    circles = [geom.Point(numpy.real(lamda), numpy.imag(lamda))
               .buffer(epsilon) for lamda in evals]
    exact_pseudo = cascaded_union(circles)
    exact_paths = pseudopy.utils.get_paths(exact_pseudo)

    N = len(paths)
    assert(N == len(exact_paths))

    # create polygons
    polys = [geom.Polygon([(numpy.real(z), numpy.imag(z))
                           for z in path.vertices])
             for path in paths]
    exact_polys = [geom.Polygon([(numpy.real(z), numpy.imag(z))
                                 for z in path.vertices])
                   for path in exact_paths]

    # match elements by measuring intersections
    M = numpy.zeros((N, N))
    for (i, j) in product(range(N), range(N)):
        M[i, j] = exact_polys[i].symmetric_difference(polys[j]).area
    for i in range(N):
        assert(numpy.min(M[i, :]) < 0.1*epsilon)