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

31

¿Cómo crear, dentro de un polígono, una cuadrícula regular de puntos espaciados x, y en postgis? Como el ejemplo:

texto alternativo

Pablo
fuente
Intenté hacer que los polígonos de recorte fusionaran este código con el código de corte "postGIS en acción", pero solo se crea un polígono ... ¿Olvidé algo? CREAR O REEMPLAZAR LA FUNCIÓN makegrid (geometry, integer, integer) RETURNS geometry AS 'SELECT st_intersection (g1.geom1, g2.geom2) AS geom FROM (SELECT $ 1 AS geom1) AS g1 INNER JOIN (Select st_setsrid (CAST (ST_MakeBox2d (st_setsrid (st_setsrid (st_setsrid (st_setsrid) ST_Point (x, y), $ 3), st_setsrid (ST_Point (x + $ 2, y + $ 2), $ 3)) como geometría), $ 3) como geom2 FROM generate_series (floor (st_xmin ($ 1)) :: int, ceiling ( st_xmax ($ 1)) :: int, $ 2) como x, generate_series (floor (st_ymin ($ 1)) :: int, ceiling (st_ymax (
aurel_nc
mira mi respuesta detallada.
Muhammad Imran Siddique

Respuestas:

30

Lo haces con generate_series.

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

No he probado el siguiente 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 desea la cuadrícula, y el segundo argumento es la distancia entre los puntos en 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

Nicklas Avén
fuente
2
Es posible que deba agregar st_setSRID () a las funciones st_point, de lo contrario, st_intersects no funciona.
JaakL
Agregué mi versión probada como respuesta separada.
JaakL
12

He recogido Nicklas Aven código de función makegrid y lo hizo un poco más genérica mediante la lectura y el uso de la srid de la geometría de polígono. De lo contrario, usar un polígono con una cuadrícula definida 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

Alexandre Neto
fuente
La respuesta aceptada no funciona con mis datos debido a errores SRID. Esta modificación funciona mejor.
Vitaly Isaev
¿Es posible que desee agregar algo cuando el polígono cruza el antimeridiano? Me imagino que podría causar un problema con xmin / xmax.
Thomas
2
No me funcionó. Usando Postgres 9.6 y PostGIS 2.3.3. Dentro de la llamada generate_series, tuve que poner esto como el segundo parámetro "ceiling (st_xmax ($ 1)) :: int" en lugar de "ceiling (st_xmax ($ 1) -st_xmin ($ 1)) :: int", y "ceiling ( st_ymax ($ 1)) :: int "en lugar de" techo (st_ymax ($ 1) -st_ymin ($ 1)) :: int "
Vitor Sapucaia
Apruebo el comentario anterior; el límite superior de generate_series debe ser el techo del máximo, y no el techo de la diferencia (max - min).
R. Bourgeon
10

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. A excepción de las geometrías muy grandes, como los países (que se encuentran en múltiples grados lat, lng), esto hará que se acumule solo 1 punto, que la mayoría de las veces ni siquiera se cruza con la geometría misma ... => ¡resultado vacío!

Mi problema fue que parece que no puedo usar generate_series () con una distancia decimal en números flotantes como esos 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. Simplemente multiplique y divida por 1000000 para obtener los decimales en el juego cuando lo necesite.

Seguramente hay una mejor solución para lograrlo. ++

Julien Garcia
fuente
Esa es una solución inteligente. ¿Has revisado los resultados? ¿Son consistentes?
Pablo
Hola. Si pablo Estoy contento con los resultados hasta ahora. Lo necesitaba para construir un polígono con una altitud relativa sobre el suelo. (Uso SRTM para calcular la altitud que quiero para cada punto de la cuadrícula). Ahora solo me falta una forma de incluir los puntos que se encuentran en el perímetro del polígono también. Actualmente, la forma renderizada está algo truncada en el borde.
Julien Garcia
funcionó, cuando todas las soluciones de los demás fallaron, ¡gracias!
Jordan Arseno
7

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, pueden ser necesarias las siguientes funciones:

¡Buena suerte!

julien
fuente
1
Tenga en cuenta que si tiene grandes áreas de polígonos complejos (p. Ej., Tengo un buffer de costa), entonces este enfoque no es del todo óptimo.
JaakL
Entonces, ¿qué propones en su lugar?
julien
4

Tres algoritmos que utilizan diferentes métodos.

Enlace de repositorio de Github

  1. El enfoque simple y mejor, utilizando la distancia real a la tierra de las coordenadas desde la dirección x e y. El algoritmo funciona con cualquier SRID, internamente funciona con WGS 1984 (EPSG: 4326) y el resultado se transforma nuevamente en SRID de entrada.

Función ================================================= ==================

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;

Use la función con una consulta simple, la geometría debe ser válida y el tipo de polígono, multi-polígonos o envolvente

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

Resultado ================================================= =====================

ingrese la descripción de la imagen aquí

  1. Segunda función basada en el algoritmo de Nicklas Avén . Lo mejoré para 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 de píxel,
    2. Nueva variable para calcular la distancia en esferoide o elipsoide.
    3. Ingrese cualquier SRID, la función transforma Geom al entorno de trabajo de Spheroid o Ellipsoid Datum, luego aplique la distancia a cada lado, obtenga el resultado y transforme para ingresar SRID.

Función ================================================= ==================

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;

Úselo con una consulta simple.

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

Resultado ================================================= ==================ingrese la descripción de la imagen aquí

  1. Función basada en generador en serie.

Función ================================================= =================

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;

Úselo con una consulta simple.

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

ingrese la descripción de la imagen aquí

Muhammad Imran Siddique
fuente
3

Entonces 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
JaakL
fuente
1
Hola, obtengo un resultado vacío con la función makegrid. El archivo de forma se importó a PostGIS usando shp2pgsql. No tengo idea de qué podría estar causando problemas, el srs está configurado en wgs84.
Michal Zimmermann
3

Aquí hay otro enfoque que sin duda es más rápido y 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 un ráster vacío, luego convierte cada píxel en un punto. Ventaja: no tenemos que verificar nuevamente si el polígono original interseca los puntos.

Opcional:

También puede agregar la alineación de 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)
obchardon
fuente
1

Una pequeña actualización potencial de las respuestas anteriores: tercer argumento como escala para wgs84 (o use 1 para las normales), y también redondeando 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
Martín
fuente
¿No sería mejor transformar la geometría en un SRID específico (como EPSG: 3857) que simplemente multiplicar por un factor de escala?
Nikolaus Krismer