将多边形细分为较小多边形的算法

Question

我有一个由平面上的连续边组成的多边形，并且想将其细分为三角形或矩形的子多边形。 我在哪里可以找到一个算法来做到这一点？ 谢谢！

Answer 1

经过多次搜索偶然发现。

感谢@Aditya Chhabra 的提交，它工作得很好，但get_squares_from_rect 由于迭代剪辑，对于小边长来说非常慢。

如果我们将所有 LineStrings 合并到一个集合中，然后一步进行剪辑和多边形化，我们可以立即做到这一点，我在 this question 中找到了这一点。

以前边长 0.0001 (EPSG:4326) 需要 > 1 分钟，现在不需要时间。

from shapely.ops import unary_union, polygonize, linemerge
from shapely.geometry import LineString
import numpy as np


def get_squares_from_rect_faster(RectangularPolygon, side_length=0.0025):
    rect_coords = np.array(RectangularPolygon.boundary.coords.xy)
    y_list = rect_coords[1]
    x_list = rect_coords[0]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells = int(np.round(width / side_length))
    ycells = int(np.round(height / side_length))

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)

    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]

    lines = horizontal_splitters + vertical_splitters
    lines.append(RectangularPolygon.boundary)
    lines = unary_union(lines)
    lines = linemerge(lines)

    return list(polygonize(lines))

Answer 2

我自己也在寻找这个问题的答案，但找不到。尝试将几块拼接在一起，结果如下。

这不一定是最理想的例程，但它为我完成了工作。如果您想提高性能，请尝试使用代码。

算法的简要说明：

1. 使用原始几何本身的边界、凸包的边界和最小旋转矩形，导出所有可能的矩形。

2. 将所有矩形分成指定边长的小方块。

3. 删除重复项使用四舍五入的质心。 (r: 四舍五入参数)

4. 保留那些“在”几何体中的方格，或者那些“相交”几何体的方格，取决于哪个更接近所需方格的总数。

已编辑

新解决方案

#### Python script for dividing any shapely polygon into smaller equal sized polygons

import numpy as np
from shapely.ops import split
import geopandas
from shapely.geometry import MultiPolygon, Polygon

def rhombus(square):
    """
    Naively transform the square into a Rhombus at a 45 degree angle
    """
    coords = square.boundary.coords.xy
    xx = list(coords[0])
    yy = list(coords[1])
    radians = 1
    points = list(zip(xx, yy))
    Rhombus = Polygon(
        [
            points[0],
            points[1],
            points[3],
            ((2 * points[3][0]) - points[2][0], (2 * points[3][1]) - points[2][1]),
            points[4],
        ]
    )
    return Rhombus


def get_squares_from_rect(RectangularPolygon, side_length=0.0025):
    """
    Divide a Rectangle (Shapely Polygon) into squares of equal area.

    `side_length` : required side of square

    """
    rect_coords = np.array(RectangularPolygon.boundary.coords.xy)
    y_list = rect_coords[1]
    x_list = rect_coords[0]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells = int(np.round(width / side_length))
    ycells = int(np.round(height / side_length))

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = RectangularPolygon
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)

    return square_polygons


def split_polygon(G, side_length=0.025, shape="square", thresh=0.9):
    """
    Using a rectangular envelope around `G`, creates a mesh of squares of required length.
    
    Removes non-intersecting polygons. 
            

    Args:
    
    - `thresh` : Range - [0,1]

        This controls - the number of smaller polygons at the boundaries.
        
        A thresh == 1 will only create (or retain) smaller polygons that are 
        completely enclosed (area of intersection=area of smaller polygon) 
        by the original Geometry - `G`.
        
        A thresh == 0 will create (or retain) smaller polygons that 
        have a non-zero intersection (area of intersection>0) with the
        original geometry - `G` 

    - `side_length` : Range - (0,infinity)
        side_length must be such that the resultant geometries are smaller 
        than the original geometry - `G`, for a useful result.

        side_length should be >0 (non-zero positive)

    - `shape` : {square/rhombus}
        Desired shape of subset geometries. 


    """
    assert side_length>0, "side_length must be a float>0"
    Rectangle    = G.envelope
    squares      = get_squares_from_rect(Rectangle, side_length=side_length)
    SquareGeoDF  = geopandas.GeoDataFrame(squares).rename(columns={0: "geometry"})
    Geoms        = SquareGeoDF[SquareGeoDF.intersects(G)].geometry.values
    if shape == "rhombus":
        Geoms = [rhombus(g) for g in Geoms]
        geoms = [g for g in Geoms if ((g.intersection(G)).area / g.area) >= thresh]
    elif shape == "square":
        geoms = [g for g in Geoms if ((g.intersection(G)).area / g.area) >= thresh]
    return geoms

