Changes between Initial Version and Version 1 of UserWikiRandomPoint


Ignore:
Timestamp:
May 31, 2011, 12:47:01 PM (12 years ago)
Author:
sorokine
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • UserWikiRandomPoint

    v1 v1  
     1= Random Point in Polygon =
     2
     3Use this function when you need to create a random point inside a
     4polygon.  The function takes geometry as a required argument and max
     5number of iterations as an optional argument (defaults to 1000
     6iterations).  The function will raise an exception if max number of
     7iterations has been exceeded.  The function works by throwing a random
     8point within the polygon bounding rectangle and then checking if the
     9point is inside the polygon.
     10
     11{{{
     12#!sql
     13CREATE OR REPLACE FUNCTION RandomPoint (
     14                geom Geometry,
     15                maxiter INTEGER DEFAULT 1000
     16        )
     17        RETURNS Geometry
     18        AS $$
     19DECLARE
     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;
     28BEGIN
     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;
     48END;
     49$$ LANGUAGE plpgsql;
     50}}}
     51
     52This function takes too many iterations if the area of the polygon is
     53small compared to the area of its bounding box.  Practical example
     54would be to place a random point inside any of the U.S. islands.  The
     55function below somewhat addresses this problem for the case of
     56multipolygons:
     57
     58{{{
     59#!sql
     60CREATE OR REPLACE FUNCTION RandomPointMulti (
     61                geom Geometry
     62        )
     63        RETURNS Geometry
     64        AS $$
     65DECLARE
     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;
     72BEGIN
     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';
     86END;
     87$$ LANGUAGE plpgsql;
     88}}}
     89
     90However, there are still cases in which the function above fails. For
     91example, to create a random point within 1 mi of the US coast one
     92would first build a 1-mile one-sided buffer along the coast line.  The
     93resulting geometry would look like two long and narrow polygons along
     94the East and West coast with a huge empty space in between.  To handle
     95a case like that a function that would use some other strategy is
     96necessary.