= Functions to Generate Random Points =
== 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.


== Random Points in Polygon ==
This set-returning function generates {{{num_points}}} uniform random points within a polygon geometry {{{geom}}}. See examples below.
{{{
#!sql
CREATE OR REPLACE FUNCTION RandomPointsInPolygon(geom geometry, num_points integer)
  RETURNS SETOF geometry AS
$BODY$DECLARE
  target_proportion numeric;
  n_ret integer := 0;
  loops integer := 0;
  x_min float8;
  y_min float8;
  x_max float8;
  y_max float8;
  srid integer;
  rpoint geometry;
BEGIN
  -- Get envelope and SRID of source polygon
  SELECT ST_XMin(geom), ST_YMin(geom), ST_XMax(geom), ST_YMax(geom), ST_SRID(geom)
    INTO x_min, y_min, x_max, y_max, srid;
  -- Get the area proportion of envelope size to determine if a
  -- result can be returned in a reasonable amount of time
  SELECT ST_Area(geom)/ST_Area(ST_Envelope(geom)) INTO target_proportion;
  RAISE DEBUG 'geom: SRID %, NumGeometries %, NPoints %, area proportion within envelope %',
                srid, ST_NumGeometries(geom), ST_NPoints(geom),
                round(100.0*target_proportion, 2) || '%';
  IF target_proportion < 0.0001 THEN
    RAISE EXCEPTION 'Target area proportion of geometry is too low (%)', 
                    100.0*target_proportion || '%';
  END IF;
  RAISE DEBUG 'bounds: % % % %', x_min, y_min, x_max, y_max;
  
  WHILE n_ret < num_points LOOP
    loops := loops + 1;
    SELECT ST_SetSRID(ST_MakePoint(random()*(x_max - x_min) + x_min,
                                   random()*(y_max - y_min) + y_min),
                      srid) INTO rpoint;
    IF ST_Contains(geom, rpoint) THEN
      n_ret := n_ret + 1;
      RETURN NEXT rpoint;
    END IF;
  END LOOP;
  RAISE DEBUG 'determined in % loops (% efficiency)', loops, round(100.0*num_points/loops, 2) || '%';
END$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 1000;
ALTER FUNCTION RandomPointsInPolygon(geometry, integer) OWNER TO postgres;
}}}

=== Examples ===
{{{
#!sql
-- Return 10 random points/rows
SELECT ST_AsText(RandomPointsInPolygon('POLYGON ((10 20, 30 60, 50 20, 10 20))', 10));

-- Or 10 random points as a single MultiPoint
SELECT ST_AsText(ST_Union(geom))
FROM RandomPointsInPolygon('POLYGON ((10 20, 30 60, 50 20, 10 20))', 10) AS geom
}}}

Now for a more complicated and practical example using a geopolitical polygon of the [http://www.aprsworld.net/gisdata/world/uncompressed/ world] table, create 10000 random points in India:
{{{
#!sql
CREATE TABLE india_random_points
(
  gid serial NOT NULL,
  CONSTRAINT india_random_points_pkey PRIMARY KEY (gid)
);
SELECT AddGeometryColumn('india_random_points', 'geom', 4326, 'POINT', 2);

-- Populate table with 10000 random points within India
INSERT INTO india_random_points(geom)
SELECT RandomPointsInPolygon(geom, 10000)
FROM world WHERE name='INDIA';
}}}