Erzeugen Sie Punkte, die innerhalb eines Polygons liegen

30

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.

user2024
quelle
3
Im Gegenteil, die Zeit ist vorhersehbar. Es ist proportional zum Verhältnis der Ausdehnungsfläche des Polygons geteilt durch die Fläche des Polygons, multipliziert mit der Zeit, die benötigt wird, um einen einzelnen Punkt zu generieren und zu testen. Die Zeit variiert ein wenig, aber die Variation ist proportional zur Quadratwurzel der Anzahl der Punkte. Bei größeren Zahlen spielt dies keine Rolle. Brechen Sie gewundene Polygone bei Bedarf in kompaktere Stücke, um dieses Flächenverhältnis auf einen niedrigen Wert zu reduzieren. Nur fraktale Polygone werden dir Probleme bereiten, aber ich bezweifle, dass du sie hast!
Whuber
@Pablo: gute Funde. Allerdings sind diese beiden Fragen Software spezifisch und betreffen beide platzieren regelmäßige Anordnungen von Punkten innerhalb von Polygonen, nicht zufällige Punkte
whuber
stimmen Sie zu, wenn der Unterschied zwischen zufälligen Punkten und der regelmäßigen Punkterzeugung innerhalb eines Polygons besteht.
Mapperz

Antworten:

20

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

Dan S.
quelle
2
+1 Einfach und effektiv. Es ist erwähnenswert, dass innerhalb eines Dreiecks gleichmäßig zufällige Punkte erzeugt werden können, ohne dass eine Zurückweisung erfolgt, da es (leicht zu berechnende) flächenerhaltende Zuordnungen zwischen einem Dreieck und einem gleichschenkligen Dreieck gibt, das ein halbes Quadrat ist, also die Hälfte wo Die y-Koordinate überschreitet die x-Koordinate. Generieren Sie zwei zufällige Koordinaten und sortieren Sie sie, um einen zufälligen Punkt im gleichschenkligen Dreieck zu erhalten. Ordnen Sie diesen dann wieder dem ursprünglichen Dreieck zu.
whuber
+1 Ich mag die Diskussion über trilineare Koordinaten, auf die in dem von Ihnen zitierten Artikel verwiesen wird, wirklich . Ich nehme an, dass dies einer Kugel zugänglich wäre, deren Oberfläche als eine Tesselation von Dreiecken dargestellt wird. In einem projizierten Flugzeug wäre es keine wirklich zufällige Verteilung, oder?
Kirk Kuykendall
@whuber - +1 zurück bei dir. Eine andere Möglichkeit (im Link, aber sie winken von Hand darüber) besteht darin, die zurückgewiesenen Punkte aus dem gleichmäßig abgetasteten Viereck über die gemeinsame Kante und zurück in das Dreieck zu reflektieren.
Dan S.
@Kirk - Der Link zum Zitieren ist insofern nicht hilfreich, als er eine Reihe falscher (nicht einheitlicher) Abtastmethoden, einschließlich trilinearer Koordinaten, vor dem "richtigen" Weg auflistet. Es sieht nicht so aus, als gäbe es einen direkten Weg, um eine einheitliche Abtastung mit trilinearen Koordinaten zu erhalten. Ich würde eine einheitliche Abtastung über die gesamte Sphäre erreichen, indem ich zufällige Einheitsvektoren in 3d in ihr Lat / Lon-Äquivalent umwandle, aber das bin nur ich. (Wir sind unsicher, ob die Abtastung auf kugelförmige Dreiecke / Polygone beschränkt ist.) (Wir sind uns auch nicht sicher, ob die Abtastung beispielsweise auf Wgs84 wirklich gleichmäßig ist. Ich denke, nur die Auswahl von Winkeln wirkt sich ein wenig auf die Pole aus.)
Dan S.
1
@Dan Verwenden Sie zum gleichmäßigen Abtasten der Kugel eine zylindrische Projektion mit gleicher Fläche (Koordinaten sind Längengrad und Kosinus des Breitengrads). Wenn Sie ohne Projektion abtasten möchten, gibt es einen schönen Trick: Generieren Sie drei unabhängige Standard-Normalvariablen (x, y, z) und projizieren Sie sie auf den Punkt (R x / n, R y / n, R * z / n ) wobei n ^ 2 = x ^ 2 + y ^ 2 + z ^ 2 und R der Erdradius ist. Bei Bedarf in (lat, lon) konvertieren (bei Arbeiten an einem Sphäroid authale Breiten verwenden). Es funktioniert, weil diese trivariate Normalverteilung sphärisch symmetrisch ist. Halten Sie sich zum Abtasten von Dreiecken an eine Projektion.
Whuber
14

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

Bildbeschreibung hier eingeben

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

Underdunkel
quelle
1
Auch 5 Jahre später immer noch sehr hilfreich!
Gestrandetes Kind
10

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;
}
Dan Walton
quelle
10

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
monkut
quelle
2
Sollte es nicht von shapely.geometry import sein ...?
PyMapr
1
Sie können auch die representative_pointMethode verwenden: shapely.readthedocs.io/en/latest/…
dain
6

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 .

mdsumner
quelle
+1 - Da R Open Source ist, wenn man replizieren möchte, kann man immer in die Quelle gehen, um zu sehen, wie sie gemacht werden. Für Punktmuster sind möglicherweise auch die Simulationswerkzeuge im spatstat-Paket von Interesse.
Andy W
5

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
}
whuber
quelle
2

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.

Algoseer
quelle
2

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.

Generieren zufälliger Punkte in ausgewählten Polygonen

markusN
quelle
1
Obligatorische Erinnerung, dass Postleitzahlen Linien und keine Polygone sind.
Richard
Können Sie näher darauf eingehen? Mir geht es auch hier um Bereiche: en.wikipedia.org/wiki/ZIP_Code#Primary_state_prefixes
markusN
Sicher: Postleitzahlen verweisen auf bestimmte Postämter und deren Postzustellungswege. Postleitzahlen sind daher Linien, keine Polygone. Sie können sich überlappen, Löcher enthalten und müssen nicht unbedingt die gesamten USA oder einen bestimmten Staat abdecken. Aus diesem Grund ist es gefährlich, sie zur Aufteilung von Flächen zu verwenden. Zensuseinheiten (wie Blockgruppen) sind eine bessere Wahl. Siehe auch: dies und das .
Richard
1
Vielen Dank! Es hängt wahrscheinlich auch vom Land ab, siehe zum Beispiel en.wikipedia.org/wiki/Postal_codes_in_Germany - Postleitzahlen sind jedoch nicht mein Kernthema, sondern wollten nur die ursprüngliche Frage "Punkte generieren, die innerhalb eines Polygons liegen" veranschaulichen und beantworten Diskutiere die Definitionen der Postleitzahl, die hier OT ist :-)
markusN
1
In beiden Punkten zur Kenntnis genommen. Ich sollte wahrscheinlich einen kleinen Blogeintrag machen, damit ich dies beim nächsten Mal genauer sagen kann :-)
Richard
1

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 =============================================== ==================== Bildbeschreibung hier eingeben

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