¿Cómo crear una cuadrícula de puntos dentro de un polígono en Postgis?

Lo haces con genera_series.

Si no desea escribir manualmente dónde se debe iniciar y detener la cuadrícula, lo más fácil es crear una función.

No he probado la siguiente información correctamente, 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 usarlo puedes hacer:

SELECT makegrid(the_geom, 1000) from mytable;

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

Si desea un punto por fila, simplemente use ST_Dump como:

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

HTH

Nicklas

He recogido el código de función makegrid Nicklas Avén y lo he hecho un poco más genérico leyendo y utilizando el srid de la geometría de polígono. De lo contrario, utilizar 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 usar la función, se hace exactamente como Nicklas Avén escribió:

SELECT makegrid(the_geom, 1000) from mytable;

o si quieres 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

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 cuadrícula requerida.

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

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

¡Buena suerte!

Las personas que usan 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 

solo devuelve enteros. Excepto en el caso de geometrías muy grandes como países (que se encuentran en latitud múltiple, grados de lng), esto hará que se acumule solo 1 punto, que la mayoría de las veces ni siquiera se interseca con la geometría en sí ... = > resultado vacío!

Mi problema fue que parece que no puedo usar generate_series () con una distancia decimal en números flotantes como los WSG84 ... Es por eso que modifiqué la función para 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 exactamente lo mismo. Solo multiplica y divide por 1000000 para obtener los decimales en el juego cuando lo necesite.

Sin duda hay una mejor solución para lograrlo. ++

Así que mi versión fija:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer)
RETURNS geometry AS
'SELECT ST_Collect(st_setsrid(ST_POINT(x,y),$3)) FROM 
  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 
where st_intersects($1,st_setsrid(ST_POINT(x,y),$3))'
LANGUAGE sql

Uso:

SELECT (ST_Dump(makegrid(the_geom, 1000, 3857))).geom as the_geom from my_polygon_table

Una pequeña actualización potencial de las respuestas anteriores: tercer argumento como escala para wgs84 (o uso 1 para las normales), y también se redondea dentro del código para que los puntos escalados en varias formas estén alineados.

Espero que esto ayude, Martin

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer)
RETURNS geometry AS



/*geometry column , integer: distance between points, integer: scale factor for distance (useful for wgs84, e.g. use there 50000 as distance and 1000000 as scale factor*/

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

LANGUAGE sql

Aquí hay otro enfoque que es ciertamente más rápido y más fácil de entender.

Por ejemplo, para una cuadrícula de 1000m por 1000m:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0))).geom 
FROM the_polygon

También se conserva el SRID original.

Este fragmento de código convierte la geometría del polígono en un ráster vacío, luego convierte cada píxel en un punto. Ventaja: no tenemos que verificar de nuevo si el polígono original se cruza con los puntos.

Opcional:

También puede agregar la alineación de la cuadrícula con el parámetro gridx y gridy. Pero como usamos el centroide de cada píxel (y no una esquina) necesitamos usar un módulo para calcular el valor correcto:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0,mod(1000/2,100),mod(1000/2,100)))).geom 
FROM the_polygon

Con mod(grid_size::integer/2,grid_precision)

Aquí está la función postgres:

CREATE OR REPLACE FUNCTION st_makegrid(geometry, float, integer)
RETURNS SETOF geometry AS
'SELECT (ST_PixelAsCentroids(ST_AsRaster($1,$2::float,$2::float,mod($2::int/2,$3),mod($2::int/2,$3)))).geom'
LANGUAGE sql;

Se puede usar con:

SELECT makegrid(the_geom,1000.0,100) as geom from the_polygon  
-- makegrid(the_geom,grid_size,alignement)

He creado dos algoritmos para la solución y he utilizado las funciones vectoriales de PostGIS.

1. Este es el mejor enfoque desde mi lado con la distancia en metros desde el píxel del lado xy el y. Ingrese la geometría, la distancia de x e y para la cuadrícula de puntos dentro de un polígono.

Funciona con casi cualquier SRID. El algoritmo interno funciona con WGS 1984 y devuelve el resultado de Geometría con Geometría SRID insertada.
He tratado de manejar todos los aspectos para coincidir con la distancia exacta de la Tierra y SRID de acuerdo con la Geometría de entrada y salida en este algoritmo.

---DROP FUNCTION public.makeGridPoint_2D_Distance(geometry, decimal, decimal);


CREATE OR REPLACE FUNCTION public.makeGridPoint_2D_Distance( 
geom public.geometry, 
x_step decimal,
y_step decimal ) RETURNS public.geometry AS $BODY$
DECLARE
xMin decimal;
xMax decimal;
yMax 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, 4326);
        ----RAISE NOTICE 'No SRID Found.';
    ELSE
        ----RAISE NOTICE 'SRID Found.';
END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    xMin := ST_XMin(geom);
    xMax := ST_XMax(geom);
    yMax := ST_YMax(geom);
    y := ST_YMin(geom);
    x := xMin; 
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
<<yloop>>
LOOP
IF (y > yMax) 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_step, radians(0))::geometry);
END CASE;

x := xMin;
<<xloop>>
LOOP
  IF (x > xMax) 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_step, 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 esta función con la consulta.

SELECT makegrid_2d(geom, 20, 20) from polygons

o si quieres un solo punto:

SELECT (ST_Dump(makegrid_2d(geom, 20, 20))).geom from polygons

2. Mi segundo enfoque es usar el algoritmo Nicklas Avén . He mejorado el algoritmo para manejar casi cualquier SRID.

estos cambios los he aplicado en el algoritmo original.

• variables separadas para la dirección x e y del tamaño de píxel,
• nueva variable para usar la distancia en esferoide o elipsoide.
• ingrese cualquier geometría SRID que se aplicará a la transformación al entorno de trabajo de Spheroid o Ellipsoid Datum, luego aplique la distancia y regrese la Geometría con SRID insertado de Geometría

Aquí está la función

CREATE OR REPLACE FUNCTION makePointGrid(geom geometry, x_step decimal, y_step decimal, spheroid boolean default true) 
RETURNS SETOF geometry AS $
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer;
input_srid integer;
BEGIN 
    CASE st_srid(geom) WHEN 0 THEN
        geom := ST_SetSRID(geom, 4326);
        RAISE NOTICE 'SRID Found Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;
CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid Not Found';
        srid := 4326;
        x_step := x_step / 100000;
        y_step := y_step / 100000;
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid';
    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
SELECT st_transform(ST_Collect(ST_SetSRID(ST_POINT(x, y), srid)), input_srid) FROM 
generate_series(x_min, x_max, x_step) as x,
generate_series(y_min, y_max, y_step) as y
WHERE st_intersects(ST_SetSRID(geom, srid), ST_SetSRID(ST_POINT(x, y),
srid));END;$BODY$ LANGUAGE plpgsql IMMUTABLE;

Utilice esta función con la consulta.

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

o si quieres un solo punto:

SELECT (ST_Dump(makePointGrid(geom, 20, 20, false))).geom from polygons