¿Crear una cuadrícula de polígonos regular en PostGIS?

56

¿Cómo crear, en la forma de un polígono, una cuadrícula regular de polígonos / cuadrados de un tamaño determinado, en postgis?

He pensado en una función como Cómo crear contenido regular Apunte la cuadrícula dentro de un polígono en Postgis? solo para los cuadrados, de modo que los cuadrados puedan ser de 5m x 5m o incluso 10m x 10m. Pero no tengo idea de cambiar esto de la manera correcta.

    
pregunta mk.archaeo 31.10.2011 - 20:41

4 respuestas

56

Aquí hay un conjunto que devuelve la función ST_CreateFishnet que crea una cuadrícula 2D de geometrías de polígonos:

CREATE OR REPLACE FUNCTION ST_CreateFishnet(
        nrow integer, ncol integer,
        xsize float8, ysize float8,
        x0 float8 DEFAULT 0, y0 float8 DEFAULT 0,
        OUT "row" integer, OUT col integer,
        OUT geom geometry)
    RETURNS SETOF record AS
$$
SELECT i + 1 AS row, j + 1 AS col, ST_Translate(cell, j * $3 + $5, i * $4 + $6) AS geom
FROM generate_series(0, $1 - 1) AS i,
     generate_series(0, $2 - 1) AS j,
(
SELECT ('POLYGON((0 0, 0 '||$4||', '||$3||' '||$4||', '||$3||' 0,0 0))')::geometry AS cell
) AS foo;
$$ LANGUAGE sql IMMUTABLE STRICT;

donde nrow y ncol son el número de filas y columnas, xsize y ysize son las longitudes del tamaño de celda, y opcionales x0 y y0 son coordenadas para la parte inferior izquierda esquina.

El resultado es row y col números, comenzando desde 1 en la esquina inferior izquierda y geom polígonos rectangulares para cada celda. Así por ejemplo:

SELECT *
FROM ST_CreateFishnet(4, 6, 10, 10) AS cells;
 row | col |         geom
-----+-----+--------------------------------
   1 |   1 | 0103000000010000000500000000...
   2 |   1 | 0103000000010000000500000000...
   3 |   1 | 0103000000010000000500000000...
   4 |   1 | 0103000000010000000500000000...
   1 |   2 | 0103000000010000000500000000...
   2 |   2 | 0103000000010000000500000000...
   ...
   3 |   6 | 0103000000010000000500000000...
   4 |   6 | 0103000000010000000500000000...
(24 rows)

O para hacer una colección de geometría única para la cuadrícula completa:

SELECT ST_Collect(cells.geom)
FROM ST_CreateFishnet(4, 6, 10, 10) AS cells;

Puede agregar las compensaciones de origen x0 / y0 (estas están predeterminadas a cero).

    
respondido por el Mike T 01.11.2011 - 02:45
50

Aquí hay una variante específica de generación, para una situación en la que necesita crear una cuadrícula para un mapa geográfico con un paso métrico constante (las celdas se pueden usar para agrupar valores, por ejemplo, densidad de rayos en una región).

La función no es muy elegante, pero no encontré ninguna solución mejor para esa tarea (incluida la función de Mike Toews arriba). Así que tienes un polígono encuadernado (por ejemplo, llegado desde una interfaz de Google Maps), tienes un valor de paso en metros:

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  grid_step integer,
  metric_srid integer = 28408 --metric SRID (this particular is optimal for the Western Russia)
)
RETURNS public.geometry AS
$body$
DECLARE
  BoundM public.geometry; --Bound polygon transformed to the metric projection (with metric_srid SRID)
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  sectors public.geometry[];
  i INTEGER;
