Thursday 27 August 2015

How to create random points in a polygon in PostGIS


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




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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...