获取位于Shapely多边形内部的所有点。

Question

我需要找出多边形内部和边上的所有格点。

输入:

from shapely.geometry import Polygon, mapping
sh_polygon = Polygon(((0,0), (2,0), (2,2), (0,2)))

输出：

(0, 0), (1, 0), (2, 0), (0, 1), (1, 1), (2, 1), (0, 2), (1, 2), (2, 2)

请建议是否有一种方法可以使用或不使用Shapely来获得预期的结果。

我编写了这段代码，它可以给出多边形内部的点，但不能给出在多边形上的点。同时，是否有更好的方法来完成同样的事情：

from shapely.geometry import Polygon, Point

def get_random_point_in_polygon(poly):
    (minx, miny, maxx, maxy) = poly.bounds
    minx = int(minx)
    miny = int(miny)
    maxx = int(maxx)
    maxy = int(maxy)
    print("poly.bounds:", poly.bounds)
    a = []
    for x in range(minx, maxx+1):
        for y in range(miny, maxy+1):
            p = Point(x, y)
            if poly.contains(p):
                a.append([x, y])
    return a


p = Polygon([(0,0), (2,0), (2,2), (0,2)])
point_in_poly = get_random_point_in_polygon(p)

print(len(point_in_poly))
print(point_in_poly)

输出：

poly.bounds: (0.0, 0.0, 2.0, 2.0)
1
[[1, 1]]

我已经简化了我的问题。实际上，我需要找到所有在一个正方形内部和边缘的点，该正方形的角落坐标为：(77,97), (141,101), (136,165), (73,160)。

- Beginner

---

你能否将最大值增加1，最小值减少1？ - SirParselot

@SirParselot 谢谢。我已经更改了它。但是，我得到的输出完全相同。 - Beginner

这是因为你有 if poly.contains(p)。只有当多边形内部有该点而不是在其上时，才会成立。如果循环在边界之间，您不必检查它是否在多边形内部。如果您删除 if poly.contains(p)，请忽略我的第一条评论。 - SirParselot

@SirParselot 我这么做是为了处理一个倾斜的多边形，例如 (77,97), (141,101), (136,165), (73,160)。你对于 (0,0), (2,0), (2,2), (0,2) 是正确的。 - Beginner

Answer 1

我会按照以下方式解决问题。

首先，定义一个点阵网格。例如，可以使用 itertools.product ：

from itertools import product
from shapely.geometry import MultiPoint

points = MultiPoint(list(product(range(5), repeat=2)))

points = MultiPoint(list(product(range(10), range(5))))

或者使用将x和y数组点的笛卡尔积转换为2D点的单个数组中的任何NumPy解决方案：

import numpy as np

x = np.linspace(0, 1, 5)
y = np.linspace(0, 1, 10)
points = MultiPoint(np.transpose([np.tile(x, len(y)), np.repeat(y, len(x))]))

然后，使用 Shapely 的 intersection 方法，我们可以得到那些既在给定多边形内部又在其边界上的点。

对于您提供的示例：

p = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
xmin, ymin, xmax, ymax = p.bounds
x = np.arange(np.floor(xmin), np.ceil(xmax) + 1)  # array([0., 1., 2.])
y = np.arange(np.floor(ymin), np.ceil(ymax) + 1)  # array([0., 1., 2.])
points = MultiPoint(np.transpose([np.tile(x, len(y)), np.repeat(y, len(x))]))
result = points.intersection(p)

而对于一个更复杂的例子：

p = Polygon([(-4.85571368308564, 37.1753007358263), 
             (-4.85520937147867, 37.174925051829), 
             (-4.85259349198842, 37.1783463712614),
             (-4.85258684662671, 37.1799609243756),
             (-4.85347524651836, 37.1804461589773),
             (-4.85343407576431, 37.182006629169),
             (-4.85516283166052, 37.1842384372115),
             (-4.85624511894443, 37.1837967179202),
             (-4.85533824429553, 37.1783762575331),
             (-4.85674599573635, 37.177038261295),
             (-4.85571368308564, 37.1753007358263)])
xmin, ymin, xmax, ymax = p.bounds  # -4.85674599573635, 37.174925051829, -4.85258684662671, 37.1842384372115
n = 1e3
x = np.arange(np.floor(xmin * n) / n, np.ceil(xmax * n) / n, 1 / n)  # array([-4.857, -4.856, -4.855, -4.854, -4.853])
y = np.arange(np.floor(ymin * n) / n, np.ceil(ymax * n) / n, 1 / n)  # array([37.174, 37.175, 37.176, 37.177, 37.178, 37.179, 37.18 , 37.181, 37.182, 37.183, 37.184, 37.185])
points = MultiPoint(np.transpose([np.tile(x, len(y)), np.repeat(y, len(x))]))
result = points.intersection(p)

- Georgy

---

非常感谢您的帮助，这真的帮助我理解了一个类似的任务，但是我的需求是绘制多边形而不是点。如果您有兴趣回答我的问题，这是我的问题链接：https://stackoverflow.com/questions/66253285/stacking-small-polygons-inside-another-bigger-one。再次感谢！ - salRad

Answer 2

有没有一个函数可以找到落在一条直线上的点？这些是你缺少的唯一点。它们只是线段定义方程的解。如果没有，自己编写算法通过暴力方法找到这些点也很容易。

对于多边形的每条边（p1，p2），执行以下操作。

p1 = (x1, y1)
p2 = (x2, y2)
xdiff = x2 - x1
ydiff = y2 - y1

# Find the line's equation, y = mx + b
m = ydiff / xdiff
b = y1 - m*x1

for xval in range(x1+1, x2):
    yval = m * xval + b
    if int(yval) == yval:
        # add (xval, yval) to your list of points

我已经留下了细节，让你来处理：确保x1 < x2（或进行其他适应），处理垂直线段等。这并不特别优雅，但它很快，易于实现和调试。

- Prune

---

那么，您是建议我使用您的方法来查找在多边形上的点，并使用我的代码来查找在多边形内部的点吗？ - Beginner

是的。或者，您可以使用我的算法找到多边形每个列x值的边界点，然后沿着每个列进行迭代。但是，既然您已经有了一个方便的调用以查找内部点，为什么要费力呢？ - Prune

如果您决定自己查找所有点，请注意需要知道多边形的内部方向，并在该方向上四舍五入计算出yval。这些yval形成每个点列的限制。 - Prune

非常感谢你。 - Beginner