Создание точек, которые лежат внутри многоугольника

30

У меня есть функция многоугольника, и я хочу иметь возможность создавать точки внутри нее. Мне нужно это для одной задачи классификации.

Генерация случайных точек до тех пор, пока один из них не окажется внутри многоугольника, не сработает, потому что это действительно непредсказуемо.

user2024
источник
3
Наоборот, время предсказуемо. Он пропорционален отношению площади экстента многоугольника, деленной на площадь многоугольника, ко времени, необходимому для генерации и проверки одной точки. Время немного меняется, но изменение пропорционально квадратному корню из числа точек. Для больших чисел это становится несущественным. Если необходимо, разбейте извилистые многоугольники на более компактные части, чтобы уменьшить это соотношение площади до низкого значения. Только фрактальные многоугольники доставят вам неприятности, но я сомневаюсь, что они у вас есть!
whuber
1
а также: gis.stackexchange.com/questions/4663/…
Пабло,
@ Пабло: хорошие находки. Однако оба эти вопроса относятся к конкретному программному обеспечению и касаются размещения регулярных массивов точек внутри многоугольников, а не случайных точек
whuber
согласен с тем, что разница между случайными точками и регулярной генерацией точек внутри многоугольника.
Mapperz

Ответы:

20

Начните с разложения многоугольника на треугольники, затем создайте точки внутри них . (Для равномерного распределения, взвесьте каждый треугольник по площади.)

Дан С.
источник
2
+1 Просто и эффективно. Стоит отметить, что равномерно случайные точки могут быть сгенерированы в пределах треугольника без отклонения вообще, потому что есть (легко вычисляемые) сохраняющие площадь отображения между любым треугольником и равнобедренным прямоугольным треугольником, который является половиной квадрата, скажем, половина где координата у превышает координату х. Создайте две случайные координаты и отсортируйте их, чтобы получить случайную точку в равнобедренном треугольнике, затем сопоставьте ее с исходным треугольником.
whuber
+1 Мне очень нравится обсуждение трилинейных координат, на которые ссылается статья, на которую вы ссылаетесь. Я полагаю, что это поддается сфере, поверхность которой представлена ​​в виде тесселяции треугольников. На спроецированной плоскости это не было бы действительно случайным распределением, не так ли?
Кирк Куйкендалл
@whuber - +1 назад на тебя. Другой способ (в ссылке, но они помахали рукой над ней) состоит в том, чтобы отразить отклоненные точки из четырехугольника с одинаковой выборкой через общий край и обратно в треугольник.
Дэн С.
@Kirk - ссылка на цитирование является сенсорной анти-полезной в том смысле, что она перечисляет кучу неправильных (неоднородных) методов выборки, включая трилинейные координаты, перед «правильным» способом. Не похоже, что есть прямой способ получить равномерную выборку с трилинейными координатами. Я бы подошел к равномерной выборке по всей сфере путем преобразования случайных единичных векторов в 3d в их эквивалент широту / долготу, но это только я. (Не уверен в том, что выборка ограничена сферическими треугольниками / многоугольниками.) (Также не уверен насчет действительно равномерной выборки, например, на wgs84: просто выбор углов будет немного смещен к полюсам, я думаю.)
Дэн С.
1
@Dan Для равномерной выборки сферы используйте цилиндрическую проекцию равной площади (координаты - долгота и косинус широты). Если вы хотите произвести выборку без использования проекции, есть красивый трюк: сгенерируйте три независимых стандартных нормальных переменной (x, y, z) и спроецируйте их в точку (R x / n, R y / n, R * z / n ) где n ^ 2 = x ^ 2 + y ^ 2 + z ^ 2 и R - радиус Земли. При необходимости преобразуйте в (широта, долгота) (используя широты при работе на сфероиде). Это работает, потому что это тривариатное нормальное распределение сферически симметрично. Для выборки треугольников придерживайтесь проекции.
whuber
14

Если вы поставите тег QGIS на этот вопрос: инструмент «Случайные точки» можно использовать с пограничным слоем.

