Regelmäßiges Polygonraster in PostGIS erstellen?

61

Wie kann man in der Form eines Polygons ein regelmäßiges Gitter aus Polygonen / Quadraten einer bestimmten Größe in Postgis erstellen?

Ich habe mir überlegt, wie man in Postgis ein regelmäßiges Punktraster innerhalb eines Polygons erstellt. Nur für Quadrate, sodass Quadrate 5 x 5 m oder sogar 10 x 10 m groß sein können. Habe aber keine Ahnung, das richtig zu ändern.

mk.archaeo
quelle
2
Die Verallgemeinerung, die Sie suchen, ist nicht klar. Wollen Sie damit sagen, dass Sie mit einem (beliebigen) einzelnen Polygon beginnen und das Flugzeug mit kongruenten Kopien davon kacheln möchten? Im Allgemeinen ist dies nicht möglich, aber möglicherweise hat dieses Polygon bestimmte Eigenschaften (möglicherweise ist es beispielsweise als Parallelogramm, Dreieck oder Sechseck bekannt).
whuber

Antworten:

60

Hier ist eine Set-Return-Funktion ST_CreateFishnet, die ein 2D-Gitter aus Polygongeometrien erstellt:

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;

wo nrowund ncoldie Anzahl der Zeilen und Spalten, xsizeund ysizedie Längen der Zellengröße, und optional x0und y0sind Koordinaten der unteren linken Ecke.

Das Ergebnis ist rowund colZahlen, beginnend mit 1 in der unteren linken Ecke, und geomrechteckige Polygone für jede Zelle. Also zum Beispiel:

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)

Oder erstellen Sie eine einzelne Geometriesammlung für das gesamte Raster:

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

4x6 Gitter

Sie können die x0/ y0origin-Offsets hinzufügen (diese sind standardmäßig auf Null gesetzt).

Mike T
quelle
1
Vielen Dank! Jetzt muss ich nur noch das Netz an die BBox des Polygons binden.
mk.archaeo
Das ist sehr hilfreich. Ich habe eine Frage. Wie kann ich Gitter in einem Polygon / einer Bbox erstellen?
Mohammed Shafeek
Gute Arbeit Mike, das ist sehr hilfreich.
Mounaim,
56

Hier ist eine spezielle Variante der Generierung für eine Situation, in der Sie ein Raster für eine geografische Karte mit einem konstanten Metrikschritt erstellen müssen (Zellen können zum Gruppieren von Werten verwendet werden, z. B. Blitzdichte in einer Region).

Die Funktion ist nicht sehr elegant, aber ich habe keine bessere Lösung für diese Aufgabe gefunden (einschließlich der obigen Funktion von Mike Toews). Sie haben also ein gebundenes Polygon (z. B. von einer Google Maps-Oberfläche angekommen) und einen Schrittwert in Metern:

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';

Wie man es benutzt:

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

Sie sehen also, dass mit generierten Polygonen formatierte Linien entlang der geografischen Parallelen und Meridiane liegen - das ist sehr praktisch.

Ein Beispiel eines Rasters mit 50 km Schritt

Hinweis: Wenn Sie so etwas wie Dichte berechnen (z. B. Karte von Blitzschlägen nach Zellen) und das Raster dynamisch generiert wird, würde ich zur Leistungssteigerung die Verwendung von temporären Tabellen zum Speichern von Zellen als Geometriepolygone empfehlen, wobei ein räumlicher Index für die Spalte steht die Zelle.

Alexander Palamarchuk
quelle
Ich wünschte, ich könnte das noch einmal hochstimmen ... das war eine perfekte Lösung! und die Fähigkeit, das Koordinatensystem anzupassen, ist fantastisch ~!
DPSSpatial
Nur ein kleiner Vorschlag, anstatt ST_GeomFromTextbeim Erstellen einer Box zu verwenden sectors, können Sie ST_MakeEnvelopedie unteren linken und oberen rechten Koordinaten der Box verwenden und angeben.
Matt
Dies bringt Potenziale
nickves
11

Sie können ein reguläres Raster erstellen, indem Sie einfach ein leeres Raster vektorisieren:

SELECT (ST_PixelAsPolygons(ST_AddBand(ST_MakeEmptyRaster(100, 100, 1.1, 1.1, 1.0), '8BSI'::text, 1, 0), 1, false)).geom
Pierre Racine
quelle
1
Das ist eine so einfache Lösung, die es so oft auf vektorielle Weise getan hat.
John Powell
6

Ich habe eine Variante der @ Alexander-Funktion erstellt, die nicht erfordert, dass wir in eine andere SRID transformieren. Dies vermeidet die Notwendigkeit, eine Projektion zu finden, die Meter als Einheiten für eine bestimmte Region verwendet. Es wird verwendet, ST_Projectum mit der angegebenen Projektion richtig voranzukommen. Ich habe auch ein width_stepund hinzugefügt height_step, um rechteckige Kacheln zuzulassen, anstatt dass sie Quadrate sein müssen.

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';

Du kannst es so benutzen:

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;
Matt
quelle
5

Hier finden Sie einen optimierten und effizienten Algorithmus zum Erstellen von Fischernetzen, regelmäßigen Gitternetzen, Polygongitternetzen und rechteckigen Gitternetzen in beliebigen Umschlägen, Polygonen oder Multipolygonen. kann fast jede SRID verarbeiten;

GitHub Repo Link

Bildbeschreibung hier eingeben

DROP FUNCTION IF EXISTS PUBLIC.I_Grid_Regular(geometry, float8, float8);
CREATE OR REPLACE FUNCTION PUBLIC.I_Grid_Regular
( geom geometry, x_side float8, y_side float8, OUT geometry )
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;
geom_cell geometry := ST_GeomFromText(FORMAT('POLYGON((0 0, 0 %s, %s %s, %s 0,0 0))',
                                        $3, $2, $3, $2), srid);
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;
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 With foo AS (
    SELECT
    ST_Translate( geom_cell, j * $2 + x_min, i * $3 + y_min ) AS cell
    FROM
        generate_series ( 0, x_series ) AS j,
        generate_series ( 0, y_series ) AS i
    ) SELECT ST_CollectionExtract(ST_Collect(ST_Transform ( ST_Intersection(cell, geom), input_srid)), 3)
    FROM foo where ST_intersects (cell, geom);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Verwenden Sie es mit einer einfachen Abfrage. Die Eingabe sollte ein gültiges Polygon, Multipolygon oder Hüllkurve sein.

select I_Grid_Regular(st_setsrid(g.geom, 4326), .0001, .0001 ), geom from polygons limit 1
Muhammad Imran Siddique
quelle