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

Could you just increase your max's by 1 and decrease your min's by 1? — SirParselot, Jun 06 '17 at 21:49
@SirParselot Thansk. I changed it. But, I am getting the exact same output. — Beginner, Jun 06 '17 at 21:52
It's because you have `if poly.contains(p)`. That is only true if the polygon has the point inside of it and not on it. You shouldn't have to check if it's inside the polygon if your loops are between the bounds. Ignore my first comment if you remove `if poly.contains(p)` — SirParselot, Jun 06 '17 at 21:54
@SirParselot I did that to deal with a polygon which is tilted, like (77,97), (141,101), (136,165), (73,160). You are correct for (0,0), (2,0), (2,2), (0,2). — Beginner, Jun 06 '17 at 22:00

score 5 · Answer 1 · answered Nov 14 '19 at 14:00

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)

Thank you so much for this, it really helped me to understand a similar task I have, but instead, I need to plot polygons not points. Here is my question if you are interested in answering: https://stackoverflow.com/questions/66253285/stacking-small-polygons-inside-another-bigger-one. Thanks again — salRad, Feb 19 '21 at 03:26

score 2 · Accepted Answer · answered Jun 06 '17 at 21:52

2

Is there not a function that will find lattice points that lie on a line? Those are the only ones you're missing. They are simply solutions to the line segment's defining equation. If not, it's easy enough to write the algorithm yourself, finding the points by brute force.

Do the following for each edge (p1, p2) of the polygon.

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

I've left details up to you: make sure that x1 < x2 (or adapt otherwise), handle a vertical segment, etc. This isn't particularly elegant, but it's fast, easy to implement, and easy to debug.

answered Jun 06 '17 at 21:52

Prune

76,765
14
60
81

So, are you suggesting me to use your method to find points that are on the polygon and use my code to find points which are inside the polygon? – Beginner Jun 06 '17 at 22:02
Yes. Alternately, you can use my algorithm to find the boundary points for each column x-value) of the polygon, and then simply iterate along each column. However, since you already have a convenient call to find the interior points, why bother with that? – Prune Jun 06 '17 at 22:04
If you *do* decide to find all the points yourself, note that you need to know which direction is the interior of the polygon, and round the computed yval in that direction. Those yval's form the limits of each column of points. – Prune Jun 06 '17 at 22:05
Thanks you so much. – Beginner Jun 06 '17 at 22:10

Get all lattice points lying inside a Shapely polygon

2 Answers2

Linked