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?

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:

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

-
For each sample:

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

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

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

-
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)\$

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

Concern was expressed in comments that the affine transformation would cause the random points to get closer along one axis than the other, due to different scaling on the two axes. But that's not a good way to think about it, because closest points are not preserved by affine transforms. (A point might have a different nearest neighbour after the transform.) Instead, remember that affine transformations are linear in both axes, so that the chance of a random point appearing in a region remains proportional to the area of the region after the transform. This means that if the points were uniformly distributed in the unit right triangle, they remain uniformly distributed in the transformed triangle. Here's a quick demo:

import matplotlib.pyplot as plt
triangle = Polygon([(0, 0), (5, 0), (0, 1)])
points = random_points_in_polygon(triangle, 500)
plt.gca().set_aspect('equal')
plt.plot(*np.array(triangle.boundary.coords).T, 'g')
coords = np.array([p.coords for p in points]).reshape(-1, 2)
plt.plot(*coords.T, 'b.')
plt.show()

If the figure doesn't convince you, then here's some code that prints the mean absolute offsets along the two axes between the random points and their nearest neighbours:

import scipy.spatial
tree = scipy.spatial.KDTree(coords)
_, neighbours = tree.query(coords, [2])
print(np.mean(np.abs(coords - coords[neighbours[:,0]]), axis=0))

Typical output shows that the means on the two axes are very close even though the \$x\$-coordinates have been scaled and the \$y\$-coordinates have not:

[0.02201801 0.02486154]

StackExchange Code Review Q#69833, answer score: 11