# Generate points that lie inside polygon

I have polygon feature and want to be able to generate points inside it. I need this for one classification task.

Generating random points until one is inside the polygon wouldn't work because it's really unpredictable the time it takes.

-
On the contrary, the time is predictable. It is proportional to the ratio of the polygon's extent area divided by the polygon's area, times the time needed to generate and test a single point. The time varies a little, but the variation is proportional to the square root of the number of points. For largish numbers, this becomes inconsequential. Break tortuous polygons into more compact pieces if necessary to get this area ratio down to a low value. Only fractal polygons will give you trouble, but I doubt you have them! – whuber Feb 22 '11 at 18:34
possible duplicate of Distribute points respect to features using Quantum Gis – Pablo Feb 22 '11 at 21:15
– Pablo Feb 22 '11 at 21:18
@Pablo: good finds. However, both of those questions are software specific and both concern placing regular arrays of points within polygons, not random points – whuber Feb 22 '11 at 21:29
agree with whuber difference is random points vs regular point generation within a polygon. – Mapperz Feb 22 '11 at 22:53

Start by decomposing the polygon into triangles, then generate points inside those. (For a uniform distribution, weight each triangle by its area.)

-
+1 Simple and effective. It's worth pointing out that uniformly random points can be generated within a triangle with no rejection at all, because there are (easily computed) area-preserving mappings between any triangle and an isosceles right triangle, which is half a square, say the half where the y coordinate exceeds the x coordinate. Generate two random coordinates and sort them to obtain a random point in the isosceles triangle, then map that back to the original triangle. – whuber Feb 22 '11 at 20:19
+1 I really like the discussion of trilinear coordinates referenced by the article you cite. I suppose this would be amenable to sphere whose surface is represented as a tesselation of triangles. On a projected plane, it wouldn't be a truly random distribution, would it? – Kirk Kuykendall Feb 22 '11 at 20:42
@whuber -- +1 back at you. Another way (in the link, but they hand-waved over it) is to reflect the rejected points from the uniformly-sampled quadrilateral across the shared edge and back into the triangle. – Dan S. Feb 22 '11 at 21:07
@Kirk -- the citation link is a touch anti-helpful in that it lists a bunch of wrong (non-uniform) sampling methods, including trilinear coordinates, before the "right" way. It doesn't look like there's a direct way to get a uniform sampling w/ trilinear coordinates. I'd approach uniform sampling over the entire sphere by converting random unit vectors in 3d to their lat/lon equivalent, but that's just me. (Unsure about sampling constrained to spherical triangles/polygons.) (Also unsure about truly uniform sampling on e.g. wgs84: just picking angles will bias a bit towards the poles, I think.) – Dan S. Feb 22 '11 at 21:15
@Dan For uniformly sampling the sphere, use a cylindrical equal-area projection (coordinates are longitude and cosine of latitude). If you want to sample without using a projection, there's a beautiful trick: generate three independent standard normal variates (x,y,z) and project them to the point (Rx/n, Ry/n, R*z/n) where n^2 = x^2 + y^2 + z^2 and R is the earth radius. Convert to (lat, lon) if need be (using authalic latitudes when working on a spheroid). It works because this trivariate normal distribution is spherically symmetric. For sampling triangles, stick to a projection. – whuber Feb 22 '11 at 21:26

As you put a QGIS tag on this question: Random Points tool can be used with a boundary layer.

If you are looking for code, the underlying plugin source code should be of help.

-
Even 5 years later, still really helpful! – Stranded Kid Jul 12 at 12:29

You could determine the extent of the polygon, then constrain the random number generation for X and Y values within those extents.

Basic process: 1) Determine maxx, maxy, minx, miny of polygon vertices, 2) Generate random points using these values as bounds 3) Test each point for intersection with your polygon, 4) Stop generating when you have enough points satisfying the intersection test

