Erzeugen Sie Punkte, die innerhalb eines Polygons liegen

Ich habe ein Polygon-Feature und möchte darin Punkte erzeugen können. Ich brauche das für eine Klassifizierungsaufgabe.

Zufällige Punkte zu generieren, bis man sich innerhalb des Polygons befindet, würde nicht funktionieren, da es wirklich unvorhersehbar ist, wie lange es dauert.

Zerlegen Sie das Polygon zunächst in Dreiecke und generieren Sie dann Punkte in diesen . (Für eine gleichmäßige Verteilung gewichten Sie jedes Dreieck nach seiner Fläche.)

Wenn Sie ein QGIS-Tag auf diese Frage setzen: Das Random Points-Tool kann mit einer Grenzfläche verwendet werden.

Wenn Sie nach Code suchen, sollte der zugrunde liegende Quellcode des Plugins hilfreich sein.

Sie können die Ausdehnung des Polygons bestimmen und dann die Zufallszahlengenerierung für X- und Y-Werte innerhalb dieser Ausdehnungen einschränken.

Grundlegende Vorgehensweise: 1) Bestimmen Sie maxx, maxy, minx, miny der Polygonscheitelpunkte. 2) Generieren Sie zufällige Punkte unter Verwendung dieser Werte als Grenzen. 3) Testen Sie jeden Punkt auf Schnittpunkte mit Ihrem Polygon Prüfung

