Creación de una cuadrícula de puntos regulares dentro de un polígono en PostGIS

Question

¿Cómo crear, dentro de un polígono, una cuadrícula regular de puntos espaciados x,y en PostGIS?

Como el ejemplo:

Answer 1

Eso se hace con generate_series.

Si no quieres escribir manualmente dónde debe empezar y parar la cuadrícula, lo más fácil es crear una función.

No he probado bien lo siguiente, pero creo que debería funcionar:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_POINT(x, y)) FROM 
generate_series(floor(ST_XMIN($1))::int, ceiling(ST_XMAX($1)-ST_XMIN($1))::int, $2) AS x,
generate_series(floor(ST_YMIN($1))::int, ceiling(ST_YMAX($1)-ST_YMIN($1))::int, $2) AS y 
WHERE st_intersects($1, ST_POINT(x, y))'
LANGUAGE sql

Para utilizarlo puedes hacerlo:

SELECT makegrid(the_geom, 1000) from mytable;

donde el primer argumento es el polígono en el que quieres la cuadrícula, y el segundo argumento es la distancia entre los puntos de la cuadrícula.

Si quieres un punto por fila sólo tienes que usar ST_Dump como:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

Answer 2

He recogido Nicklas Avén código de la función makegrid y lo hizo un poco más genérico mediante la lectura y el uso de la srid de la geometría del polígono. De lo contrario, el uso de un polígono con un srid definido, daría un error.

La función:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_SetSRID(ST_POINT(x, y), ST_SRID($1))) FROM
generate_series(floor(ST_XMIN($1))::int, ceiling(ST_XMAX($1) - ST_XMIN($1))::int, $2) AS x,
generate_series(floor(ST_YMIN($1))::int, ceiling(ST_YMAX($1) - ST_YMIN($1))::int, $2) AS y
WHERE st_intersects($1, ST_SetSRID(ST_POINT(x,y), ST_SRID($1)))'
LANGUAGE sql

Para utilizar la función se hace exactamente como Nicklas Avén escribió:

SELECT makegrid(the_geom, 1000) from mytable;

o si quiere un punto por fila:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

Espero que esto sea útil para alguien.

Alex

Answer 3

Las personas que utilizan una geometría wgs84 probablemente tendrán problemas con esta función ya que

generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1))::int,$2) as y 

sólo devuelven números enteros. Excepto en el caso de geometrías muy grandes, como los países (que se encuentran en varios grados lat, lng), esto hará que se recoja sólo un punto que la mayoría de las veces ni siquiera se cruza con la propia geometría... => ¡resultado vacío!

Mi problema fue que no puedo usar generate_series() con una distancia decimal en números flotantes como los WSG84... Es por eso que he ajustado la función para conseguir que funcione de todos modos :

SELECT ST_Collect(st_setsrid(ST_POINT(x/1000000::float,y/1000000::float),st_srid($1))) FROM 
  generate_series(floor(st_xmin($1)*1000000)::int, ceiling(st_xmax($1)*1000000)::int,$2) as x ,
  generate_series(floor(st_ymin($1)*1000000)::int, ceiling(st_ymax($1)*1000000)::int,$2) as y 
WHERE st_intersects($1,ST_SetSRID(ST_POINT(x/1000000::float,y/1000000::float),ST_SRID($1)))

Básicamente es lo mismo. Sólo multiplicar y dividir por 1000000 para obtener los decimales en el juego cuando lo necesito.

Seguro que hay una solución mejor para conseguirlo. ++

Answer 4

Tres algoritmos que utilizan diferentes métodos.

Repo de Github: https://github.com/imran-5/Postgis-Custom

1. El enfoque más sencillo y mejor, utilizando la distancia terrestre real de las coordenadas de la dirección x e y. El algoritmo funciona con cualquier SRID, internamente funciona con WGS 1984(EPSG:4326), y el resultado se transforma de nuevo a SRID de entrada.

Function ===================================================================

