Changes between Version 2 and Version 3 of UserWikiRandomPoint


Ignore:
Timestamp:
Jun 7, 2011, 5:11:39 AM (10 years ago)
Author:
Mike Taves
Comment:

add set returning function

Legend:

Unmodified
Added
Removed
Modified
  • UserWikiRandomPoint

    v2 v3  
    1 = Random Point in Polygon =
     1= Functions to Generate Random Points =
     2== Random Point in Polygon ==
    23
    34Use this function when you need to create a random point inside a
     
    9495the East and West coast with a huge empty space in between.  To handle
    9596a case like that a function that would use some other strategy is
    96 necessary.
     97necessary.
     98
     99
     100== Random Points in Polygon ==
     101This set-returning function generates {{{num_points}}} uniform random points within a polygon geometry {{{geom}}}. See examples below.
     102{{{
     103#!sql
     104CREATE OR REPLACE FUNCTION RandomPointsInPolygon(geom geometry, num_points integer)
     105  RETURNS SETOF geometry AS
     106$BODY$DECLARE
     107  target_proportion numeric;
     108  n_ret integer := 0;
     109  loops integer := 0;
     110  x_min float8;
     111  y_min float8;
     112  x_max float8;
     113  y_max float8;
     114  srid integer;
     115  rpoint geometry;
     116BEGIN
     117  -- Get envelope and SRID of source polygon
     118  SELECT ST_XMin(geom), ST_YMin(geom), ST_XMax(geom), ST_YMax(geom), ST_SRID(geom)
     119    INTO x_min, y_min, x_max, y_max, srid;
     120  -- Get the area proportion of envelope size to determine if a
     121  -- result can be returned in a reasonable amount of time
     122  SELECT ST_Area(geom)/ST_Area(ST_Envelope(geom)) INTO target_proportion;
     123  RAISE DEBUG 'geom: SRID %, NumGeometries %, NPoints %, area proportion within envelope %',
     124                srid, ST_NumGeometries(geom), ST_NPoints(geom),
     125                round(100.0*target_proportion, 2) || '%';
     126  IF target_proportion < 0.0001 THEN
     127    RAISE EXCEPTION 'Target area proportion of geometry is too low (%)',
     128                    100.0*target_proportion || '%';
     129  END IF;
     130  RAISE DEBUG 'bounds: % % % %', x_min, y_min, x_max, y_max;
     131 
     132  WHILE n_ret < num_points LOOP
     133    loops := loops + 1;
     134    SELECT ST_SetSRID(ST_MakePoint(random()*(x_max - x_min) + x_min,
     135                                   random()*(y_max - y_min) + y_min),
     136                      srid) INTO rpoint;
     137    IF ST_Contains(geom, rpoint) THEN
     138      n_ret := n_ret + 1;
     139      RETURN NEXT rpoint;
     140    END IF;
     141  END LOOP;
     142  RAISE DEBUG 'determined in % loops (% efficiency)', loops, round(100.0*num_points/loops, 2) || '%';
     143END$BODY$
     144  LANGUAGE plpgsql VOLATILE
     145  COST 100
     146  ROWS 1000;
     147ALTER FUNCTION RandomPointsInPolygon(geometry, integer) OWNER TO postgres;
     148}}}
     149
     150=== Examples ===
     151{{{
     152#!sql
     153-- Return 10 random points/rows
     154SELECT ST_AsText(RandomPointsInPolygon('POLYGON ((10 20, 30 60, 50 20, 10 20))', 10));
     155
     156-- Or 10 random points as a single MultiPoint
     157SELECT ST_AsText(ST_Union(geom))
     158FROM RandomPointsInPolygon('POLYGON ((10 20, 30 60, 50 20, 10 20))', 10) AS geom
     159}}}
     160
     161Now 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:
     162{{{
     163#!sql
     164CREATE TABLE india_random_points
     165(
     166  gid serial NOT NULL,
     167  CONSTRAINT india_random_points_pkey PRIMARY KEY (gid)
     168);
     169SELECT AddGeometryColumn('india_random_points', 'geom', 4326, 'POINT', 2);
     170
     171-- Populate table with 10000 random points within India
     172INSERT INTO india_random_points(geom)
     173SELECT RandomPointsInPolygon(geom, 10000)
     174FROM world WHERE name='INDIA';
     175}}}