¿Cómo crear, dentro de un polígono, una cuadrícula regular de puntos espaciados x, y en postgis? Como el ejemplo:
¿Cómo crear, dentro de un polígono, una cuadrícula regular de puntos espaciados x, y en postgis? Como el ejemplo:
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:
¡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.
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
estos cambios los he aplicado en el algoritmo original.
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
Lea otras preguntas en las etiquetas postgis