Generuj punkty, które leżą wewnątrz wielokąta

Mam funkcję wielokąta i chcę być w stanie generować w niej punkty. Potrzebuję tego do jednego zadania klasyfikacyjnego.

Generowanie losowych punktów, dopóki jeden nie znajdzie się w wielokącie, nie zadziałałoby, ponieważ to naprawdę nieprzewidywalny czas, jaki zajmuje.

Zacznij od rozłożenia wielokąta na trójkąty, a następnie wygeneruj punkty wewnątrz nich . (Aby uzyskać równomierny rozkład, zważ każdy trójkąt według jego powierzchni.)

Gdy umieścisz znacznik QGIS w tym pytaniu: narzędzia Losowe punkty można używać z warstwą graniczną.

Jeśli szukasz kodu, pomocny może być kod źródłowy wtyczki.

Możesz określić zasięg wielokąta, a następnie ograniczyć generowanie liczb losowych dla wartości X i Y w tych zakresach.

Podstawowy proces: 1) Określ maxx, maxy, minx, miny wierzchołków wielokąta, 2) Wygeneruj losowe punkty, używając tych wartości jako granic 3) Przetestuj każdy punkt pod kątem przecięcia z wielokątem, 4) Przestań generować, gdy masz wystarczającą liczbę punktów spełniających przecięcie test

Oto algorytm (C #) dla testu skrzyżowania:

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

Istnieje kilka dobrych bibliotek, które wykonują większość ciężkich zadań dla Ciebie.

Przykład użycia [foremnego] [1] w pythonie.

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)

Lub użyj, .representative_point()aby uzyskać punkt w obiekcie (jak wspomniano przez dain):

Zwraca tanio obliczony punkt, który z pewnością znajduje się w obiekcie geometrycznym.

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

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

Jeśli R jest opcją, patrz ?spsamplew sppakiecie. Wielokąty można odczytać z dowolnego formatu obsługiwanego przez GDAL wbudowanego w pakiet rgdal, a następnie spsampledziała bezpośrednio na importowanym obiekcie z różnymi opcjami próbkowania.

Chciałbym zaoferować rozwiązanie, które wymaga niewiele w zakresie analizy GIS. W szczególności nie wymaga triangulacji żadnych wielokątów.

Poniższy algorytm, podany w pseudokodzie, odnosi się do kilku prostych operacji oprócz podstawowych możliwości obsługi list (tworzenie, znajdowanie długości, dołączanie, sortowanie, wyodrębnianie podlist i konkatenację) oraz generowanie losowych liczb zmiennoprzecinkowych w przedziale [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.

Wszystkie są dostępne w prawie każdym środowisku programowania GIS lub graficznym (i łatwe do kodowania, jeśli nie). Clipnie mogą zwracać zdegenerowanych wielokątów (tzn. tych o zerowym obszarze).

Procedura SimpleRandomSampleskutecznie uzyskuje listę punktów losowo rozmieszczonych w wielokącie. Jest to opakowanie dlaSRS , która rozbija wielokąt na mniejsze kawałki, aż każdy kawałek będzie wystarczająco zwarty, aby można było z niego efektywnie pobrać próbki. Aby to zrobić, wykorzystuje wstępnie obliczoną listę liczb losowych, aby zdecydować, ile punktów należy przeznaczyć na każdy kawałek.

SRS można „dostroić”, zmieniając parametr t. Jest to maksymalny stosunek ramki granicznej: pola powierzchni wielokąta, który może być tolerowany. Zmniejszenie (ale większe niż 1) spowoduje, że większość wielokątów zostanie podzielona na wiele części; jej powiększenie może spowodować odrzucenie wielu punktów próbnych dla niektórych wielokątów (faliste, z fragmentami lub pełne dziur). Gwarantuje to, że maksymalny czas na próbkowanie oryginalnego wielokąta jest przewidywalny.

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

Następna procedura wywołuje się rekurencyjnie, jeśli to konieczne. Tajemnicze wyrażenie t*N + 5*Sqrt(t*N)konserwatywnie szacuje górną granicę liczby punktów, które będą potrzebne, uwzględniając zmienność szans. Prawdopodobieństwo, że to się nie powiedzie, wynosi tylko 0,3 na milion wywołań procedur. Zwiększ od 5 do 6 lub nawet 7, aby zmniejszyć to prawdopodobieństwo, jeśli chcesz.

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
}

Jeśli twój wielokąt jest wypukły i znasz wszystkie wierzchołki, możesz rozważyć wykonanie „losowego” wypukłego ważenia wierzchołków, aby wypróbować nowy punkt, który z pewnością będzie leżał wewnątrz wypukłego kadłuba (w tym przypadku wielokąta).

Powiedzmy na przykład, że masz wielokąt wypukły N z wierzchołkami

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

Następnie wygeneruj losowo N wypukłych wag

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

Losowo próbkowany punkt jest następnie podawany przez

Y= ∑ w_i*V_i

Może być inny sposób próbkowania N wypukłych wag

• Wybieraj liczby N-1 równomiernie losowo w obrębie zakresu (bez zamiany), sortuj je i znormalizuj przedziały N między nimi, aby uzyskać wagi.
• Możesz także próbkować z rozkładu Dirichleta, który jest często używany jako koniugat przed rozkładem wielomianowym, który jest podobny do wypukłych wag w twoim przypadku.

Kiedy twój wielokąt nie jest bardzo niezbyt wypukły, możesz najpierw rozważyć przekształcenie go w wypukły kadłub. Powinno to przynajmniej w dużym stopniu ograniczyć liczbę punktów leżących poza twoim wielokątem.

Zadanie jest bardzo łatwe do rozwiązania w GRASS GIS (jedno polecenie) za pomocą v.random .

Poniżej przykład, jak dodać 3 losowe punkty do wybranych wielokątów (tutaj obszary kodu pocztowego miasta Raleigh, NC) ze strony podręcznika. Zmieniając instrukcję SQL „where”, można wybrać wielokąt (y).

Link do odpowiedzi

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

Trzy algorytmy wykorzystujące różne podejścia.

Git Repo Link

1. Oto proste i najlepsze podejście, wykorzystujące rzeczywistą odległość współrzędnych od kierunku xiy. Algorytm wewnętrzny wykorzystuje WGS 1984 (4326) i przekształca wynik na wstawiony SRID.

Funkcja ================================================= ==================

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;

Użyj funkcji z prostym zapytaniem, geometria musi być poprawna i wielokąta, wielokąta lub obwiedni

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

Wynik ================================================= =====================

2. Druga funkcja oparta na Nicklas Avén algorytmie . Próbowałem obsłużyć dowolny identyfikator SRID.

   Zastosowałem następujące zmiany w algorytmie.

   1. Oddzielna zmienna dla kierunku xiy dla rozmiaru piksela,
   2. Nowa zmienna do obliczania odległości w sferoidie lub elipsoidzie.
   3. Wprowadź dowolny SRID, przekształć funkcję Geom w środowisko pracy sferoidy lub elipsoidy Datum, następnie zastosuj odległość z każdej strony, uzyskaj wynik i przekształć na wejściowy SRID.

Funkcja ================================================= ==================

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;

Użyj go za pomocą prostego zapytania.

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

Wynik ================================================= ==================

3. Funkcja oparta na generatorze szeregowym.

Funkcja ================================================= =================

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;

Użyj go za pomocą prostego zapytania.

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