29

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.

  • 3
    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
  • 3
    possible duplicate of Distribute points respect to features using Quantum Gis – Pablo Feb 22 '11 at 21:15
  • 1
    and also: gis.stackexchange.com/questions/4663/… – 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
20

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

  • 2
    +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
  • 1
    @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
14

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

enter image description here

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

  • 1
    Even 5 years later, still really helpful! – Stranded Kid Jul 12 '16 at 12:29
10

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

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
6

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

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.

2

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.

Generating random points in selected polygons

  • 1
    Obligatory reminder that zip codes are lines, not polygons. – Richard Mar 8 '18 at 17:25
  • Can you elaborate? To me also here it refers to areas: en.wikipedia.org/wiki/ZIP_Code#Primary_state_prefixes – markusN Mar 12 '18 at 22:49
  • Sure: zip codes reference particular post offices and their mail delivery routes. As a result, zip codes are lines, not polygons. They can overlap with each other, contain holes, and don't necessarily cover the entire US or any given state. Using them to divide up area is dangerous for this reason. Census units (like block groups) are a better choice. See also: this and this. – Richard Mar 12 '18 at 23:14
  • 1
    Thanks! It probably also depends on the country, see for example en.wikipedia.org/wiki/Postal_codes_in_Germany - however, ZIP codes are not my core topic, just wanted to illustrate and answer the original question "Generate points that lie inside polygon" rather than discuss ZIP code definitions which is OT here :-) – markusN Mar 13 '18 at 13:35
  • 1
    Noted on both counts. I should probably make a little blog entry so I can say this all more succinctly next time :-) – Richard Mar 13 '18 at 16:19
1

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====================================================================== enter image description here

  1. 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===================================================================enter image description here

  1. 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==========================================================================

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.