Here is an algorithm (C#) for the intersection test:

``````bool PointIsInGeometry(PointCollection points, MapPoint point)
{
int i;
int j = points.Count - 1;
bool output = false;

for (i = 0; i < points.Count; i++)
{
if (points[i].X < point.X && points[j].X >= point.X || points[j].X < point.X && points[i].X >= point.X)
{
if (points[i].Y + (point.X - points[i].X) / (points[j].X - points[i].X) * (points[j].Y - points[i].Y) < point.Y)
{
output = !output;
}
}
j = i;
}
return output;
}
``````
-

If R is an option, see `?spsample` in the `sp` package. The polygons can be read in from any GDAL-supported format built into the rgdal package, and then `spsample` works directly on imported object with a variety of sampling options.

-
+1 - Since R is open source if one wants to replicate you can always go into the source to see how they are done. For point patterns one may also be interested in the simulation tools in the spatstat package. – Andy W Feb 24 '11 at 14:05

I would like to offer a solution that requires very little in terms of GIS analysis. In particular, it does not require triangulating any polygons.

The following algorithm, given in pseudocode, refers to some simple operations in addition to basic list handling capabilities (create, find length, append, sort, extract sublists, and concatenate) and generation of random floats in the interval [0, 1):

``````Area:        Return the area of a polygon (0 for an empty polygon).
BoundingBox: Return the bounding box (extent) of a polygon.
Width:       Return the width of a rectangle.
Height:      Return the height of a rectangle.
Left:        Split a rectangle into two halves and return the left half.
Right:       ... returning the right half.
Top:         ... returning the top half.
Bottom:      ... returning the bottom half.
Clip:        Clip a polygon to a rectangle.
RandomPoint: Return a random point in a rectangle.
Search:      Search a sorted list for a target value.  Return the index
of the last element less than the target.
In:          Test whether a point is inside a polygon.
``````

These are all available in almost any GIS or graphics programming environment (and easy to code if not). `Clip` must not return degenerate polygons (that is, those with zero area).

Procedure `SimpleRandomSample` efficiently obtains a list of points randomly distributed within a polygon. It is a wrapper for `SRS`, which breaks the polygon into smaller pieces until each piece is sufficiently compact to be sampled efficiently. To do this, it uses a precomputed list of random numbers to decide how many points to allocate to each piece.

SRS can be "tuned" by changing the parameter `t`. This is the maximum bounding box:polygon area ratio that can be tolerated. Making it small (but greater than 1) will cause most polygons to be split into many pieces; making it large can cause many trial points to be rejected for some polygons (sinuous, with slivers, or full of holes). This guarantees that the maximum time to sample the original polygon is predictable.

``````Procedure SimpleRandomSample(P:Polygon, N:Integer) {
U = Sorted list of N independent uniform values between 0 and 1
Return SRS(P, BoundingBox(P), U)
}
``````

The next procedure calls itself recursively if necessary. The mysterious expression `t*N + 5*Sqrt(t*N)` conservatively estimates an upper limit on how many points will be needed, accounting for chance variability. The likelihood that this will fail is only 0.3 per million procedure calls. Increase 5 to 6 or even 7 to reduce this likelihood if you like.

``````Procedure SRS(P:Polygon, B:Rectangle, U:List) {
N = Length(U)
If (N == 0) {Return empty list}
aP = Area(P)
If (aP <= 0) {
Error("Cannot sample degenerate polygons.")
Return empty list
}
t = 2
If (aP*t < Area(B)) {
# Cut P into pieces
If (Width(B) > Height(B)) {
B1 = Left(B); B2 = Right(B)
} Else {
B1 = Bottom(B); B2 = Top(B)
}
P1 = Clip(P, B1); P2 = Clip(P, B2)
K = Search(U, Area(P1) / aP)
V = Concatenate( SRS(P1, B1, U[1::K]), SRS(P2, B2, U[K+1::N]) )
} Else {
# Sample P
V = empty list
maxIter = t*N + 5*Sqrt(t*N)
While(Length(V) < N and maxIter > 0) {
Decrement maxIter
Q = RandomPoint(B)
If (Q In P) {Append Q to V}
}
If (Length(V) < N) {
Error("Too many iterations.")
}
}
Return V
}
``````
-

There are some good libraries out there that does most of the heavy lifting for you.

Example using shapely in python.

``````import random
from shapely import Polygon, Point

def get_random_point_in_polygon(poly):
(minx, miny, maxx, maxy) = poly.bounds
while True:
p = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if poly.contains(p):
return p

p = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
point_in_poly = get_random_point_in_polygon(mypoly)
``````
-