用于在较小多边形中细分多边形的算法

Question

我有一个由平面上的连续边构成的多边形,并且希望将其细分为三角形或矩形的子多边形.我在哪里可以找到算法来做到这一点？谢谢 !

Answer 1

在计算几何中,您要解决的问题称为三角测量.

有一些算法可以解决这个问题,给出了具有不同属性的三角测量.您需要确定哪一个最合适.

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)