Functions to Generate Random Points

Random Point in Polygon

Use this function when you need to create a random point inside a polygon. The function takes geometry as a required argument and max number of iterations as an optional argument (defaults to 1000 iterations). The function will raise an exception if max number of iterations has been exceeded. The function works by throwing a random point within the polygon bounding rectangle and then checking if the point is inside the polygon.

CREATE OR REPLACE FUNCTION RandomPoint (
                geom Geometry,
                maxiter INTEGER DEFAULT 1000
        )
        RETURNS Geometry
        AS $$
DECLARE
        i INTEGER := 0;
        x0 DOUBLE PRECISION;
        dx DOUBLE PRECISION;
        y0 DOUBLE PRECISION;
        dy DOUBLE PRECISION;
        xp DOUBLE PRECISION;
        yp DOUBLE PRECISION;
        rpoint Geometry;
BEGIN
        -- find envelope
        x0 = ST_XMin(geom);
        dx = (ST_XMax(geom) - x0);
        y0 = ST_YMin(geom);
        dy = (ST_YMax(geom) - y0);
        
        WHILE i < maxiter LOOP
                i = i + 1;
                xp = x0 + dx * random();
                yp = y0 + dy * random();
                rpoint = ST_SetSRID( ST_MakePoint( xp, yp ), ST_SRID(geom) );
                EXIT WHEN ST_Within( rpoint, geom );
        END LOOP;
        
        IF i >= maxiter THEN
                RAISE EXCEPTION 'RandomPoint: number of interations exceeded %', maxiter;
        END IF; 
        
        RETURN rpoint;
END; 
$$ LANGUAGE plpgsql;

This function takes too many iterations if the area of the polygon is small compared to the area of its bounding box. Practical example would be to place a random point inside any of the U.S. islands. The function below somewhat addresses this problem for the case of multipolygons:

CREATE OR REPLACE FUNCTION RandomPointMulti (
                geom Geometry
        )
        RETURNS Geometry
        AS $$
DECLARE
        maxiter INTEGER := 100000;
        i INTEGER := 0; 
        n INTEGER := 0; -- total number of geometries in collection
        g INTEGER := 0; -- geometry number in collection to find random point in
        total_area DOUBLE PRECISION; -- total area
        cgeom Geometry;
BEGIN
        total_area = ST_Area(geom);
        n = ST_NumGeometries(geom);
        
        WHILE i < maxiter LOOP
                i = i + 1;
                g = floor(random() * n)::int;
                cgeom = ST_GeometryN(geom, g); -- weight the probability of selecting a subpolygon by its relative area
                IF random() < ST_Area(cgeom)/total_area THEN
                        RETURN RandomPoint( cgeom );
                END IF;
        END LOOP;
        
        RAISE EXCEPTION 'RandomPointMulti: too many iterations';
END; 
$$ LANGUAGE plpgsql;

However, there are still cases in which the function above fails. For example, to create a random point within 1 mi of the US coast one would first build a 1-mile one-sided buffer along the coast line. The resulting geometry would look like two long and narrow polygons along the East and West coast with a huge empty space in between. To handle a case like that a function that would use some other strategy is necessary.

Random Points in Polygon

This set-returning function generates num_points uniform random points within a polygon geometry geom. See examples below.

CREATE OR REPLACE FUNCTION RandomPointsInPolygon(geom geometry, num_points integer)
  RETURNS SETOF geometry AS
$BODY$DECLARE
  target_proportion numeric;
  n_ret integer := 0;
  loops integer := 0;
  x_min float8;
  y_min float8;
  x_max float8;
  y_max float8;
  srid integer;
  rpoint geometry;
BEGIN
  -- Get envelope and SRID of source polygon
  SELECT ST_XMin(geom), ST_YMin(geom), ST_XMax(geom), ST_YMax(geom), ST_SRID(geom)
    INTO x_min, y_min, x_max, y_max, srid;
  -- Get the area proportion of envelope size to determine if a
  -- result can be returned in a reasonable amount of time
  SELECT ST_Area(geom)/ST_Area(ST_Envelope(geom)) INTO target_proportion;
  RAISE DEBUG 'geom: SRID %, NumGeometries %, NPoints %, area proportion within envelope %',
                srid, ST_NumGeometries(geom), ST_NPoints(geom),
                round(100.0*target_proportion, 2) || '%';
  IF target_proportion < 0.0001 THEN
    RAISE EXCEPTION 'Target area proportion of geometry is too low (%)', 
                    100.0*target_proportion || '%';
  END IF;
  RAISE DEBUG 'bounds: % % % %', x_min, y_min, x_max, y_max;
  
  WHILE n_ret < num_points LOOP
    loops := loops + 1;
    SELECT ST_SetSRID(ST_MakePoint(random()*(x_max - x_min) + x_min,
                                   random()*(y_max - y_min) + y_min),
                      srid) INTO rpoint;
    IF ST_Contains(geom, rpoint) THEN
      n_ret := n_ret + 1;
      RETURN NEXT rpoint;
    END IF;
  END LOOP;
  RAISE DEBUG 'determined in % loops (% efficiency)', loops, round(100.0*num_points/loops, 2) || '%';
END$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 1000;
ALTER FUNCTION RandomPointsInPolygon(geometry, integer) OWNER TO postgres;

Examples

-- Return 10 random points/rows
SELECT ST_AsText(RandomPointsInPolygon('POLYGON ((10 20, 30 60, 50 20, 10 20))', 10));

-- Or 10 random points as a single MultiPoint
SELECT ST_AsText(ST_Union(geom))
FROM RandomPointsInPolygon('POLYGON ((10 20, 30 60, 50 20, 10 20))', 10) AS geom

Now for a more complicated and practical example using a geopolitical polygon of the ​world table, create 10000 random points in India:

CREATE TABLE india_random_points
(
  gid serial NOT NULL,
  CONSTRAINT india_random_points_pkey PRIMARY KEY (gid)
);
SELECT AddGeometryColumn('india_random_points', 'geom', 4326, 'POINT', 2);

-- Populate table with 10000 random points within India
INSERT INTO india_random_points(geom)
SELECT RandomPointsInPolygon(geom, 10000)
FROM world WHERE name='INDIA';