Changes between Initial Version and Version 1 of UserWikiRandomPoint

Timestamp:

May 31, 2011, 12:47:01 PM (13 years ago)

Author:

sorokine

Comment:

—

UserWikiRandomPoint

@@ +1 @@
+= 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.