Generar puntos que se encuentran dentro del polígono.

30

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.

usuario2024
fuente
3
Por el contrario, el tiempo es predecible. Es proporcional a la proporción del área de extensión del polígono dividida por el área del polígono, multiplicada por el tiempo necesario para generar y probar un solo punto. El tiempo varía un poco, pero la variación es proporcional a la raíz cuadrada del número de puntos. Para números grandes, esto se vuelve intrascendente. Rompa los polígonos tortuosos en piezas más compactas si es necesario para reducir esta relación de área a un valor bajo. Solo los polígonos fractales te darán problemas, ¡pero dudo que los tengas!
whuber
3
posible duplicado de puntos
Pablo
@Pablo: buenos hallazgos. Sin embargo, ambas preguntas son específicas del software y ambas se refieren a colocar matrices regulares de puntos dentro de polígonos, no puntos aleatorios
whuber
de acuerdo con la diferencia de whuber es puntos aleatorios vs generación de puntos regulares dentro de un polígono
Mapperz

Respuestas:

20

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

Dan S.
fuente
2
+1 Simple y efectivo. Vale la pena señalar que se pueden generar puntos aleatorios uniformes dentro de un triángulo sin ningún rechazo, porque hay mapeos (fácilmente calculables) que preservan el área entre cualquier triángulo y un triángulo rectángulo isósceles, que es medio cuadrado, digamos la mitad donde la coordenada y excede la coordenada x. Genere dos coordenadas aleatorias y ordénelas para obtener un punto aleatorio en el triángulo isósceles, luego vuelva a asignarlo al triángulo original.
whuber
+1 Me gusta mucho la discusión de las coordenadas trilineales a las que hace referencia el artículo que cita. Supongo que esto sería susceptible a la esfera cuya superficie se representa como una teselación de triángulos. En un plano proyectado, no sería una distribución verdaderamente aleatoria, ¿verdad?
Kirk Kuykendall
@whuber - +1 de vuelta a ti. Otra forma (en el enlace, pero lo agitaron manualmente) es reflejar los puntos rechazados del cuadrilátero muestreado uniformemente a través del borde compartido y de vuelta al triángulo.
Dan S.
@Kirk: el enlace de citas es un poco útil en el sentido de que enumera un montón de métodos de muestreo incorrectos (no uniformes), incluidas las coordenadas trilineales, antes de la forma "correcta". No parece que haya una forma directa de obtener un muestreo uniforme con coordenadas trilineales. Me acercaría al muestreo uniforme en toda la esfera mediante la conversión de vectores unitarios aleatorios en 3d a su equivalente lat / lon, pero solo soy yo. (No estoy seguro acerca del muestreo restringido a triángulos / polígonos esféricos). (También no estoy seguro acerca del muestreo verdaderamente uniforme, por ejemplo, en wgs84: solo elegir ángulos sesgará un poco hacia los polos, creo).
Dan S.
1
@Dan Para muestrear uniformemente la esfera, use una proyección cilíndrica de igual área (las coordenadas son la longitud y el coseno de la latitud). Si desea muestrear sin usar una proyección, hay un truco hermoso: generar tres variables normales estándar independientes (x, y, z) y proyectarlas al punto (R x / n, R y / n, R * z / n ) donde n ^ 2 = x ^ 2 + y ^ 2 + z ^ 2 y R es el radio de la tierra. Convierte a (lat, lon) si es necesario (usando latitudes authalic cuando trabajas en un esferoide). Funciona porque esta distribución normal trivariada es esféricamente simétrica. Para muestrear triángulos, adhiérase a una proyección.
whuber
14

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

ingrese la descripción de la imagen aquí

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

bajo oscuro
fuente
1
Incluso 5 años después, ¡sigue siendo realmente útil!
Stranded Kid
10

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;
}
Dan Walton
fuente
10

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
Monkut
fuente
2
¿No debería ser de shapely.geometry import ...?
PyMapr
1
También puede usar el representative_pointmétodo: shapely.readthedocs.io/en/latest/…
desde el
6

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.

mdsumner
fuente
+1: dado que R es de código abierto, si uno quiere replicar, siempre puede ir a la fuente para ver cómo se hace. Para los patrones de puntos, uno también puede estar interesado en las herramientas de simulación en el paquete spatstat.
Andy W
5

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
}
whuber
fuente
2

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.

algoseer
fuente
2

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.

Generando puntos aleatorios en polígonos seleccionados

markusN
fuente
1
Recordatorio obligatorio de que los códigos postales son líneas, no polígonos.
Richard
¿Puedes elaborar? Para mí también aquí se refiere a áreas: en.wikipedia.org/wiki/ZIP_Code#Primary_state_prefixes
markusN
Claro: los códigos postales hacen referencia a oficinas postales particulares y sus rutas de entrega de correo. Como resultado, los códigos postales son líneas, no polígonos. Pueden superponerse entre sí, contener agujeros y no necesariamente cubren todo Estados Unidos o un estado determinado. Usarlos para dividir el área es peligroso por esta razón. Las unidades de censo (como los grupos de bloques) son una mejor opción. Ver también: esto y esto .
Richard
1
¡Gracias! Probablemente también depende del país, consulte por ejemplo en.wikipedia.org/wiki/Postal_codes_in_Germany ; sin embargo, los códigos postales no son mi tema central, solo quería ilustrar y responder la pregunta original "Generar puntos que se encuentran dentro del polígono" en lugar de discuta las definiciones del código postal que es OT aquí :-)
markusN
1
Destacado en ambos aspectos. Probablemente debería hacer una pequeña entrada en el blog para poder decirlo de manera más sucinta la próxima vez :-)
Richard
1

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 ================================================= ===================== ingrese la descripción de la imagen aquí

  1. 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 ================================================= ==================ingrese la descripción de la imagen aquí

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

ingrese la descripción de la imagen aquí

Muhammad Imran Siddique
fuente