python - Get all lattice points lying inside a Shapely polygon

Question

I need to find all the lattice points inside and on a polygon.

Input:

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

Output:

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

Please suggest if there is a way to get the expected result with or without using Shapely.

I have written this piece of code that gives points inside the polygon, but it doesn't give points on it. Also is there a better way to do the same thing:

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)

Output:

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

I have simplified my problem. Actually, I need to find all points inside and on a square with corners: (77,97), (141,101), (136,165), (73,160).

Answer 1

I would approach the problem as follows.

First, define a grid of lattice points. One could use, for example, 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))))

or any NumPy solution from Cartesian product of x and y array points into single array of 2D points:

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

Then, using intersection method of Shapely we can get those lattice points that lie both inside and on the boundary of the given polygon.

For your given example:

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)

And for a bit more sophisticated example:

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)