Generating points that lie inside polygon using QGIS

Question

I am using QGIS.

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

The 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

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

Example using [shapely][1] in python.

import random
from shapely.geometry 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)

Or use .representative_point() to get a point within the object (as mentioned by dain):

Returns a cheaply computed point that is guaranteed to be within the geometric object.

poly.representative_point().wkt
'POINT (-1.5000000000000000 0.0000000000000000)'

  [1]: https://shapely.readthedocs.io

Answer 4

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 5

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 6

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
}

Answer 7

If your polygon is convex and you know all the vertices, you might want to consider doing a "random" convex weighting of the vertices to sample a new point which is guaranteed to lie inside the convex hull (polygon in this case).

For example say you have a N sided convex polygon with vertices

V_i, i={1,..,N}

Then generate randomly N convex weights

 w_1,w_2,..,w_N such that ∑ w_i = 1; w_i>=0

The randomly sampled point is then given by

Y= ∑ w_i*V_i

There can be different way to sample N convex weights

• Pick N-1 numbers uniformly randomly within a range (without replacement) , sort them and normalize the N intervals between them to get the weights.
• You can also sample from the Dirichlet distribution which is often used as a conjugate prior for the multinomial distribution, which is similar to the convex weights in your case.

When your polygon is not very severely non-convex you might consider first converting it to a convex hull. This should at least limit the number of points lying outside your polygon to a large extent.

Answer 8

The task is very easy to solve in GRASS GIS (one command) using v.random.

Below an example on how to add 3 random points into selected polygons (here ZIP code areas of the city of Raleigh, NC) from the manual page. By modifying the SQL "where" statement the polygon(s) can be selected.

Answer 9

Answer link

https://gis.stackexchange.com/a/307204/103524

Three algorithms using different approaches.

Git Repo: https://github.com/imran-5/Postgis-Custom

1. Here is a simple and best approach, using the actual distance of coordinates from the x and y direction. The internal algorithm use the WGS 1984 (4326) and result transform to inserted SRID.

Function ===================================================================

CREATE OR REPLACE FUNCTION public.I_Grid_Point_Distance(geom public.geometry, x_side decimal, y_side decimal)
RETURNS public.geometry AS $BODY$
DECLARE
x_min decimal;
x_max decimal;
y_max decimal;
x decimal;
y decimal;
returnGeom public.geometry[];
i integer := -1;
srid integer := 4326;
input_srid integer;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
    geom := ST_SetSRID(geom, srid);
        ----RAISE NOTICE 'No SRID Found.';
    ELSE
        ----RAISE NOTICE 'SRID Found.';
END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_min := ST_XMin(geom);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    y := ST_YMin(geom);
    x := x_min;
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
<<yloop>>
LOOP
IF (y > y_max) THEN
    EXIT;
END IF;

CASE i WHEN 0 THEN 
    y := ST_Y(returnGeom[0]);
ELSE 
    y := ST_Y(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), y_side, radians(0))::geometry);
END CASE;

x := x_min;
<<xloop>>
LOOP
  IF (x > x_max) THEN
      EXIT;
  END IF;
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
    x := ST_X(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), x_side, radians(90))::geometry);
END LOOP xloop;
END LOOP yloop;
RETURN
ST_CollectionExtract(st_transform(ST_Intersection(st_collect(returnGeom), geom), input_srid), 1);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE;

Use the function with a simple query, geometry must be valid and polygon, multi-polygons, or envelope

SELECT I_Grid_Point_Distance(geom, 50, 61) from polygons limit 1;

Result======================================================================

2. Second function based on Nicklas Avén algorithm. Have tried to handle any SRID.

   I have applied the following changes in the algorithm.

   1. Separate variable for x and y direction for pixel size,
   2. New variable for calculates the distance in spheroid or ellipsoid.
   3. Input any SRID, function transform Geom to the working environment of Spheroid or Ellipsoid Datum, then apply the distance to each side, get the result and transform to input SRID.

Function ===================================================================

CREATE OR REPLACE FUNCTION I_Grid_Point(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$ 
DECLARE
x_max decimal; 
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer; 
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
        srid := 4326;
        x_side := x_side / 100000;
        y_side := y_side / 100000;
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);
RETURN QUERY
WITH res as (SELECT ST_SetSRID(ST_MakePoint(x, y), srid) point FROM
generate_series(x_min, x_max, x_side) as x,
generate_series(y_min, y_max, y_side) as y
WHERE st_intersects(geom, ST_SetSRID(ST_MakePoint(x, y), srid))
) select ST_TRANSFORM(ST_COLLECT(point), input_srid) from res;
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Use it with a simple query.

SELECT I_Grid_Point(geom, 22, 15, false) from polygons;

Result===================================================================

3. Function based on the series generator.

Function==================================================================

CREATE OR REPLACE FUNCTION I_Grid_Point_Series(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer;
x_series DECIMAL;
y_series DECIMAL;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);

    x_series := CEIL ( @( x_max - x_min ) / x_side);
    y_series := CEIL ( @( y_max - y_min ) / y_side );
RETURN QUERY
SELECT st_collect(st_setsrid(ST_MakePoint(x * x_side + x_min, y*y_side + y_min), srid)) FROM
generate_series(0, x_series) as x,
generate_series(0, y_series) as y
WHERE st_intersects(st_setsrid(ST_MakePoint(x*x_side + x_min, y*y_side + y_min), srid), geom);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Use it with a simple query.

SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons; Result==========================================================================

Answer 10

Many of the answers already given point in the right direction, but the answer seems fairly simple: generate an array of regularly spaced (spacing equal to your required precision) points within the polygon extent then select those points that intersect the polygon. Use the selection of regularly spaced points to sample from. Select the number of these points that you require at random. The only issue is that planar projections of elipsoidal objects (e.g the Earth) will cause points farther from the equator to be slightly closer to each other in the x dimension than those closer to the equator. If the slight convergence away from the equator is an issue, one can calculate an array of points where the spacing is different for each row and diverges to compensate as the rows of points move away from the equator.