Generate sample coordinates inside a Polygon

I have a Polygon named as poly. I attempted to randomly select 5 coordinate points that lies inside the polygon.

import numpy as np
from shapely.geometry import Polygon, Point

poly = Polygon([(141.4378366,-25.95915986), (165.4279876,-29.43400298), (163.1382942,-47.65345814), (133.1675418,-42.99807751)])

minx, miny, maxx, maxy = poly.bounds

longs = np.arange(minx, maxx, 0.002); lats = np.arange(miny, maxy, 0.002)
longs = np.tile(longs,3).ravel(); lats = np.repeat(lats,3).ravel()
coords = np.array([(x,y) for x,y in zip(longs,lats)])

points = [Point(xy) for xy in coords]

check = [xy.within(poly) for xy in points]
pointsInside = coords[check]

ranIdx = np.random.choice(len(pointsInside),5,replace=False)
result = pointsInside[ranIdx]

print result


I think my code is ineffective. Are there any ideas for a straight and elegant implementation?

• – YXD Nov 14 '14 at 14:32
• @MrE In my problem just selecting 5 points inside the polygon is sufficient and it does not require that they are exactly uniformly and randomly distributed. – Borys Nov 14 '14 at 14:36

Style Nitpicks

Put spaces after commas / Keep Lines within 120 Chars

poly = Polygon([(141.4378366,-25.95915986), (165.4279876,-29.43400298), (163.1382942,-47.65345814), (133.1675418,-42.99807751)])


becomes

poly = Polygon([(141.4378366, -25.95915986), (165.4279876, -29.43400298), (163.1382942, -47.65345814),
(133.1675418, -42.99807751)])


Use Pythonic Underscores

minx, miny, maxx, maxy = poly.bounds


becomes

min_x, min_y, max_x, max_y = poly.bounds


Eschew Semicolons (Personal Opinion)

You can use semicolons like this:

longs = np.arange(minx, maxx, 0.002); lats = np.arange(miny, maxy, 0.002)


But personally, I just wouldn't in Python. Eww.

An Alternative Algorithm

EDIT: Having read @MrE's and @MartinR's comments, I now propose this rejection sampling method. Although, this could miss frequently in a polygon with a large bounding box relative to its area; .e.g. an 8-point Christmas star with a small inner circle.

def random_points_within(poly, num_points):
min_x, min_y, max_x, max_y = poly.bounds

points = []

while len(points) < num_points:
random_point = Point([random.uniform(min_x, max_x), random.uniform(min_y, max_y)])
if (random_point.within(poly)):
points.append(random_point)

return points


Another alternative that never misses, but may not be well distributed.

This was my original idea, but it didn't look so good in the morning.

Firstly, work out how to select one point within a polygon.

def random_point_within(poly):
min_x, min_y, max_x, max_y = poly.bounds

x = random.uniform(min_x, max_x)
x_line = LineString([(x, min_y), (x, max_y)])
x_line_intercept_min, x_line_intercept_max = x_line.intersection(poly).xy[1].tolist()
y = random.uniform(x_line_intercept_min, x_line_intercept_max)

return Point([x, y])


Then simply call that in a list comprehension to generate however many points you desire.

points = [random_point_within(poly) for i in range(5)]
checks = [point.within(poly) for point in points]


My approach is to select x randomly within the polygon, then constrain y.

This approach doesn't require NumPy and always produces a point within the polygon.

• I liked your code besides the way you have used linestring ratherthan finding the way to directy capture the x and y positions... – Borys Nov 14 '14 at 13:58
• I might be mistaken, but I think that your random point selection is not uniform with respect to the polygon. E.g. for the triangle with vertices (0,0), (1,0), (0,1) the chance to get a point with x > 0.5 is 50%, but the part of the triangle with x > 0.5 is only 25% of the total triangle area. – Martin R Nov 14 '14 at 14:14
• The simplest method to get uniformly distributed samples is rejection sampling: 1. Find the bounding box of the polygon. 2. Generate a uniform sample in the bounding box 3. If sample is inside polygon, return sample, else go to 2. – YXD Nov 14 '14 at 14:35
• @MrE Could you explain me how to generate samples within the bounding box, not necessarily uniformly random though. – Borys Nov 14 '14 at 14:39
• I personally prefer seeing 80-100 chars as a limit. 120 as an absolute max for rare cases or if heavily nested. – flakes Nov 14 '14 at 17:34

