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

30

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.

użytkownik2024
źródło
3
Przeciwnie, czas jest przewidywalny. Jest proporcjonalny do stosunku powierzchni zasięgu wielokąta podzielonej przez powierzchnię wielokąta, razy czas potrzebny do wygenerowania i przetestowania pojedynczego punktu. Czas nieco się różni, ale zmiana jest proporcjonalna do pierwiastka kwadratowego z liczby punktów. W przypadku dużych liczb staje się to nieistotne. W razie potrzeby połam kręte wielokąty na bardziej zwarte kawałki, aby zmniejszyć ten współczynnik powierzchni do niskiej wartości. Tylko fraktalne wielokąty sprawią ci kłopoty, ale wątpię, czy je masz!
whuber
@Pablo: dobre znaleziska. Jednak oba te pytania są specyficzne dla oprogramowania i oba dotyczą umieszczania regularnych tablic punktów w wielokątach, a nie losowych punktów
whuber
zgadzają się z różnicą whuber to losowe punkty vs regularne generowanie punktów w obrębie wielokąta.
Mapperz

Odpowiedzi:

20

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

Dan S.
źródło
2
+1 Proste i skuteczne. Warto zauważyć, że jednolicie losowe punkty mogą być generowane w obrębie trójkąta bez żadnego odrzucenia, ponieważ istnieją (łatwo obliczalne) mapowania zachowujące obszar między dowolnym trójkątem a prawym trójkątem równoramiennym, który jest pół kwadratu, powiedz połowę gdzie współrzędna y przekracza współrzędną x. Wygeneruj dwie losowe współrzędne i posortuj je, aby uzyskać losowy punkt w trójkącie równoramiennym, a następnie zamapuj go z powrotem do oryginalnego trójkąta.
whuber
+1 Naprawdę podoba mi się dyskusja o współrzędnych trójliniowych, do których odwołuje się cytowany artykuł. Podejrzewam, że byłoby to podatne na kulę, której powierzchnia jest reprezentowana przez mozaikę trójkątów. W rzutowanym samolocie nie byłby to naprawdę losowy rozkład, prawda?
Kirk Kuykendall
@ whuber - +1 z powrotem na ciebie. Innym sposobem (w łączu, ale machają nad nim ręcznie) jest odzwierciedlenie odrzuconych punktów z czworoboku o jednakowej próbce na wspólnej krawędzi i z powrotem w trójkąt.
Dan S.
@Kirk - link cytowania jest dotykowym anty-pomocnym, ponieważ zawiera listę niewłaściwych (nierównomiernych) metod próbkowania, w tym współrzędnych trójliniowych, przed „właściwą” drogą. Nie wygląda na to, że istnieje bezpośredni sposób uzyskania jednolitego próbkowania z współrzędnymi trójliniowymi. Podchodziłbym do jednolitego próbkowania na całej kuli, konwertując losowe wektory jednostek w 3d na ich odpowiednik lat / lon, ale to tylko ja. (Nie jestem pewien co do próbkowania ograniczonego do sferycznych trójkątów / wielokątów.) (Również nie jestem pewien co do naprawdę jednolitego próbkowania na np. Wgs84: myślę, że tylko wybieranie kątów będzie nieco przesunięte w kierunku biegunów.)
Dan S.
1
@Dan Aby równomiernie próbkować kulę, użyj cylindrycznego rzutu o równej powierzchni (współrzędne to długość i cosinus szerokości). Jeśli chcesz próbkować bez użycia projekcji, istnieje piękna sztuczka: wygeneruj trzy niezależne standardowe zmienne normalne (x, y, z) i rzutuj je do punktu (R x / n, R y / n, R * z / n ) gdzie n ^ 2 = x ^ 2 + y ^ 2 + z ^ 2, a R jest promieniem ziemi. Konwertuj na (lat, lon), jeśli to konieczne (przy użyciu szerokości geograficznej podczas pracy na sferoidie). Działa, ponieważ ten trójwartościowy rozkład normalny jest sferycznie symetryczny. W przypadku próbkowania trójkątów trzymaj się występu.
whuber
14

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

wprowadź opis zdjęcia tutaj

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

podmrok
źródło
1
Nawet 5 lat później, nadal bardzo pomocny!
Stranded Kid
10

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;
}
Dan Walton
źródło
10

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
monkut
źródło
2
Czy nie powinno to pochodzić z importu shapely.geometry ...?
PyMapr
1
Możesz także użyć representative_pointmetody: shapely.readthedocs.io/en/latest/…
Dain
6

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.

mdsumner
źródło
+1 - Ponieważ R jest open source, jeśli ktoś chce się replikować, zawsze możesz przejść do źródła, aby zobaczyć, jak się to robi. W przypadku wzorów punktowych można również zainteresować się narzędziami do symulacji w pakiecie spatstat.
Andy W
5

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
}
Whuber
źródło
2

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.

algoseer
źródło
2

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

Generowanie losowych punktów w wybranych wielokątach

markusN
źródło
1
Obowiązkowe przypomnienie, że kody pocztowe to linie, a nie wielokąty.
Richard
Czy możesz rozwinąć? Dla mnie również tutaj odnosi się do obszarów: en.wikipedia.org/wiki/ZIP_Code#Primary_state_prefixes
markusN
Pewnie: kody pocztowe odnoszą się do poszczególnych urzędów pocztowych i ich dróg dostarczania poczty. W rezultacie kody pocztowe to linie, a nie wielokąty. Mogą nakładać się na siebie, zawierać dziury i niekoniecznie obejmować całe USA lub dany stan. Używanie ich do podziału obszaru jest niebezpieczne z tego powodu. Jednostki spisowe (jak grupy bloków) są lepszym wyborem. Zobacz także: to i to .
Richard
1
Dzięki! Prawdopodobnie zależy to również od kraju, patrz np. En.wikipedia.org/wiki/Postal_codes_in_Germany - jednak kody pocztowe nie są moim głównym tematem, chciałem tylko zilustrować i odpowiedzieć na pierwotne pytanie „Generuj punkty, które leżą wewnątrz wielokąta”, a nie przedyskutuj definicje kodów pocztowych, które są OT tutaj :-)
markusN
1
Zwrócono uwagę na obie liczby. Powinienem chyba zrobić mały wpis na blogu, aby następnym razem bardziej zwięźle to powiedzieć :-)
Richard
1

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 ================================================= ===================== wprowadź opis zdjęcia tutaj

  1. 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 ================================================= ==================wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

Muhammad Imran Siddique
źródło