Generar puntos que se encuentran dentro del polígono.

Tengo una función de polígono y quiero poder generar puntos dentro de ella. Necesito esto para una tarea de clasificación.

Generar puntos aleatorios hasta que uno esté dentro del polígono no funcionaría porque es realmente impredecible el tiempo que lleva.

Comience descomponiendo el polígono en triángulos, luego genere puntos dentro de ellos . (Para una distribución uniforme, pondera cada triángulo por su área).

Al poner una etiqueta QGIS en esta pregunta: la herramienta Puntos aleatorios se puede usar con una capa límite.

Si está buscando código, el código fuente del complemento subyacente debería ser de ayuda.

Puede determinar la extensión del polígono, luego restringir la generación de números aleatorios para los valores X e Y dentro de esas extensiones.

Proceso básico: 1) Determine los vértices maxx, maxy, minx, miny de polígono, 2) Genere puntos aleatorios usando estos valores como límites 3) Pruebe cada punto para la intersección con su polígono, 4) Deje de generar cuando tenga suficientes puntos que satisfagan la intersección prueba

Aquí hay un algoritmo (C #) para la prueba de intersección:

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

Hay algunas buenas bibliotecas que hacen la mayor parte del trabajo pesado por usted.

Ejemplo usando [shapely] [1] en 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)

O use .representative_point()para obtener un punto dentro del objeto (como lo menciona dain):

Devuelve un punto calculado de forma económica que se garantiza que está dentro del objeto geométrico.

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

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

Si R es una opción, ver ?spsampleen el sppaquete. Los polígonos se pueden leer desde cualquier formato compatible con GDAL integrado en el paquete rgdal, y luego spsamplefunciona directamente en objetos importados con una variedad de opciones de muestreo.

Me gustaría ofrecer una solución que requiera muy poco en términos de análisis SIG. En particular, no requiere triangular ningún polígono.

El siguiente algoritmo, dado en pseudocódigo, se refiere a algunas operaciones simples además de las capacidades básicas de manejo de listas (crear, encontrar longitud, agregar, ordenar, extraer sublistas y concatenar) y la generación de flotantes aleatorios en el intervalo [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.

Todos estos están disponibles en casi cualquier entorno de programación gráfica o SIG (y fácil de codificar si no). Clipno debe devolver polígonos degenerados (es decir, aquellos con área cero).

Procedimiento SimpleRandomSampleobtiene eficientemente una lista de puntos distribuidos aleatoriamente dentro de un polígono. Es un envoltorio para SRS, que rompe el polígono en piezas más pequeñas hasta que cada pieza sea lo suficientemente compacta como para tomar muestras de manera eficiente. Para hacer esto, utiliza una lista precalculada de números aleatorios para decidir cuántos puntos asignar a cada pieza.

SRS se puede "ajustar" cambiando el parámetro t. Este es el cuadro de límite máximo: relación de área de polígono que se puede tolerar. Hacerlo pequeño (pero mayor que 1) hará que la mayoría de los polígonos se dividan en muchos pedazos; hacerlo grande puede hacer que se rechacen muchos puntos de prueba para algunos polígonos (sinuosos, con astillas o llenos de agujeros). Esto garantiza que el tiempo máximo para muestrear el polígono original es predecible.

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

El siguiente procedimiento se llama a sí mismo recursivamente si es necesario. La expresión misteriosa t*N + 5*Sqrt(t*N)estima conservadoramente un límite superior sobre cuántos puntos se necesitarán, lo que explica la variabilidad del azar. La probabilidad de que esto falle es de solo 0.3 por millón de llamadas a procedimientos. Aumente de 5 a 6 o incluso 7 para reducir esta probabilidad si lo desea.

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
}

Si su polígono es convexo y conoce todos los vértices, es posible que desee considerar hacer una ponderación convexa "aleatoria" de los vértices para muestrear un nuevo punto que se garantiza que se encuentra dentro del casco convexo (polígono en este caso).

Por ejemplo, supongamos que tiene un polígono convexo de N lados con vértices

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

Luego genera aleatoriamente N pesos convexos

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

El punto muestreado al azar viene dado por

Y= ∑ w_i*V_i

Puede haber diferentes formas de muestrear N pesos convexos

• Elija números N-1 de manera uniforme y aleatoria dentro de un rango (sin reemplazo), ordénelos y normalice los N intervalos entre ellos para obtener los pesos.
• También puede tomar muestras de la distribución de Dirichlet, que a menudo se usa como conjugado antes de la distribución multinomial, que es similar a los pesos convexos en su caso.

Cuando su polígono no es muy no convexo, puede considerar primero convertirlo en un casco convexo. Esto debería al menos limitar en gran medida el número de puntos que se encuentran fuera de su polígono.

La tarea es muy fácil de resolver en GRASS GIS (un comando) usando v.random .

A continuación se muestra un ejemplo sobre cómo agregar 3 puntos aleatorios en polígonos seleccionados (aquí áreas de código postal de la ciudad de Raleigh, NC) desde la página del manual. Al modificar la instrucción SQL "donde" se pueden seleccionar los polígonos.

Enlace de respuesta

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

Tres algoritmos que utilizan diferentes enfoques.

Git Repo Link

1. Aquí hay un enfoque simple y mejor, utilizando la distancia real de las coordenadas desde la dirección x e y. El algoritmo interno usa el WGS 1984 (4326) y el resultado se transforma en SRID insertado.

Función ================================================= ==================

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;

Utilice la función con una consulta simple, la geometría debe ser válida y polígono, polígonos múltiples o envolvente

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

Resultado ================================================= =====================

2. Segunda función basada en el algoritmo de Nicklas Avén . He intentado manejar cualquier SRID.

   He aplicado los siguientes cambios en el algoritmo.

   1. Variable separada para la dirección x e y para el tamaño de píxel,
   2. Nueva variable para calcular la distancia en esferoide o elipsoide.
   3. Ingrese cualquier SRID, la función transforma Geom al entorno de trabajo de Spheroid o Ellipsoid Datum, luego aplique la distancia a cada lado, obtenga el resultado y transforme para ingresar SRID.

Función ================================================= ==================

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;

Úselo con una consulta simple.

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

Resultado ================================================= ==================

3. Función basada en generador en serie.

Función ================================================= =================

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;

Úselo con una consulta simple.

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