введите описание изображения здесь

Если вы ищете код, исходный код плагина должен помочь.

Подземье
источник
1
Даже 5 лет спустя, все еще очень полезно!
Многожильный Малыш
10

Вы можете определить экстент многоугольника, а затем ограничить генерацию случайных чисел для значений X и Y в пределах этих экстентов.

Основной процесс: 1) Определите maxx, maxy, minx, miny вершин многоугольника, 2) Создайте случайные точки, используя эти значения в качестве границ 3) Проверьте каждую точку на пересечение с вашим многоугольником, 4) Прекратите генерировать, когда у вас достаточно точек, удовлетворяющих пересечению тест

Вот алгоритм (C #) для теста пересечения:

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;
}
Дэн Уолтон
источник
10

Есть несколько хороших библиотек, которые делают большую часть тяжелой работы за вас.

Пример использования [shapely] [1] в 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)

Или используйте, .representative_point()чтобы получить точку внутри объекта (как упомянуто Дейном):

Возвращает дешевую вычисленную точку, которая гарантированно находится внутри геометрического объекта.

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

  [1]: https://shapely.readthedocs.io
monkut
источник
2
Разве это не должно быть из импорта shapely.geometry ...?
PyMapr
1
Можно также использовать representative_pointметод: shapely.readthedocs.io/en/latest/...
Дейн
6

Если R вариант, см. ?spsampleВ spпакете. Полигоны могут быть считаны из любого формата, поддерживаемого GDAL, встроенного в пакет rgdal, и затем spsampleработают непосредственно с импортированным объектом с различными вариантами выборки.

mdsumner
источник
+1 - Так как R является открытым исходным кодом, если кто-то хочет повторить, вы всегда можете зайти в источник, чтобы увидеть, как они сделаны. Для точечных шаблонов также могут быть заинтересованы инструменты моделирования в пакете spatstat.
Энди W
5

Я хотел бы предложить решение, которое требует очень мало с точки зрения ГИС-анализа. В частности, он не требует триангуляции каких-либо полигонов.

