How to create random points in a polygon in PostGIS

Question

I have a polygon layer, and would like to generate random points in each polygon. An initial searching in Google resulted this website which presents a function for this purpose:

http://trac.osgeo.org/postgis/wiki/UserWikiRandomPoint

The function I'm interested to use is called RandomPointsInPolygon()

When I use this function and execute my query, it doesn't give an error but also never finishes the process. I tried giving a WHERE clause and limiting my query, this time process executed without any error but no points were stored in the relevant field in my table.

Does anybody know what the problem is? and/or what could be an easy way to solve my problem? using this function or any other options that you may have know?

The code for function:

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;

My query:

INSERT INTO grid(point_geom) SELECT RandomPointsInPolygon(h.geom, 1) FROM grid AS h;

Here I have a point geometry column column named point_geom to store the result of query, and in the same table I have polygons stored in the geom field. I would like to generate points inside each polygon of this grid table. In the query above I tried to create 1 random point for each polgyon and store it in the relevant field.

Answer 1

To create 1 random point for each polgyon in a table I have used this table and code:

enter image description here

DO $$
DECLARE
 r RECORD;
BEGIN
 FOR r IN SELECT id_0 FROM "Grid" LOOP
-- RAISE NOTICE 'affected row id: %', r.id_0;
 UPDATE "Grid" SET "point_geom" = (SELECT RandomPointsInPolygon(geom, 1) FROM "Grid" WHERE "id_0" = r.id_0) WHERE "id_0" = r.id_0;
 END LOOP;
END$$;

Running all the code at once, I've obtained the following:

enter image description here

enter image description here

enter image description here

Answer 2

It's probably because you have MultiPolygons in your geometry, which are collections of polygons. The author followed up with a function RandomPointMulti that handled the known issue with MultiPolygons.

If you're not sure what geometry type you're working with (Points, Linestring, Polygon, MultiPoints, MultiLinestring, or MultiPolygon), I highly recommend you become familiar with the great Boundless Geo documentation on them.

• Catlover: No, I have simple polygons. Any other suggestion?

Answer 3

Have you set the SRID on the geometries? This code works for me (I have multipolygons, hence the use of ST_Dump):

SELECT RandomPointsInPolygon((ST_Dump(ST_SetSRID(geom, 4283))).geom, 50)

On my home PC it can do ~20 million points in under an hour based on ~350,000 not too complex polygons.