Hier ist ein Algorithmus (C #) für den Schnittpunkttest:

bool PointIsInGeometry(PointCollection points, MapPoint point)
{
int i;
int j = points.Count - 1;
bool output = false;

for (i = 0; i < points.Count; i++)
{
    if (points[i].X < point.X && points[j].X >= point.X || points[j].X < point.X && points[i].X >= point.X)
    {
        if (points[i].Y + (point.X - points[i].X) / (points[j].X - points[i].X) * (points[j].Y - points[i].Y) < point.Y)
        {
            output = !output;
        }
    }
    j = i;
}
return output;
}

Es gibt einige gute Bibliotheken, die das meiste für Sie tun.

Beispiel mit [shapely] [1] in Python.

import random
from shapely.geometry import Polygon, Point

def get_random_point_in_polygon(poly):
     minx, miny, maxx, maxy = poly.bounds
     while True:
         p = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
         if poly.contains(p):
             return p

p = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
point_in_poly = get_random_point_in_polygon(mypoly)

Oder verwenden Sie .representative_point(), um einen Punkt innerhalb des Objekts zu erhalten (wie von dain erwähnt):

Gibt einen kostengünstig berechneten Punkt zurück, der sich garantiert innerhalb des geometrischen Objekts befindet.

poly.representative_point().wkt
'POINT (-1.5000000000000000 0.0000000000000000)'

  [1]: https://shapely.readthedocs.io

Wenn R eine Option ist, siehe ?spsampleim spPaket. Die Polygone können aus jedem von GDAL unterstützten Format eingelesen werden, das in das rgdal-Paket integriert ist, und können dann spsamplemit einer Vielzahl von Abtastoptionen direkt auf importierte Objekte angewendet werden .

Ich möchte eine Lösung anbieten, die sehr wenig GIS-Analyse erfordert. Insbesondere müssen keine Polygone trianguliert werden.

Der folgende Algorithmus, der im Pseudocode angegeben ist, bezieht sich auf einige einfache Operationen zusätzlich zu den grundlegenden Funktionen zur Listenverarbeitung (Erstellen, Suchen von Längen, Anhängen, Sortieren, Extrahieren von Unterlisten und Verketten) und Generieren von zufälligen Gleitkommazahlen im Intervall [0, 1):

Area:        Return the area of a polygon (0 for an empty polygon).
BoundingBox: Return the bounding box (extent) of a polygon.
Width:       Return the width of a rectangle.
Height:      Return the height of a rectangle.
Left:        Split a rectangle into two halves and return the left half.
Right:       ... returning the right half.
Top:         ... returning the top half.
Bottom:      ... returning the bottom half.
Clip:        Clip a polygon to a rectangle.
RandomPoint: Return a random point in a rectangle.
Search:      Search a sorted list for a target value.  Return the index  
             of the last element less than the target.
In:          Test whether a point is inside a polygon.

Diese sind alle in fast jeder GIS- oder Grafikprogrammierungsumgebung verfügbar (und wenn nicht, einfach zu codieren). Clipdarf keine degenerierten Polygone zurückgeben (dh solche mit einer Fläche von Null).

Die Prozedur SimpleRandomSampleerhält effizient eine Liste von Punkten, die zufällig innerhalb eines Polygons verteilt sind. Es ist eine Hülle für SRS, die das Polygon in kleinere Stücke zerlegt, bis jedes Stück kompakt genug ist, um effizient abgetastet zu werden. Dazu verwendet es eine vorberechnete Liste von Zufallszahlen, um zu entscheiden, wie viele Punkte jedem Stück zugewiesen werden sollen.

SRS kann durch Ändern des Parameters "abgestimmt" werden t. Dies ist das maximale Verhältnis von Begrenzungsrahmen zu Polygonfläche, das toleriert werden kann. Wenn Sie es klein (aber größer als 1) machen, werden die meisten Polygone in viele Teile geteilt. Wenn Sie es groß machen, werden möglicherweise viele Testpunkte für einige Polygone verworfen (gewunden, mit Splittern oder voller Löcher). Dies garantiert, dass die maximale Zeit zum Abtasten des ursprünglichen Polygons vorhersehbar ist.

Procedure SimpleRandomSample(P:Polygon, N:Integer) {
    U = Sorted list of N independent uniform values between 0 and 1
    Return SRS(P, BoundingBox(P), U)
}

Die nächste Prozedur ruft sich bei Bedarf rekursiv auf. Der mysteriöse Ausdruck t*N + 5*Sqrt(t*N)schätzt konservativ eine Obergrenze für die Anzahl der benötigten Punkte unter Berücksichtigung der Zufallsvariabilität. Die Wahrscheinlichkeit, dass dies fehlschlägt, beträgt nur 0,3 pro Million Prozeduraufrufe. Erhöhen Sie 5 bis 6 oder sogar 7, um diese Wahrscheinlichkeit zu verringern, wenn Sie möchten.

Procedure SRS(P:Polygon, B:Rectangle, U:List) {
    N = Length(U)
    If (N == 0) {Return empty list}
    aP = Area(P)
    If (aP <= 0) {
        Error("Cannot sample degenerate polygons.")
        Return empty list
    }
    t = 2
    If (aP*t < Area(B)) {
        # Cut P into pieces
        If (Width(B) > Height(B)) {
            B1 = Left(B); B2 = Right(B)
        } Else {
            B1 = Bottom(B); B2 = Top(B)
        }
        P1 = Clip(P, B1); P2 = Clip(P, B2)
        K = Search(U, Area(P1) / aP)
        V = Concatenate( SRS(P1, B1, U[1::K]), SRS(P2, B2, U[K+1::N]) )
    } Else {
        # Sample P
        V = empty list
        maxIter = t*N + 5*Sqrt(t*N)
        While(Length(V) < N and maxIter > 0) {
            Decrement maxIter
            Q = RandomPoint(B)
            If (Q In P) {Append Q to V}
        }
       If (Length(V) < N) {
            Error("Too many iterations.")
       }
    }
    Return V
}

Wenn Ihr Polygon konvex ist und Sie alle Scheitelpunkte kennen, sollten Sie eine "zufällige" konvexe Gewichtung der Scheitelpunkte in Betracht ziehen, um einen neuen Punkt abzutasten, der garantiert innerhalb der konvexen Hülle liegt (in diesem Fall Polygon).

Angenommen, Sie haben ein N-seitiges konvexes Polygon mit Eckpunkten

V_i, i={1,..,N}

Dann erzeugen Sie zufällig N konvexe Gewichte

 w_1,w_2,..,w_N such that ∑ w_i = 1; w_i>=0

Der zufällig ausgewählte Punkt ist dann gegeben durch

Y= ∑ w_i*V_i

Es gibt verschiedene Möglichkeiten, N konvexe Gewichte abzutasten

• Wählen Sie N-1-Zahlen innerhalb eines Bereichs gleichmäßig zufällig aus (ohne Ersetzung), sortieren Sie sie und normalisieren Sie die N-Intervalle zwischen ihnen, um die Gewichte zu erhalten.
• Sie können auch eine Stichprobe aus der Dirichlet-Verteilung ziehen, die häufig als Konjugat vor der multinomialen Verteilung verwendet wird, die den konvexen Gewichten in Ihrem Fall ähnlich ist.

Wenn Ihr Polygon nicht sehr stark nicht konvex ist, sollten Sie es zunächst in eine konvexe Hülle umwandeln. Dies sollte zumindest die Anzahl der Punkte, die außerhalb Ihres Polygons liegen, stark einschränken.

Die Aufgabe ist in GRASS GIS (ein Befehl) mit v.random sehr einfach zu lösen .

Im Folgenden finden Sie ein Beispiel zum Hinzufügen von 3 zufälligen Punkten zu ausgewählten Polygonen (hier Postleitzahlenbereiche der Stadt Raleigh, NC) auf der Handbuchseite. Durch Ändern der SQL-Anweisung "where" können die Polygone ausgewählt werden.

Link beantworten

https://gis.stackexchange.com/a/307204/103524

Drei Algorithmen mit unterschiedlichen Ansätzen.

Git Repo Link

1. Hier ist ein einfacher und bester Ansatz, bei dem der tatsächliche Abstand der Koordinaten von der x- und y-Richtung verwendet wird. Der interne Algorithmus verwendet die WGS 1984 (4326) und die Ergebnistransformation zur eingefügten SRID.

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 sein und ein Polygon, mehrere Polygone oder eine Hüllkurve

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

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

2. Zweite Funktion basierend auf dem Nicklas Avén- Algorithmus. Habe versucht mit irgendwelchen SRID umzugehen.

   Ich habe folgende Änderungen am Algorithmus vorgenommen.

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

3. 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 ====================================== ========================