Следующий алгоритм, приведенный в псевдокоде, относится к некоторым простым операциям в дополнение к базовым возможностям обработки списка (создание, поиск длины, добавление, сортировка, извлечение подсписков и объединение) и генерации случайных чисел с плавающей точкой в ​​интервале [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.

Все они доступны практически в любой среде ГИС или графического программирования (и, если нет, легко кодируются). Clipне должен возвращать вырожденные полигоны (то есть те, которые имеют нулевую площадь).

Процедура SimpleRandomSampleэффективно получает список точек, случайно распределенных внутри многоугольника. Это обертка SRS, которая разбивает многоугольник на более мелкие кусочки до тех пор, пока каждый кусочек не станет достаточно компактным для эффективного отбора проб. Чтобы сделать это, он использует предварительно вычисленный список случайных чисел, чтобы решить, сколько очков выделить для каждой части.

SRS можно «настроить», изменив параметр t. Это максимальное ограничение: отношение площадей полигонов, которое может быть допустимым. Если сделать его маленьким (но больше 1), большинство полигонов будет разбито на множество частей; увеличение его размера может привести к отклонению многих пробных точек для некоторых полигонов (извилистых, с осколками или дырявыми). Это гарантирует, что максимальное время выборки исходного полигона предсказуемо.

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

Следующая процедура вызывает себя рекурсивно, если это необходимо. Таинственное выражение t*N + 5*Sqrt(t*N)консервативно оценивает верхний предел того, сколько точек потребуется, учитывая случайную изменчивость. Вероятность того, что это не удастся, составляет всего 0,3 на миллион вызовов процедур. Увеличьте 5 до 6 или даже 7, чтобы уменьшить эту вероятность, если хотите.

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
источник
2

Если ваш многоугольник является выпуклым, и вы знаете все вершины, вы можете рассмотреть возможность "случайного" выпуклого взвешивания вершин, чтобы выбрать новую точку, которая гарантированно лежит внутри выпуклой оболочки (в данном случае многоугольника).

Например, скажем, у вас есть N-сторонний выпуклый многоугольник с вершинами

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

Тогда генерируем случайным образом N выпуклых весов

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

Случайно выбранная точка затем определяется как

Y=  w_i*V_i

Может быть другой способ выбрать N выпуклых весов

  • Выбирайте числа N-1 равномерно случайным образом в пределах диапазона (без замены), сортируйте их и нормализуйте N интервалов между ними, чтобы получить веса.
  • Вы также можете получить выборку из распределения Дирихле, которое часто используется в качестве сопряженного априора для многочленного распределения, которое аналогично выпуклым весам в вашем случае.

Когда ваш многоугольник не очень сильно невыпуклый, вы можете сначала преобразовать его в выпуклый корпус. Это должно как минимум ограничить количество точек, лежащих за пределами вашего многоугольника в значительной степени.

algoseer
источник
2

Эту задачу очень легко решить в GRASS GIS (одна команда) с помощью v.random .

Ниже приведен пример того, как добавить 3 случайных точки в выбранные многоугольники (здесь - области с индексами города Роли, Северная Каролина) со страницы руководства. Изменяя оператор SQL «где», можно выбрать многоугольник (ы).

Генерация случайных точек в выбранных полигонах

markusN
источник
1
Обязательное напоминание о том, что почтовые индексы - это линии, а не полигоны.
Ричард
Можете ли вы уточнить? Для меня также здесь это относится к областям: en.wikipedia.org/wiki/ZIP_Code#Primary_state_prefixes
markusN
Конечно: почтовые индексы ссылаются на конкретные почтовые отделения и их маршруты доставки почты. В результате почтовые индексы являются линиями, а не многоугольниками. Они могут накладываться друг на друга, содержать дыры и не обязательно покрывать весь США или любой данный штат. Использование их для разделения области опасно по этой причине. Единицы переписи (например, группы блоков) - лучший выбор. Смотрите также: это и это .
Ричард
1
Благодарность! Вероятно, это также зависит от страны, см., Например, en.wikipedia.org/wiki/Postal_codes_in_Germany - однако почтовые индексы не являются моей основной темой, они просто хотели проиллюстрировать и ответить на оригинальный вопрос «Создание точек, которые лежат внутри многоугольника», а не обсудите определения почтового индекса, который является ОТ здесь :-)
markusN
1
Отмечено по обоим пунктам. Я, вероятно, должен сделать небольшую запись в блоге, чтобы я мог кратко сказать это в следующий раз :-)
Ричард
1

Ссылка для ответа

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

Три алгоритма, использующие разные подходы.

Git Repo Link

  1. Вот простой и лучший подход, использующий фактическое расстояние координат от направления x и y. Внутренний алгоритм использует WGS 1984 (4326) и преобразовывает результат во вставленный SRID.

Функция ======================================================= ==================

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;

Используйте функцию с простым запросом, геометрия должна быть правильной и полигон, многогранник или конверт

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

Результат ================================================= ===================== введите описание изображения здесь

  1. Вторая функция на основе Никласа Авена алгоритме . Пытался справиться с любым SRID.

    Я применил следующие изменения в алгоритме.

    1. Отдельная переменная для направления x и y для размера пикселя,
    2. Новая переменная для расчета расстояния в сфероиде или эллипсоиде.
    3. Введите любой SRID, функция преобразования Geom в рабочую среду Spheroid или Ellipsoid Datum, затем примените расстояние к каждой стороне, получите результат и преобразуйте его во входной SRID.

Функция ======================================================= ==================

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;

Используйте это с простым запросом.

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

Результат ================================================= ==================введите описание изображения здесь

  1. Функция основана на последовательном генераторе.

Функция ================================================= =================

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;

Используйте это с простым запросом.

SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons; Результат ================================================= =========================

введите описание изображения здесь

Мухаммед Имран Сиддик
источник