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

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 convierte la geometría del polígono en una trama vacía, y luego convierte cada píxel en un punto. Ventaja: no tenemos que comprobar de nuevo si el polígono original interseca los puntos.

Opcional:

También puede añadir la alineación de la rejilla con los parámetros 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 de 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 utilizar con:

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

Answer 2

Así que mi versión arreglada:

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

Answer 3

Una pequeña actualización potencial a las respuestas anteriores - tercer argumento como escala para wgs84 (o usar 1 para los normales), y también redondear dentro del código para que los puntos escalados en múltiples 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

Answer 4

Cree primero la función:

CREATE OR REPLACE FUNCTION st_polygrid(geometry, integer) RETURNS geometry AS
$$
SELECT
    ST_SetSRID(ST_Collect(ST_POINT(x,y)), ST_SRID($1))
FROM
  generate_series(floor(ST_XMin($1))::numeric, ceiling(ST_xmax($1))::numeric, $2) as x,
  generate_series(floor(ST_ymin($1))::numeric, ceiling(ST_ymax($1))::numeric,$2) as y 
WHERE
  ST_Intersects($1,ST_SetSRID(ST_POINT(x,y), ST_SRID($1)))
$$
  LANGUAGE sql VOLATILE;

A continuación, cree la cuadrícula de puntos utilizando esa función:

CREATE TABLE some.other_table AS
SELECT
  a.gid,
  st_polygrid(a.the_geom, 1000) AS the_geom
FROM
  some.table a

cambiando el '1000' al espacio de puntos deseado.

Answer 5

No necesitaba una función, sólo una consulta que pudiera ejecutar. Basado en algunas de las respuestas aquí, hice lo siguiente:

WITH geom AS (
    SELECT st_setsrid('polygon-string'::geometry, 4326) AS geometry
),
points AS (
    SELECT generate_series(ST_XMIN(ST_Transform(geometry, 3857))::int, ST_XMAX(ST_Transform(geometry, 3857))::int, 1000) AS x,
    generate_series(ST_YMIN(ST_Transform(geometry, 3857))::int, ST_YMAX(ST_Transform(geometry, 3857))::int, 1000) AS y
    FROM geom
)
SELECT ST_Collect(ST_Transform(st_setsrid(ST_POINT(x, y), 3857), 4326)) AS points
FROM points, geom
WHERE st_intersects(geometry, ST_Transform(st_setsrid(ST_POINT(x, y), 3857), 4326));

1. El geom CTE convierte la cadena en geometría y se asegura de que su SRID sea 4326
2. El points El CTE formatea la geometría en 3857, y crea dos listas de puntos separados por 1000m
3. La consulta final transforma los puntos de nuevo en 4326, filtra los puntos que no están dentro de la forma del polígono y los recoge en una colección de geometría