Rejection sampling was proposed in comments on the other answer. The problem with rejection sampling is that the area of a polygon can be an arbitrarily small fraction of its bounding box, for example:

def ε_poly(ε):
"Return a polygon that occupies a fraction ε of its bounding box."
assert 0 < ε <= 1
return Polygon([(0, 0), (1, 0), (ε, ε), (0, 1)])


Rejection sampling will take on average $$\1\over ε\$$ attempts to generate one sample point inside ε_poly(ε), and $$\1\over ε\$$ can be arbitrarily large.

A more robust approach is as follows:

1. Triangulate the polygon and calculate the area of each triangle.

2. For each sample:

1. Pick the triangle $$\t\$$ containing the sample, using random selection weighted by the area of each triangle.

2. Pick a random point uniformly in the triangle, as follows:

1. Pick a random point $$\x, y\$$ uniformly in the unit square.

2. If $$\x + y > 1\$$, use the point $$\1-x, 1-y\$$ instead. The effect of this is to ensure that the point is chosen uniformly in the unit right triangle with vertices $$\(0, 0), (0, 1), (1, 0)\$$

3. Apply the appropriate affine transformation to transform the unit right triangle to the triangle $$\t\$$.

Here's one possible implementation, using shapely.ops.triangulate and shapely.affinity.affine_transform:

import random
from shapely.affinity import affine_transform
from shapely.geometry import Point, Polygon
from shapely.ops import triangulate

def random_points_in_polygon(polygon, k):
"Return list of k points chosen uniformly at random inside the polygon."
areas = []
transforms = []
for t in triangulate(polygon):
areas.append(t.area)
(x0, y0), (x1, y1), (x2, y2), _ = t.exterior.coords
transforms.append([x1 - x0, x2 - x0, y2 - y0, y1 - y0, x0, y0])
points = []
for transform in random.choices(transforms, weights=areas, k=k):
x, y = [random.random() for _ in range(2)]
if x + y > 1:
p = Point(1 - x, 1 - y)
else:
p = Point(x, y)
points.append(affine_transform(p, transform))
return points


This is fine if $$\k\$$ is small, as indicated by the OP. (If $$\k\$$ is large, then you would want to vectorize the construction of the points in NumPy.)

used this function with raster based 10m spacing

may be useful:

def random_points_within_poygon_and_10m_sbuffer(poly, num_points):

min_x, min_y, max_x, max_y = poly.bounds
min_lon = min_x
max_lon = max_x
min_lat = min_y
max_lat = max_y
# obtain utm bounds and dimentions
utmZone = int((np.floor((min_lon+ 180)/6)% 60) + 1)
fieldProj = Proj("+proj=utm +zone="+str(utmZone)+", +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
UTMx_min, UTMy_max = fieldProj(min_lon, max_lat)
UTMx_max, UTMy_max = fieldProj(max_lon, max_lat)
UTMx_min, UTMy_min = fieldProj(min_lon, min_lat)
UTMx_min = math.floor(UTMx_min*0.2)*5.0 + 2.5 # make it standard grid to match img
UTMy_max = math.ceil(UTMy_max*0.2)*5.0 - 2.5  # make it standard grid to match img

utm_npix_x = int(0.1*(math.ceil((UTMx_max - UTMx_min)*0.1)*10))
utm_npix_y = int(0.1*(math.ceil((UTMy_max - UTMy_min)*0.1)*10))
#set spacing raster
spacing_utm_grid = np.arange(utm_npix_y*utm_npix_x).reshape((utm_npix_y, utm_npix_x))
spacing_utm_grid[:,:] = 0

points = []
while len(points) < num_points:
pair_coordinates = [random.uniform(min_x, max_x), random.uniform(min_y, max_y)]
random_point = Point(pair_coordinates)
if (random_point.within(poly)):
utm_x, utm_y = fieldProj(pair_coordinates[0], pair_coordinates[1])
utm_x_pix = math.floor((utm_x - UTMx_min)/10)
utm_y_pix = math.floor((UTMy_max -utm_y)/10)
if(spacing_utm_grid[utm_y_pix, utm_x_pix]==0):
points.append(random_point)
spacing_utm_grid[utm_y_pix,utm_x_pix]=1

return points

• Welcome to Code Review! You have presented an alternative solution, but haven't reviewed the code. Please explain your reasoning (how your solution works and why it is better than the original) so that the author and other readers can learn from your thought process. Please read Why are alternative solutions not welcome? – Sᴀᴍ Onᴇᴌᴀ Aug 23 '18 at 20:22