| 1 | = Random Point in Polygon = |
| 2 | |
| 3 | Use this function when you need to create a random point inside a |
| 4 | polygon. The function takes geometry as a required argument and max |
| 5 | number of iterations as an optional argument (defaults to 1000 |
| 6 | iterations). The function will raise an exception if max number of |
| 7 | iterations has been exceeded. The function works by throwing a random |
| 8 | point within the polygon bounding rectangle and then checking if the |
| 9 | point is inside the polygon. |
| 10 | |
| 11 | {{{ |
| 12 | #!sql |
| 13 | CREATE OR REPLACE FUNCTION RandomPoint ( |
| 14 | geom Geometry, |
| 15 | maxiter INTEGER DEFAULT 1000 |
| 16 | ) |
| 17 | RETURNS Geometry |
| 18 | AS $$ |
| 19 | DECLARE |
| 20 | i INTEGER := 0; |
| 21 | x0 DOUBLE PRECISION; |
| 22 | dx DOUBLE PRECISION; |
| 23 | y0 DOUBLE PRECISION; |
| 24 | dy DOUBLE PRECISION; |
| 25 | xp DOUBLE PRECISION; |
| 26 | yp DOUBLE PRECISION; |
| 27 | rpoint Geometry; |
| 28 | BEGIN |
| 29 | -- find envelope |
| 30 | x0 = ST_XMin(geom); |
| 31 | dx = (ST_XMax(geom) - x0); |
| 32 | y0 = ST_YMin(geom); |
| 33 | dy = (ST_YMax(geom) - y0); |
| 34 | |
| 35 | WHILE i < maxiter LOOP |
| 36 | i = i + 1; |
| 37 | xp = x0 + dx * random(); |
| 38 | yp = y0 + dy * random(); |
| 39 | rpoint = ST_SetSRID( ST_MakePoint( xp, yp ), ST_SRID(geom) ); |
| 40 | EXIT WHEN ST_Within( rpoint, geom ); |
| 41 | END LOOP; |
| 42 | |
| 43 | IF i >= maxiter THEN |
| 44 | RAISE EXCEPTION 'RandomPoint: number of interations exceeded ', maxiter; |
| 45 | END IF; |
| 46 | |
| 47 | RETURN rpoint; |
| 48 | END; |
| 49 | $$ LANGUAGE plpgsql; |
| 50 | }}} |
| 51 | |
| 52 | This function takes too many iterations if the area of the polygon is |
| 53 | small compared to the area of its bounding box. Practical example |
| 54 | would be to place a random point inside any of the U.S. islands. The |
| 55 | function below somewhat addresses this problem for the case of |
| 56 | multipolygons: |
| 57 | |
| 58 | {{{ |
| 59 | #!sql |
| 60 | CREATE OR REPLACE FUNCTION RandomPointMulti ( |
| 61 | geom Geometry |
| 62 | ) |
| 63 | RETURNS Geometry |
| 64 | AS $$ |
| 65 | DECLARE |
| 66 | maxiter INTEGER := 100000; |
| 67 | i INTEGER := 0; |
| 68 | n INTEGER := 0; -- total number of geometries in collection |
| 69 | g INTEGER := 0; -- geometry number in collection to find random point in |
| 70 | total_area DOUBLE PRECISION; -- total area |
| 71 | cgeom Geometry; |
| 72 | BEGIN |
| 73 | total_area = ST_Area(geom); |
| 74 | n = ST_NumGeometries(geom); |
| 75 | |
| 76 | WHILE i < maxiter LOOP |
| 77 | i = i + 1; |
| 78 | g = floor(random() * n)::int; |
| 79 | cgeom = ST_GeometryN(geom, g); -- weight the probability of selecting a subpolygon by its relative area |
| 80 | IF random() < ST_Area(cgeom)/total_area THEN |
| 81 | RETURN RandomPoint( cgeom ); |
| 82 | END IF; |
| 83 | END LOOP; |
| 84 | |
| 85 | RAISE EXCEPTION 'RandomPointMulti: too many iterations'; |
| 86 | END; |
| 87 | $$ LANGUAGE plpgsql; |
| 88 | }}} |
| 89 | |
| 90 | However, there are still cases in which the function above fails. For |
| 91 | example, to create a random point within 1 mi of the US coast one |
| 92 | would first build a 1-mile one-sided buffer along the coast line. The |
| 93 | resulting geometry would look like two long and narrow polygons along |
| 94 | the East and West coast with a huge empty space in between. To handle |
| 95 | a case like that a function that would use some other strategy is |
| 96 | necessary. |