确定分割shapely几何图形的“左”和“右”侧

Question

我的问题是：如何确定一个已经分割旋转的矩形几何的Aside和Bside哪个是任意LineString的“左”和“右”侧，该LineString将该几何体分割？

对于这个问题，“左”和“右”的定义是指从节点到节点“行走”时LineString分裂器的左手和右手侧。

我已经创建了这个函数，可以将任何任意shapely几何（非集合）分成两个部分-“左”和“右”。

import shapely.geometry as geo
import shapely.ops as ops

def splitLR(geom, splitter):
    """Split a geometry into a 'left' and 'right' side using the shapely API"""
    if not isinstance(splitter, geo.LineString):
        raise TypeError("The splitter must be a LineString")
    if not splitter.is_simple:
        raise ValueError("Only simple splitter objects allowed")
    if hasattr(geom, "__iter__"):
        raise ValueError("Geometry collections not allowed")
    geom_extents = geo.GeometryCollection([geom, splitter]).minimum_rotated_rectangle
    sides = ops.split(geom_extents, splitter)
    try:
        Aside, Bside = sides
    except TypeError:
        # only 1 result - rotated rectangle wasn't split
        if len(ops.split(geom,splitter)) == 1:
            # geom isn't split by splitter
            raise ValueError("the splitter does not appear to split the geometry")
        else:
            # splitter too small for algorithm
            raise ValueError("the splitter must extend beyond minimum_rotated_rectangle "
                             "of the combined geometry")
    # determine which is Lside and Rside here
    Lside,Rside = get_LRsides(Aside, Bside, splitter)
    return tuple(side.intersection(geom) for side in (Lside, Rside))

以上的想法在这里链接的笔记本中有所说明（与上面相同的链接）。

http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/ll3mchnx0jwzjnf/determine%20left-right%20split.ipynb

总结一下：A面和B面是包围geom几何图形和splitter线串的minimum_rotated_rectangle的两个面。当执行side.intersection(geom)时，结果是原始给定几何图形geom在该面中包含的部分。

注意事项：

• 对于奇怪形状的“土豆”类型物体，这种交集可能导致一个或两个方向上出现多个对象（请参见nbviewer示例）
• 我在这里创建了自己的函数（而不是使用ops.split），因为ops.split函数只返回一个“袋子”中的拆分对象，并且没有办法确定它们在哪一侧（据我所知）

目前，我对get_LRsides的调用仅执行此函数，显然是无用的：

def get_LRsides(Aside, Bside, splitter):
    """Determine the 'left' and 'right' sides of an already split geometry"""
    return Aside,Bside

如何成功地将A和B标记为“左”和“右”？

- Rick

Answer 1

这可能行得通：

1. 用分裂器的端点和Aside中的一个点形成一个LinearRing
2. 应用object.is_ccw
3. 如果它返回True，则Aside在分裂器的左侧。

https://shapely.readthedocs.io/en/stable/manual.html#object.is_ccw

- yang5

Answer 2

我不是工具箱中最聪明的GIS工具，一直在努力想弄清楚如何将“Aside中的一个点”添加到已接受答案中的LinearRing中。这是我能想出来的替代方法。

我将相交分隔线串的两个点投影到第一个多边形的外部，从而为我提供了沿着该外部线串的距离。如果该距离增加，则第一个多边形是分割的“正确”一侧。

以下是示例代码：

from shapely.ops import split
from shapely.geometry import Polygon, LineString, Point


def rhs_split(poly, splitter):
    intersect_splitter = splitter.intersection(poly)
    geomcollect = split(poly, splitter)
    polya, polyb = geomcollect.geoms[0], geomcollect.geoms[1]
    # Test directionality
    pt0 = Point(intersect_splitter.coords[0])
    pt1 = Point(intersect_splitter.coords[1])
    start_dist = polya.exterior.project(pt0)
    end_dist = polya.exterior.project(pt1)
    return polya if end_dist > start_dist else polyb


print(
    rhs_split(
        Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
        LineString([(0.5, -0.1), (0.5, 1.1)]),
    )
)
# POLYGON ((0.5 1, 1 1, 1 0, 0.5 0, 0.5 1))
print(
    rhs_split(
        Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
        LineString([(0.5, 1.1), (0.5, -0.1)]),
    )
)
# POLYGON ((0 0, 0 1, 0.5 1, 0.5 0, 0 0))

- daryl