Python ops.polygonize函数代码示例

def doPolygonize():
  blocks = polygonize(lines)
  writeBlocks(blocks, args[0] + '-blocks.geojson')

  blocks = polygonize(lines)
  bounds = Polygon([
    [minlng, minlat],
    [minlng, maxlat],
    [maxlng, maxlat],
    [maxlng, minlat],
    [minlng, minlat]
  ])
  # Geometry transform function based on pyproj.transform
  project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:3785'),
    pyproj.Proj(init='EPSG:4326'))
  print bounds
  print transform(project, bounds)

  print 'finding holes'
  for index, block in enumerate(blocks):
    if index % 1000 == 0:
      print "diff'd  %s" % (index)
    if not block.is_valid:
      print explain_validity(block)
      print transform(project, block)
    else:
      bounds = bounds.difference(block)
  print bounds

def create_union(in_ply1, in_ply2, result_geojson):
    """
    Create union polygon
    :param in_ply1: first input shapely polygon
    :param in_ply2: second input shapely polygon
    :param result_geojson: output geojson file including full file path
    :return: shapely MultiPolygon
    """
    # union the polygon outer linestrings together
    outer_bndry = in_ply1.boundary.union(in_ply2.boundary)

    # rebuild linestrings into polygons
    output_poly_list = polygonize(outer_bndry)

    out_geojson = dict(type='FeatureCollection', features=[])

    # generate geojson file output
    for (index_num, ply) in enumerate(output_poly_list):
        feature = dict(type='Feature', properties=dict(id=index_num))
        feature['geometry'] = ply.__geo_interface__
        out_geojson['features'].append(feature)

    # create geojson file on disk
    json.dump(out_geojson, open(result_geojson, 'w'))

    # create shapely MultiPolygon
    ply_list = []
    for fp in polygonize(outer_bndry):
        ply_list.append(fp)

    out_multi_ply = MultiPolygon(ply_list)

    return out_multi_ply

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 alpha_shape(points, alpha):
    from shapely.ops import cascaded_union, polygonize
    from scipy.spatial import Delaunay
    import numpy as np
    import math
    import shapely.geometry as geometry
    """
    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

    def add_edge(edges, edge_points, coords, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
                # already added
            return
        edges.add( (i, j) )
        edge_points.append(coords[ [i, j] ])
    coords = np.array([point.coords[0]
                       for point in points])
    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

    def get_boundingpolygon(self):
        """
            TODO: Implement ncell bbox

            order of lines creation (assumes 0,0 is x)
            -----3-----
            |         |
            4         2
            |         |
            x----1-----
        """

        if self._ndim == 2: # CGRID
            nx,ny = self._xarray.shape
            one = MultiLineString([((self._xarray[i][0],self._yarray[i][0]),(self._xarray[i+1][0],self._yarray[i+1][0])) for i in range(nx-1)])
            two = MultiLineString([((self._xarray[nx-1][j],self._yarray[nx-1][j]),(self._xarray[nx-1][j+1],self._yarray[nx-1][j+1])) for j in range(ny-1)])
            three = MultiLineString([((self._xarray[i][ny-1],self._yarray[i][ny-1]),(self._xarray[i-1][ny-1],self._yarray[i-1][ny-1])) for i in reversed(list(range(1,nx)))])
            four = MultiLineString([((self._xarray[0][j],self._yarray[0][j]),(self._xarray[0][j-1],self._yarray[0][j-1])) for j in reversed(list(range(1,ny)))])
            m = one.union(two).union(three).union(four)
        else: # RGRID
            nx,ny = self._xarray.shape[0], self._yarray.shape[0]
            one = LineString([(self._xarray[i], self._yarray[0]) for i in range(nx)])
            two = LineString([(self._xarray[-1], self._yarray[i]) for i in range(ny)])
            three = LineString([(self._xarray[i], self._yarray[-1]) for i in reversed(list(range(nx)))])
            four = LineString([(self._xarray[0], self._yarray[i]) for i in reversed(list(range(ny)))])
            m = MultiLineString([one,two,three,four])

        polygons = list(polygonize(m))
        # -- polygonize returns a list of polygons, including interior features, the largest in area "should" be the full feature
        assert len(polygons) > 0, "Could not determine a polygon"
        polygon = sorted(polygons, key=lambda x: x.area)[-1]
        return polygon

def alpha_shape(x, y, alpha):
    coords = np.c_[x, y]
    tri = scipy.spatial.Delaunay(coords)
    
    # edges = []
    edge_coords = []
    for ia, ib, ic in tri.vertices:
        ab = math.sqrt((coords[ia,0]-coords[ib,0])**2 + (coords[ia,1]-coords[ib,1])**2)
        bc = math.sqrt((coords[ib,0]-coords[ic,0])**2 + (coords[ib,1]-coords[ic,1])**2)
        ac = math.sqrt((coords[ia,0]-coords[ic,0])**2 + (coords[ia,1]-coords[ic,1])**2)
        
        semiperim = 0.5*(ab+bc+ac)
        area = math.sqrt(semiperim*(semiperim-ab)*(semiperim-bc)*(semiperim-ac))
        radius = 0.25*ab*bc*ac/area
        
        if radius < 1.0/alpha:
            # edges.append((tuple(coords[ia]), tuple(coords[ib])))
            # edges.append((tuple(coords[ib]), tuple(coords[ic])))
            # edges.append((tuple(coords[ia]), tuple(coords[ic])))
            edge_coords.append(coords[[ia, ib]])
            edge_coords.append(coords[[ib, ic]])
            edge_coords.append(coords[[ia, ic]])
            
    # edges = set(edges)
    m = MultiLineString(edge_coords)
    return cascaded_union(list(polygonize(m))), edge_coords

def difference(line1, line2, close_ends=False):
    """ Create polygons from two LineString objects.

    Parameters:
        line1 (LineString) : a line representing the initial condition.
        line2 (LineString) : a line representing the final condition.
        close_ends (bool) : option to close open line ends with vertical line segments.

    Returns:
        intersections (Point array) : the intersections between the LineString objects.
        polygons (Polygon array) : the polygons between the lines.
        signs (int array) : contains values of +1 or -1 to identify polygons as cut or fill.

    """
    if close_ends==True:
        line1, line2 = close(line1, line2)

    intersections = line1.intersection(line2)

    segs1 = cut_by_points([line1], intersections)
    segs2 = cut_by_points([line2], intersections)

    polygons = polygonize([segs1, segs2])

    signs = sign(linemerge(segs1), linemerge(segs2))

    # can't pass the polygonize generator to my class so convert the polygons into an array
    polygontxt = []
    areas = []
    for i, poly in enumerate(polygons):
        polygontxt.append(poly)
        areas.append(poly.area*signs[i])
    cutfill = pnd.Series(asarray(areas), name='area')

    return intersections, polygontxt, cutfill

    def _cont_to_polys(self, cont_lats, cont_lons):
        polys = []

        start_idx = 0
        splits = []
        while True:
            try:
                split_idx = cont_lats[start_idx:].index(99.99)
                splits.append(start_idx + split_idx)
                start_idx += split_idx + 1
            except ValueError:
                break

        splits = [ -1 ] + splits + [ len(cont_lats) + 1 ]
        poly_lats = [ cont_lats[splits[i] + 1:splits[i + 1]] for i in xrange(len(splits) - 1) ]
        poly_lons = [ cont_lons[splits[i] + 1:splits[i + 1]] for i in xrange(len(splits) - 1) ]

        # Intersect with the US boundary shape file.
        for plat, plon in zip(poly_lats, poly_lons):
            cont = LineString(zip(plon, plat))

            if plat[0] != plat[-1] or plon[0] != plon[-1]:
                # If the line is not a closed contour, then it intersects with the edge of the US. Extend
                #   the ends a little bit to make sure it's outside the edge.
                dln = np.diff(plon)
                dlt = np.diff(plat)

                pre = [ (plon[0] - 0.5 * dln[0], plat[0] - 0.5 * dlt[0]) ]
                post = [ (plon[-1] + 0.5 * dln[-1], plat[-1] + 0.5 * dlt[-1]) ]

                cont.coords = pre + list(cont.coords) + post 

            # polygonize() will split the country into two parts: one inside the outlook and one outside.
            #   Construct test_ln that is to the right of (inside) the contour and keep only the polygon
            #   that contains the line
            test_ln = cont.parallel_offset(0.05, 'right')

            for poly in polygonize(self._conus.boundary.union(cont)):
                if (poly.crosses(test_ln) or poly.contains(test_ln)) and self._conus.contains(poly.buffer(-0.01)):
                    polys.append(poly)

        # Sort the polygons by area so we intersect the big ones with the big ones first.
        polys.sort(key=lambda p: p.area, reverse=True)

        # If any polygons intersect, replace them with their intersection.
        intsct_polys = []
        while len(polys) > 0:
            intsct_poly = polys.pop(0)
            pops = []
            for idx, poly in enumerate(polys):
                if intsct_poly.intersects(poly):
                    intsct_poly = intsct_poly.intersection(poly)
                    pops.append(idx)

            for pop_idx in pops[::-1]:
                polys.pop(pop_idx)

            intsct_polys.append(intsct_poly)
        return intsct_polys

def procesaLineaInterna(featuresExternas, featuresInternas, featuresCentroide, featureDefn):
	#print "Procedemos a procesar las lineas internas"
	
	centroides = []
	for centroide in featuresCentroide:
		#obtenemos la altura y el rotulo del estilo de cada centroide
		for n in centroide.GetStyleString().split(','):
			if n.startswith('s'):
				altura = float(n.replace('s:', '').replace('g', ''))
			elif n.startswith('t'):
				rotulo = n.split('"')[1]
		punto = centroide.GetGeometryRef()
		x = punto.GetX()
		y = punto.GetY()
		longitudRotulo = len(rotulo)
		factor = 0.15 * (altura * 3.3333)
		desfaseX = longitudRotulo * factor - 0.05
		punto.SetPoint(point = 0, x = x + desfaseX, y = y - 0.20)
		
		centroides.append((rotulo, punto))
	
	featuresProceso = featuresExternas + featuresInternas
	
	outFeature = []
	if len(featuresProceso) > 1:
		geometry_out = None
		for inFeature in featuresProceso:
			geometry_in = inFeature.GetGeometryRef()
			if geometry_out is None:
				geometry_out = geometry_in
				geometry_out = ogr.ForceToMultiLineString(geometry_out)
			else:
				geometry_out = geometry_out.Union(geometry_in) 
		
		lineasInternasShapely = loads(geometry_out.ExportToWkt())
		polygonsShapely = polygonize(lineasInternasShapely)
	
		polygonGeom = []
		for polygon in polygonsShapely:
			polygonGeom.append(ogr.CreateGeometryFromWkt(dumps(polygon)))
		
		for pol in polygonGeom:
			for cen in centroides:
				if pol.Contains(cen[1]):
					feature = ogr.Feature(featureDefn)
					feature.SetGeometry(pol)
					feature.SetField('rotulo', cen[0])
					outFeature.append(feature.Clone())
					feature.Destroy()
	else:
		feature = ogr.Feature(featureDefn)
		geometryPoly = ogr.BuildPolygonFromEdges(ogr.ForceToMultiLineString(featuresProceso[0].GetGeometryRef()), dfTolerance = 0)
		feature.SetGeometry(geometryPoly)
		feature.SetField('rotulo', centroides[0][0])
		outFeature.append(feature.Clone())
		feature.Destroy()
	
	
	return outFeature

def get_convex_hull(points):
    """ Find the convex hull of points, return as Shapely polygon """
    log.info("Finding convex hull of %i points", len(points))
    hull = ConvexHull(points)
    log.info("Found convex hull with %i facets", len(hull.simplices))
    hull_edges = []
    for simplex in hull.simplices:
        hull_edges.append(
            zip(points[simplex, 0], points[simplex, 1]))
    polygon = list(polygonize(hull_edges))
    assert(len(polygon) == 1)
    return polygon[0]

def get_shapes(filelike):
    shapes = []
    node_cache = {}
    way_cache = {}
    for thing in iter_osm_file(filelike):
        if type(thing) == Node:
            pt = (thing.lon, thing.lat)
            shape = Point(pt)

            node_cache[thing.id] = pt

            if thing.tags:
                shapes.append((thing, shape))

        elif type(thing) == Way:
            points = []
            for nd in thing.nds:
                node_loc = node_cache.get(nd)
                if node_loc:
                    points.append(node_loc)
                else:
                    raise Exception("Way %s references node %s which is not parsed yet." % (thing.id, nd))

            if way_is_polygon(thing):
                shape = Polygon(points)
            else:
                shape = LineString(points)

            way_cache[thing.id] = points

            if any(thing.tags):
                # Only include tagged things at this point. Otherwise,
                # the shapes that are part of multipolygon relations
                # will be included twice.
                shapes.append((thing, shape))

        elif type(thing) == Relation:
            if any([t.key == 'type' and t.value == 'multipolygon' for t in thing.tags]):
                parts = []
                for member in thing.members:
                    if member.type == 'way':
                        shape = way_cache.get(member.ref)
                        if not shape:
                            raise Exception("Relation %s references way %s which is not parsed yet." % (thing.id, member.ref))

                        parts.append(shape)

                # Polygonize will return all the polygons created, so the
                # inner parts of the multipolygons will be returned twice
                # we only want the first one
                shapes.append((thing, next(polygonize(parts))))

    return shapes

 def test_polygonize(self):
     lines = [
         LineString(((0, 0), (1, 1))),
         LineString(((0, 0), (0, 1))),
         LineString(((0, 1), (1, 1))),
         LineString(((1, 1), (1, 0))),
         LineString(((1, 0), (0, 0))),
         LineString(((5, 5), (6, 6))),
         Point(0, 0),
         ]
     result = list(polygonize(lines))
     self.assertTrue(all([isinstance(x, Polygon) for x in result]))

def save(datasource, filename):
    """ Save a Datasource instance to a named OGR datasource.
    """
    ext = splitext(filename)[1]
    
    out_driver = ogr.GetDriverByName(drivers.get(ext))
    out_source = out_driver.CreateDataSource(filename)
    
    if out_source is None:
        raise Exception('Failed creation of %s - is there one already?' % filename)
    
    out_layer = out_source.CreateLayer('default', datasource.srs, ogr.wkbMultiPolygon)
    
    for field in datasource.fields:
        field_defn = ogr.FieldDefn(field.name, field.type)
        field_defn.SetWidth(field.width)
        out_layer.CreateField(field_defn)
    
    for i in datasource._indexes():
    
        segments = datasource.db.execute("""SELECT x1, y1, x2, y2
                                            FROM segments
                                            WHERE (src1_id = ? OR src2_id = ?)
                                              AND removed = 0""", (i, i))
    
        lines = [datasource.memo_line(x1, y1, x2, y2) for (x1, y1, x2, y2) in segments]
        
        try:
            poly = polygonize(lines).next()
    
        except StopIteration:
            lost_area = datasource.shapes[i].area
            lost_portion = lost_area / (datasource.tolerance ** 2)
            
            if lost_portion < 4:
                # It's just small.
                print >> stderr, 'Skipped small feature #%(i)d' % locals()
                continue
    
            # This is a bug we don't understand yet.
            raise Exception('Failed to get a meaningful polygon out of large feature #%(i)d' % locals())
        
        feat = ogr.Feature(out_layer.GetLayerDefn())
        
        for (j, field) in enumerate(datasource.fields):
            feat.SetField(field.name, datasource.values[i][j])
        
        geom = ogr.CreateGeometryFromWkb(dumps(poly))
        
        feat.SetGeometry(geom)
    
        out_layer.CreateFeature(feat)

def intersect(a, b):
    ap, bp = [orient_offset_polygon(p) for p in [a,b]]
    overlap = ap.intersection(bp)

    if overlap.is_empty:
        return []

    if not hasattr(overlap, 'exterior'):
        return [Poly(p.exterior.coords, p.centroid.coords[0], 0)
                for p in polygonize(overlap)]
    
    return [Poly(overlap.exterior.coords,
                 overlap.centroid.coords[0], 0)]

 def processAlgorithm(self, progress):
     try:
         from shapely.ops import polygonize
         from shapely.geometry import Point, MultiLineString
     except ImportError:
         raise GeoAlgorithmExecutionException(
                 'Polygonize algorithm requires shapely module!')
     vlayer = dataobjects.getObjectFromUri(
             self.getParameterValue(self.INPUT))
     output = self.getOutputFromName(self.OUTPUT)
     vprovider = vlayer.dataProvider()
     if self.getParameterValue(self.FIELDS):
         fields = vprovider.fields()
     else:
         fields = QgsFields()
     if self.getParameterValue(self.GEOMETRY):
         fieldsCount = fields.count()
         fields.append(QgsField('area', QVariant.Double, 'double', 16, 2))
         fields.append(QgsField('perimeter', QVariant.Double,
                                'double', 16, 2))
     allLinesList = []
     features = vector.features(vlayer)
     current = 0
     total = 40.0 / float(len(features))
     for inFeat in features:
         inGeom = inFeat.geometry()
         if inGeom.isMultipart():
             allLinesList.extend(inGeom.asMultiPolyline())
         else:
             allLinesList.append(inGeom.asPolyline())
         current += 1
         progress.setPercentage(int(current * total))
     progress.setPercentage(40)
     allLines = MultiLineString(allLinesList)
     allLines = allLines.union(Point(0, 0))
     progress.setPercentage(45)
     polygons = list(polygonize([allLines]))
     progress.setPercentage(50)
     writer = output.getVectorWriter(fields, QGis.WKBPolygon, vlayer.crs())
     outFeat = QgsFeature()
     current = 0
     total = 50.0 / float(len(polygons))
     for polygon in polygons:
         outFeat.setGeometry(QgsGeometry.fromWkt(polygon.wkt))
         if self.getParameterValue(self.GEOMETRY):
             outFeat.setAttributes([None] * fieldsCount + [polygon.area,
                                   polygon.length])
         writer.addFeature(outFeat)
         current += 1
         progress.setPercentage(50 + int(current * total))
     del writer