|Version 2 (modified by 11 years ago) ( diff ),|
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.