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.