Wie erstelle ich ein regelmäßiges Punktraster in einem Polygon in Postgis?

30

Wie kann man innerhalb eines Polygons ein regelmäßiges Punktraster mit den Abständen x, y in postgis erstellen? Wie das Beispiel:

Alt-Text

Pablo
quelle
Ich habe versucht, Schnittpolygone zu erstellen, die diesen Code mit dem Würfelcode "postGIS in Aktion" zusammenführen. Es wurde jedoch nur ein Polygon erstellt. Habe ich etwas vergessen? CREATE OR REPLACE FUNCTION makegrid (Geometrie, Ganzzahl, Ganzzahl) RETURNS Geometrie 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_Point (x, y), $ 3), st_setsrid (ST_Point (x + $ 2, y + $ 2), $ 3)) als Geometrie), $ 3) als Geom2 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 (
aurel_nc
Schau dir meine ausführliche Antwort an.
Muhammad Imran Siddique

Antworten:

29

Das machen Sie mit generate_series.

Wenn Sie nicht manuell schreiben möchten, wo das Raster beginnen und enden soll, erstellen Sie am einfachsten eine Funktion.

Ich habe das Folgende nicht richtig getestet, aber ich denke, es sollte funktionieren:

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

Um es zu benutzen, können Sie folgendes tun:

SELECT makegrid(the_geom, 1000) from mytable;

Dabei ist das erste Argument das Polygon, in dem Sie das Gitter haben möchten, und das zweite Argument der Abstand zwischen den Punkten im Gitter.

Wenn Sie einen Punkt pro Zeile möchten, verwenden Sie einfach ST_Dump wie folgt:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

HTH

Nicklas

Nicklas Avén
quelle
1
Möglicherweise müssen Sie st_setSRID () zu den st_point-Funktionen hinzufügen, andernfalls funktioniert st_intersects nicht.
JaakL
Meine getestete Version wurde als separate Antwort hinzugefügt.
JaakL
12

Ich habe den makegrid-Funktionscode von Nicklas Avén aufgegriffen und ihn durch Lesen und Verwenden des srid aus der Polygongeometrie etwas allgemeiner gestaltet. Andernfalls würde die Verwendung eines Polygons mit einem definierten srid einen Fehler ergeben.

Die Funktion:

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

Die Verwendung der Funktion erfolgt genau wie von Nicklas Avén geschrieben:

SELECT makegrid(the_geom, 1000) from mytable;

oder wenn Sie einen Punkt pro Zeile wollen:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

Hoffe, dass dies für jemanden nützlich sein wird.

Alex

Alexandre Neto
quelle
Die akzeptierte Antwort funktioniert aufgrund von SRID-Fehlern nicht mit meinen Daten. Diese Änderung funktioniert besser.
Vitaly Isaev
Vielleicht möchten Sie etwas hinzufügen, wenn das Polygon den Antimeridian überschreitet? Ich kann mir vorstellen, dass dies ein Problem mit xmin / xmax verursachen würde.
Thomas
2
Es hat bei mir nicht funktioniert. Verwenden von Postgres 9.6 und PostGIS 2.3.3. Im Aufruf generate_series musste ich dies als zweiten Parameter "ceiling (st_xmax ($ 1)) :: int" anstelle von "ceiling (st_xmax ($ 1) -st_xmin ($ 1)) :: int" und "ceiling ( st_ymax ($ 1)) :: int "anstelle von" ceiling (st_ymax ($ 1) -st_ymin ($ 1)) :: int "
Vitor Sapucaia
Ich stimme dem vorherigen Kommentar zu. Die Obergrenze von generate_series sollte die Obergrenze des Maximums und nicht die Obergrenze der Differenz (Maximum - Minimum) sein.
R. Bourgeon
10

Personen, die eine wgs84-Geometrie verwenden, werden wahrscheinlich Probleme mit dieser Funktion haben, da

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 

Nur ganze Zahlen zurückgeben. Mit Ausnahme von sehr großen Geometrien wie Ländern (die auf mehreren Längen- und Breitengraden liegen) wird hierdurch nur ein Punkt erfasst, der die Geometrie in den meisten Fällen nicht einmal selbst schneidet ... => leeres Ergebnis!

Mein Problem war, dass ich bei Gleitkommazahlen wie WSG84 scheinbar keine generate_series () mit einem Dezimalabstand verwenden kann ... Deshalb habe ich die Funktion so optimiert, dass sie trotzdem funktioniert:

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)))

Grundsätzlich genau das gleiche. Multiplizieren und dividieren Sie einfach mit 1000000, um die Dezimalstellen ins Spiel zu bringen, wenn ich sie brauche.

Dafür gibt es sicherlich eine bessere Lösung. ++

