= 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.
{{{
#!sql
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:
{{{
#!sql
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.