Vjylku

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

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.

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;

}

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

Example using shapely 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)

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.

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

}

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.

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 link

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

Three algorithms using different approaches.

Git Repo Link

1. Here is a simple and best approach, using the actual distance of coordinates from 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 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 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==========================================================================