Julien Garcia
quelle
Das ist eine clevere Lösung. Haben Sie die Ergebnisse überprüft? Sind sie konsistent?
Pablo
Hallo. Ja, Pablo. Ich bin mit den bisherigen Ergebnissen zufrieden. Ich brauchte es, um ein Polygon mit relativer Höhe über dem Boden zu bauen. (Ich verwende SRTM, um die gewünschte Höhe für jeden Gitterpunkt zu berechnen.) Mir fehlt jetzt nur noch eine Möglichkeit, die Punkte einzubeziehen, die auch auf dem Umfang des Polygons liegen. Derzeit ist die gerenderte Form am Rand etwas abgeschnitten.
Julien Garcia
hat funktioniert, als alle anderen Lösungen fehlgeschlagen sind, danke!
Jordan Arseno
7

Dieser Algorithmus sollte in Ordnung sein:

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);
       }
    }
}

Dabei ist "Polygon" das Polygon und "Auflösung" die erforderliche Rasterauflösung.

Zur Implementierung in PostGIS sind möglicherweise folgende Funktionen erforderlich:

Viel Glück!

julien
quelle
1
Beachten Sie, dass dieser Ansatz nicht optimal ist, wenn Sie große komplexe Polygonbereiche haben (z. B. Küstenpuffer).
JaakL
Was schlagen Sie stattdessen vor?
Juli
4

Drei Algorithmen mit unterschiedlichen Methoden.

Github Repo Link

  1. Der einfache und beste Ansatz unter Verwendung der tatsächlichen Erdentfernung der Koordinaten von der x- und y-Richtung. Der Algorithmus funktioniert mit jeder SRID, intern arbeitet er mit WGS 1984 (EPSG: 4326) und führt die Rücktransformation zur Eingabe der SRID durch.

Funktion =============================================== =================

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;

Verwenden Sie die Funktion mit einer einfachen Abfrage. Die Geometrie muss gültig und vom Typ Polygon, Multipolygon oder Umschlag sein

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

Ergebnis =============================================== ====================

Bildbeschreibung hier eingeben

  1. Zweite Funktion basierend auf dem Nicklas Avén- Algorithmus. Ich habe es verbessert, um mit jeder SRID umzugehen.

    habe folgende änderungen im algorithmus aktualisiert.

    1. Separate Variable für die x- und y-Richtung für die Pixelgröße,
    2. Neue Variable zur Berechnung des Abstands in Sphäroid oder Ellipsoid.
    3. Geben Sie eine beliebige SRID ein, transformieren Sie Geom in die Arbeitsumgebung von Spheroid oder Ellipsoid Datum, wenden Sie dann den Abstand zu jeder Seite an, erhalten Sie das Ergebnis und transformieren Sie die SRID.

Funktion =============================================== =================

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;

Verwenden Sie es mit einer einfachen Abfrage.

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

Ergebnis ====================================== =================Bildbeschreibung hier eingeben

  1. Funktion basiert auf Seriengenerator.

Funktion ====================================== =================

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;

Verwenden Sie es mit einer einfachen Abfrage.

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

Bildbeschreibung hier eingeben

Muhammad Imran Siddique
quelle
3

Hier ist ein weiterer Ansatz, der sicherlich schneller und verständlicher ist.

Zum Beispiel für ein 1000 x 1000 m-Raster:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0))).geom 
FROM the_polygon

Auch die ursprüngliche SRID bleibt erhalten.

Dieses Snippet konvertiert die Geometrie des Polygons in ein leeres Raster und konvertiert dann jedes Pixel in einen Punkt. Vorteil: Wir müssen nicht noch einmal prüfen, ob das ursprüngliche Polygon die Punkte schneidet.

Wahlweise:

Sie können die Rasterausrichtung auch mit den Parametern gridx und gridy hinzufügen. Da wir jedoch den Schwerpunkt jedes Pixels (und nicht eine Ecke) verwenden, müssen wir ein Modulo verwenden, um den richtigen Wert zu berechnen:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0,mod(1000/2,100),mod(1000/2,100)))).geom 
FROM the_polygon

Mit mod(grid_size::integer/2,grid_precision)

Hier ist die Postgres-Funktion:

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;

Kann verwendet werden mit:

SELECT makegrid(the_geom,1000.0,100) as geom from the_polygon  
-- makegrid(the_geom,grid_size,alignement)
Obchardon
quelle
2

Also meine feste Version:

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

Verwendung:

SELECT (ST_Dump(makegrid(the_geom, 1000, 3857))).geom as the_geom from my_polygon_table
JaakL
quelle
1
Hi, ich bekomme leeres Ergebnis mit makegrid Funktion. Das Shapefile wurde mit shp2pgsql in PostGIS importiert. Habe keine Ahnung, was Probleme verursachen könnte, der srs ist auf wgs84 eingestellt.
Michal Zimmermann
1

Ein kleines mögliches Update der vorherigen Antworten - drittes Argument als Maßstab für wgs84 (oder verwenden Sie 1 für normale) und auch das Runden innerhalb des Codes, damit die skalierten Punkte auf mehreren Formen ausgerichtet werden.

Hoffe das hilft, 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
Martin
quelle
Wäre es nicht besser, die Geometrie in eine bestimmte SRID (wie EPSG: 3857) zu transformieren, als nur mit einem Skalierungsfaktor zu multiplizieren?
Nikolaus Krismer