Hasilkan titik yang terletak di dalam poligon

Saya memiliki fitur poligon dan ingin dapat menghasilkan poin di dalamnya. Saya memerlukan ini untuk satu tugas klasifikasi.

Menghasilkan titik acak sampai ada di dalam poligon tidak akan berhasil karena itu benar-benar tidak dapat diprediksi waktu yang diperlukan.

Mulailah dengan mendekomposisi poligon menjadi segitiga, lalu hasilkan poin di dalamnya . (Untuk distribusi yang seragam, timbang setiap segitiga berdasarkan wilayahnya.)

Saat Anda memasang tag QGIS pada pertanyaan ini: Alat Poin Acak dapat digunakan dengan lapisan batas.

Jika Anda mencari kode, kode sumber plugin yang mendasarinya harus membantu.

Anda bisa menentukan tingkat poligon, lalu membatasi pembangkitan angka acak untuk nilai X dan Y dalam luasan tersebut.

Proses dasar: 1) Tentukan maxx, maxy, minx, miny dari simpul poligon, 2) Hasilkan titik acak menggunakan nilai-nilai ini sebagai batas 3) Uji setiap titik untuk persimpangan dengan poligon Anda, 4) Berhenti menghasilkan ketika Anda memiliki cukup poin memenuhi persimpangan uji

Berikut ini adalah algoritma (C #) untuk tes persimpangan:

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

Ada beberapa perpustakaan bagus di luar sana yang melakukan sebagian besar pekerjaan berat untuk Anda.

Contoh menggunakan [rupawan] [1] dalam 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)

Atau gunakan .representative_point()untuk mendapatkan titik di dalam objek (seperti yang disebutkan oleh dain):

Mengembalikan titik komputasi yang murah yang dijamin berada dalam objek geometrik.

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

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

Jika R adalah opsi, lihat ?spsampledalam sppaket. Poligon dapat dibaca dari format apa pun yang didukung GDAL yang dibangun ke dalam paket rgdal, dan kemudian spsamplebekerja langsung pada objek yang diimpor dengan berbagai opsi pengambilan sampel.

Saya ingin menawarkan solusi yang membutuhkan sangat sedikit dalam hal analisis GIS. Secara khusus, tidak memerlukan triangulasi poligon apa pun.

Algoritme berikut, diberikan dalam pseudocode, merujuk pada beberapa operasi sederhana selain kemampuan penanganan daftar dasar (membuat, menemukan panjang, menambahkan, mengurutkan, mengekstrak sublists, dan menyatukan) dan generasi float acak dalam interval [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.

Ini semua tersedia di hampir semua SIG atau lingkungan pemrograman grafis (dan mudah dikodekan jika tidak). Cliptidak boleh mengembalikan poligon yang merosot (yaitu, yang tidak memiliki luas).

Prosedur secara SimpleRandomSampleefisien mendapatkan daftar titik yang didistribusikan secara acak dalam poligon. Ini adalah pembungkus untuk SRS, yang memecah poligon menjadi potongan-potongan kecil sampai masing-masing potongan cukup kompak untuk disampel secara efisien. Untuk melakukan ini, ia menggunakan daftar angka acak yang telah dihitung untuk memutuskan berapa banyak poin yang dialokasikan untuk masing-masing bagian.

SRS dapat "disetel" dengan mengubah parameter t. Ini adalah kotak batas maksimum: rasio area poligon yang dapat ditoleransi. Menjadikannya kecil (tetapi lebih besar dari 1) akan menyebabkan sebagian besar poligon terpecah menjadi banyak; membuatnya besar dapat menyebabkan banyak poin percobaan ditolak untuk beberapa poligon (berliku-liku, dengan sliver, atau penuh lubang). Ini menjamin bahwa waktu maksimum untuk mengambil sampel poligon asli dapat diprediksi.

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

Prosedur selanjutnya menyebut dirinya secara rekursif jika perlu. Ekspresi misterius secara t*N + 5*Sqrt(t*N)konservatif memperkirakan batas atas berapa banyak poin yang akan dibutuhkan, yang memperhitungkan variabilitas kesempatan. Kemungkinan bahwa ini akan gagal hanya 0,3 per juta panggilan prosedur. Tingkatkan 5 hingga 6 atau bahkan 7 untuk mengurangi kemungkinan ini jika Anda mau.

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
}

Jika poligon Anda adalah cembung dan Anda tahu semua simpul, Anda mungkin ingin mempertimbangkan untuk melakukan pembobotan cembung "acak" dari simpul untuk sampel titik baru yang dijamin terletak di dalam cembung cembung (poligon dalam kasus ini).

Misalnya, Anda memiliki poligon cembung sisi N dengan simpul

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

Kemudian menghasilkan bobot N cembung secara acak

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

Titik sampel acak kemudian diberikan oleh

Y= ∑ w_i*V_i

Mungkin ada cara berbeda untuk sampel bobot N cembung

• Pilih nomor N-1 secara acak secara acak dalam suatu rentang (tanpa penggantian), urutkan dan normalkan interval N di antara mereka untuk mendapatkan bobot.
• Anda juga dapat mengambil sampel dari distribusi Dirichlet yang sering digunakan sebagai konjugasi sebelum untuk distribusi multinomial, yang mirip dengan bobot cembung dalam kasing Anda.

Ketika poligon Anda tidak terlalu non-cembung, Anda dapat mempertimbangkan untuk mengubahnya menjadi lambung cembung. Setidaknya ini harus membatasi jumlah titik yang berada di luar poligon Anda untuk sebagian besar.

Tugas ini sangat mudah diselesaikan dalam GRASS GIS (satu perintah) menggunakan v.random .

Di bawah ini adalah contoh tentang cara menambahkan 3 titik acak ke dalam poligon yang dipilih (di sini area kode ZIP kota Raleigh, NC) dari halaman manual. Dengan memodifikasi pernyataan SQL "where", poligon dapat dipilih.

Tautan jawaban

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

Tiga algoritma menggunakan pendekatan yang berbeda.

Git Repo Link

1. Berikut ini adalah pendekatan sederhana dan terbaik, menggunakan jarak aktual koordinat dari arah x dan y. Algoritma internal menggunakan WGS 1984 (4326) dan hasil mentransformasikan ke SRID yang dimasukkan.

Fungsi ================================================= ==================

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;

Gunakan fungsi dengan kueri sederhana, geometri harus valid dan poligon, multi-poligon, atau amplop

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

Hasil ================================================= =====================

2. Fungsi kedua berdasarkan algoritma Nicklas Avén . Sudah mencoba menangani SRID apa pun.

   Saya telah menerapkan perubahan berikut dalam algoritma.

   1. Variabel terpisah untuk arah x dan y untuk ukuran piksel,
   2. Variabel baru untuk menghitung jarak dalam spheroid atau ellipsoid.
   3. Masukkan SRID, fungsi ubah Geom ke lingkungan kerja Spheroid atau Ellipsoid Datum, lalu terapkan jarak ke setiap sisi, dapatkan hasilnya dan ubah ke input SRID.

Fungsi ================================================= ==================

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;

Gunakan dengan kueri sederhana.

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

Hasil ================================================= ==================

3. Fungsi berdasarkan generator seri.

Fungsi ================================================= =================

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;

Gunakan dengan kueri sederhana.

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