BEGIN
  BoundM := ST_Transform($1, $3); --From WGS84 (SRID 4326) to the metric projection, to operate with step in meters
  Xmin := ST_XMin(BoundM);
  Xmax := ST_XMax(BoundM);
  Ymax := ST_YMax(BoundM);

  Y := ST_YMin(BoundM); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  --Better if generating polygons exceeds the bound for one step. You always can crop the result. But if not you may get not quite correct data for outbound polygons (e.g. if you calculate frequency per sector)
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      i := i + 1;
      sectors[i] := ST_GeomFromText('POLYGON(('||X||' '||Y||', '||(X+$2)||' '||Y||', '||(X+$2)||' '||(Y+$2)||', '||X||' '||(Y+$2)||', '||X||' '||Y||'))', $3);

      X := X + $2;
    END LOOP xloop;
    Y := Y + $2;
  END LOOP yloop;

  RETURN ST_Transform(ST_Collect(sectors), ST_SRID($1));
END;
$body$
LANGUAGE 'plpgsql';

Cómo usarlo:

SELECT cell FROM 
(SELECT (
ST_Dump(makegrid_2d(ST_GeomFromText('Polygon((35.099577 45.183417,47.283415 45.183417,47.283415 49.640445,35.099577 49.640445,35.099577 45.183417))',
 4326), -- WGS84 SRID
 10000) -- cell step in meters
)).geom AS cell) AS q_grid

Para que pueda ver que las líneas formateadas por polígonos generados se encuentran a lo largo de los paralelos y meridianos geográficos, eso es muy conveniente.

Consejo: Si calcula algo como densidad (p. Ej., Mapa de trazos de rayos por celdas), y la cuadrícula se genera dinámicamente Para aumentar el rendimiento, sugiero usar tablas temporales para almacenar celdas como polígonos de geometría, con una índice espacial en la columna representa la celda.

    
respondido por el Alexander Palamarchuk 16.04.2012 - 17:13
9

Puede crear una cuadrícula regular simplemente vectorizando un ráster vacío:

SELECT (ST_PixelAsPolygons(ST_AddBand(ST_MakeEmptyRaster(100, 100, 1.1, 1.1, 1.0), '8BSI'::text, 1, 0), 1, false)).geom
    
respondido por el Pierre Racine 06.07.2017 - 15:30
5

He creado una variante de la función @ Alexander que no requiere que nos transformemos en otro SRID. Esto evita el problema de tener que encontrar una proyección que utilice medidores como unidades para una región en particular. Utiliza ST_Project para caminar adecuadamente usando la proyección dada. También agregué width_step y height_step para permitir mosaicos rectangulares en lugar de requerir que sean cuadrados.

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  width_step integer,
  height_step integer
)
RETURNS public.geometry AS
$body$
DECLARE
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  NextX DOUBLE PRECISION;
  NextY DOUBLE PRECISION;
  CPoint public.geometry;
  sectors public.geometry[];
  i INTEGER;
  SRID INTEGER;
BEGIN
  Xmin := ST_XMin(bound_polygon);
  Xmax := ST_XMax(bound_polygon);
  Ymax := ST_YMax(bound_polygon);
  SRID := ST_SRID(bound_polygon);

  Y := ST_YMin(bound_polygon); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
      NextX := ST_X(ST_Project(CPoint, $2, radians(90))::geometry);
      NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);

      i := i + 1;
      sectors[i] := ST_MakeEnvelope(X, Y, NextX, NextY, SRID);

      X := NextX;
    END LOOP xloop;
    CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
    NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);
    Y := NextY;
  END LOOP yloop;

  RETURN ST_Collect(sectors);
END;
$body$
LANGUAGE 'plpgsql';

Puedes usarlo así:

SELECT ST_AsGeoJSON(cell) FROM (
  SELECT (
    ST_Dump(
      makegrid_2d(
        ST_GeomFromText(
          'Polygon((35.099577 45.183417,47.283415 45.183417,47.283415 49.640445,35.099577 49.640445,35.099577 45.183417))',
          4326
        ),
         10000, -- width step in meters
         10000  -- height step in meters
       ) 
    )
  ) .geom AS cell
)q;
    
respondido por el Matt 04.10.2017 - 17:47

Lea otras preguntas en las etiquetas