Generate points that lie inside polygon

Question

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.

Answer 1

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

Answer 2

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.

Answer 3

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;
}

Answer 4

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.

Answer 5

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
}