Python shapely.ops.polygonize() Examples

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 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 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 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 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 find_fit_corners(tile1, tile2, left):
    # Draw missed corners polygons
    c = can_be_corner(tile1, tile2, left)
    ret = None
    if c is not None:
        x1 = c[0][0]
        x2 = c[1][0]
        y1 = c[0][1]
        y2 = c[1][1]
        pol = fit_third_order(x1, x2, y1, y2, c[0][2], c[1][2])
        if pol is not None:
            pol_points = []
            if left:
                pol_points.append((min(y1,y2), min(x1, x2)))
            else:
                pol_points.append((min(y1,y2), max(x1, x2)))
            pol_points_a = []
            for x_i in range(min(x1, x2), max(x1, x2)+1):
                y_i = int(poly_3(*pol, x_i))
                pol_points_a.append((int(y_i), int(x_i)))
            if all([0<=x[0]<=img_side_len for x in pol_points_a])\
                    and all([0<=x[1]<=img_side_len for x in pol_points_a]):
                pol_points.extend(pol_points_a)
            if left:
                pol_points.append((max(y1, y2), max(x1, x2)))
                pol_points.append((img_side_len, 0))
            else:
                pol_points.append((max(y1,y2), min(x1, x2)))
                pol_points.append((img_side_len, img_side_len))
            
            square_points = [(0, 0), (img_side_len, 0),
                             (img_side_len, img_side_len),
                             (0, img_side_len)]
            square_poly = Polygon(square_points)
            line = LineString(pol_points)

            merged = linemerge([square_poly.boundary, line])
            borders = unary_union(merged)
            polygons = []
            for poly in polygonize(borders):
                polygons.append(poly)
            if len(polygons) > 1:
                return mask_for_polygons([polygons[1]],
                                         (img_side_len, img_side_len)) 

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

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

    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 alpha_shape(self, coords, alpha):
        import shapely.geometry as geometry
        from shapely.ops import cascaded_union, polygonize
        from scipy.spatial import Delaunay
        """
        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] ])
        
        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 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)