# Reading geometric data

geo_filepath = "/data/geojson/pc_14.geojson"
GeoDF        = geopandas.read_file(geo_filepath)

# Selecting random shapely-geometry

G = np.random.choice(GeoDF.geometry.values)


squares   = split_polygon(G,shape='square',thresh=0.5,side_length=0.025)
rhombuses = split_polygon(G,shape='rhombus',thresh=0.5,side_length=0.025)

以前的解决方案：

import numpy as np
import geopandas
from shapely.ops import split
from shapely.geometry import MultiPolygon, Polygon, Point, MultiPoint

def get_rect_from_geom(G, r=2):
    """
    Get rectangles from a geometry.
    r = rounding factor.
    small r ==> more rounding off ==> more rectangles
    """
    coordinate_arrays = G.exterior.coords.xy
    coordinates = list(
        zip(
            [np.round(c, r) for c in coordinate_arrays[0]],
            [np.round(c, r) for c in coordinate_arrays[1]],
        )
    )
    Rectangles = []
    for c1 in coordinates:
        Coords1 = [a for a in coordinates if a != c1]
        for c2 in Coords1:
            Coords2 = [b for b in Coords1 if b != c2]
            x1, y1 = c1[0], c1[1]
            x2, y2 = c2[0], c2[1]
            K1 = [k for k in Coords2 if k == (x1, y2)]
            K2 = [k for k in Coords2 if k == (x2, y1)]
            if (len(K1) > 0) & (len(K2) > 0):
                rect = [list(c1), list(K1[0]), list(c2), list(K2[0])]
                Rectangles.append(rect)
    return Rectangles


def get_squares_from_rect(rect, side_length=0.0025):
    """
    Divide a rectangle into equal area squares
    side_length = required side of square
    """
    y_list = [r[1] for r in rect]
    x_list = [r[0] for r in rect]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells, ycells = int(np.round(width / side_length)), int(
        np.round(height / side_length)
    )

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = Polygon(rect)
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)
    return [np.stack(SQPOLY.exterior.coords.xy, axis=1) for SQPOLY in square_polygons]


def round_centroid(g, r=10):
    """
    Get Centroids.
    Round off centroid coordinates to `r` decimal points.
    """
    C = g.centroid.coords.xy
    return (np.round(C[0][0], r), np.round(C[1][0], r))


def subdivide_polygon(g, side_length=0.0025, r=10):
    """
    1. Create all possible rectangles coordinates from the geometry, its minimum rotated rectangle, and its convex hull.
    2. Divide all rectangles into smaller squares.

    small r ==> more rounding off ==> fewer overlapping squares. (these are dropped as duplicates)
    large r ==> may lead to a few overlapping squares.
    """
    #   Number of squares required.
    num_squares_reqd = g.area // (side_length ** 2)

    # Some of these combinations can be dropped to improve performance.
    Rectangles = []
    Rectangles.extend(get_rect_from_geom(g))
    Rectangles.extend(get_rect_from_geom(g.minimum_rotated_rectangle))
    Rectangles.extend(get_rect_from_geom(g.convex_hull))

    Squares = []
    for r in range(len(Rectangles)):
        rect = Rectangles[r]
        Squares.extend(get_squares_from_rect(rect, side_length=side_length))
    SquarePolygons = [Polygon(square) for square in Squares]
    GDF = geopandas.GeoDataFrame(SquarePolygons).rename(columns={0: "geometry"})
    GDF.loc[:, "centroid"] = GDF.geometry.apply(round_centroid, r=r)
    GDF = GDF.drop_duplicates(subset=["centroid"])

    wgeoms = GDF[GDF.within(g)].geometry.values
    igeoms = GDF[GDF.intersects(g)].geometry.values

    w = abs(num_squares_reqd - len(wgeoms))
    i = abs(num_squares_reqd - len(igeoms))

    print(w, i)
    if w <= i:
        return wgeoms
    else:
        return igeoms



geoms = subdivide(g)

Answer 3

在computational geometry中，你要解决的问题叫做triangulation。

有一些算法可以解决这个问题，给出具有不同属性的三角剖分。您需要决定哪一个最合适。