Извлечь все координаты внутри многоугольника из корректного многоугольника Python

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

Answer 1

Разве нет функции, которая найдет точки решетки, лежащие на прямой? Это единственные, кого ты пропустил. Они являются просто решениями определяющего уравнения отрезка. Если нет, то достаточно просто написать алгоритм самостоятельно, находя точки грубой силой.

Сделайте следующее для каждого ребра (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

Answer 2

Я бы подошел к проблеме следующим образом.

Сначала определите сетку из точек решетки. Можно использовать, например, 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))))

или любое решение NumPy из декартова произведения точек массива x и y в один массив 2D точек:

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

Затем, используя intersection Методом Shapely мы можем получить те точки решетки, которые лежат как внутри, так и на границе данного многоугольника.

Для вашего примера:

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)