Cómo crear puntos aleatorios en un polígono en PostGIS

Question

Tengo una capa de polígonos, y me gustaría generar puntos aleatorios en cada polígono. Una búsqueda inicial en Google resultó este sitio web que presenta una función para este propósito:

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

La función que me interesa utilizar se llama PuntosAleatoriosEnPolígono()

Cuando uso esta función y ejecuto mi consulta, no da error pero tampoco termina nunca el proceso. Probé dando una cláusula WHERE y limitando mi consulta, esta vez el proceso se ejecutó sin ningún error pero no se almacenó ningún punto en el campo correspondiente de mi tabla.

¿Alguien sabe cuál es el problema? y/o ¿cuál podría ser una forma fácil de resolver mi problema? ¿usando esta función o cualquier otra opción que conozca?

El código de la función:

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;

Mi pregunta:

INSERT INTO grid(punto_geom) SELECT RandomPointsInPolygon(h.geom, 1) FROM rejilla COMO h;

Aquí tengo una columna de geometría de puntos llamada punto_geom para almacenar el resultado de la consulta, y en la misma tabla tengo polígonos almacenados en el geom campo. Me gustaría generar puntos dentro de cada polígono de este cuadrícula mesa. En la consulta anterior he intentado crear 1 punto aleatorio para cada poligono y almacenarlo en el campo correspondiente.

Answer 1

Para crear 1 punto aleatorio por cada polgyon en una tabla he utilizado esta tabla y este código:

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$$;

Ejecutando todo el código a la vez, he obtenido lo siguiente:

Answer 2

Probablemente sea porque tiene MultiPolígonos en su geometría, que son colecciones de polígonos. El autor siguió con una función RandomPointMulti que trataba el problema conocido con MultiPolygons.

Si no está seguro de con qué tipo de geometría está trabajando (Puntos, Cadena de líneas, Polígono, MultiPuntos, MultiCadena de líneas o MultiPolígono), le recomiendo encarecidamente que se familiarice con la magnífica herramienta Documentación de Boundless Geo en ellos.

Answer 3

¿Ha configurado el SRID en las geometrías? Este código me funciona (tengo multipolígonos, de ahí el uso de ST_Dump):

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

En mi PC doméstico puede hacer ~20 millones de puntos en menos de una hora a partir de ~350.000 polígonos no demasiado complejos.