CREATE OR REPLACE FUNCTION public.I_Grid_Point_Distance(geom public.geometry, x_side decimal, y_side decimal)
RETURNS public.geometry AS $BODY$
DECLARE
x_min decimal;
x_max decimal;
y_max decimal;
x decimal;
y decimal;
returnGeom public.geometry[];
i integer := -1;
srid integer := 4326;
input_srid integer;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
    geom := ST_SetSRID(geom, srid);
        ----RAISE NOTICE 'No SRID Found.';
    ELSE
        ----RAISE NOTICE 'SRID Found.';
END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_min := ST_XMin(geom);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    y := ST_YMin(geom);
    x := x_min;
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
<<yloop>>
LOOP
IF (y > y_max) THEN
    EXIT;
END IF;

CASE i WHEN 0 THEN 
    y := ST_Y(returnGeom[0]);
ELSE 
    y := ST_Y(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), y_side, radians(0))::geometry);
END CASE;

x := x_min;
<<xloop>>
LOOP
  IF (x > x_max) THEN
      EXIT;
  END IF;
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
    x := ST_X(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), x_side, radians(90))::geometry);
END LOOP xloop;
END LOOP yloop;
RETURN
ST_CollectionExtract(st_transform(ST_Intersection(st_collect(returnGeom), geom), input_srid), 1);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE;

Utilice la función con una consulta simple, la geometría debe ser válida y de tipo polígono, multipolígono o envolvente

SELECT I_Grid_Point_Distance(geom, 50, 61) from polygons limit 1;

Result======================================================================

2. Segunda función basada en Nicklas Avén algoritmo. He implementado una forma de manejar cualquier SRID.

   han actualizado los siguientes cambios en el algoritmo.

   1. Variable separada para la dirección x e y para el tamaño del píxel,
   2. Nueva variable para calcular la distancia en esferoide o elipsoide.
   3. Introducir cualquier SRID, función transformar Geom al entorno de trabajo de Esferoide o Elipsoide Datum, luego aplicar la distancia a cada lado, obtener el resultado y transformarlo en el SRID de entrada.

Function ===================================================================

CREATE OR REPLACE FUNCTION I_Grid_Point(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$ 
DECLARE
x_max decimal; 
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer; 
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
        srid := 4326;
        x_side := x_side / 100000;
        y_side := y_side / 100000;
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);
RETURN QUERY
WITH res as (SELECT ST_SetSRID(ST_MakePoint(x, y), srid) point FROM
generate_series(x_min, x_max, x_side) as x,
generate_series(y_min, y_max, y_side) as y
WHERE st_intersects(geom, ST_SetSRID(ST_MakePoint(x, y), srid))
) select ST_TRANSFORM(ST_COLLECT(point), input_srid) from res;
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Utilízalo con una simple consulta.

SELECT I_Grid_Point(geom, 22, 15, false) from polygons;

Result===================================================================

3. Función basada en un generador en serie.

Function==================================================================

CREATE OR REPLACE FUNCTION I_Grid_Point_Series(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer;
x_series DECIMAL;
y_series DECIMAL;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);

    x_series := CEIL ( @( x_max - x_min ) / x_side);
    y_series := CEIL ( @( y_max - y_min ) / y_side );
RETURN QUERY
SELECT st_collect(st_setsrid(ST_MakePoint(x * x_side + x_min, y*y_side + y_min), srid)) FROM
generate_series(0, x_series) as x,
generate_series(0, y_series) as y
WHERE st_intersects(st_setsrid(ST_MakePoint(x*x_side + x_min, y*y_side + y_min), srid), geom);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Utilízalo con una simple consulta.

SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons; Result==========================================================================

Answer 5

Este algoritmo debería estar bien:

createGridInPolygon(polygon, resolution) {
    for(x=polygon.xmin; x<polygon.xmax; x+=resolution) {
       for(y=polygon.ymin; y<polygon.ymax; y+=resolution) {
          if(polygon.contains(x,y)) createPoint(x,y);
       }
    }
}

donde "polígono" es el polígono y "resolución" es la resolución de malla requerida.

Para implementarlo en PostGIS, pueden ser necesarias las siguientes funciones:

• ST_XMin , ST_XMax , ST_YMin y ST_YMax para obtener las coordenadas mínimas y máximas del polígono,
• ST_Contiene para comprobar si el polígono contiene el punto,
• y Punto_ST para crear un punto.

Buena suerte.