Générer des points situés à l'intérieur du polygone

J'ai une entité polygonale et je veux pouvoir générer des points à l'intérieur. J'en ai besoin pour une tâche de classification.

Générer des points aléatoires jusqu'à ce que l'on soit à l'intérieur du polygone ne fonctionnerait pas, car le temps que cela prend est vraiment imprévisible.

Commencez par décomposer le polygone en triangles, puis générez des points à l'intérieur de ceux-ci . (Pour une distribution uniforme, pondérez chaque triangle par sa surface.)

Lorsque vous mettez une balise QGIS sur cette question: l'outil Points aléatoires peut être utilisé avec une couche limite.

Si vous recherchez du code, le code source du plugin sous-jacent devrait vous être utile.

Vous pouvez déterminer l'étendue du polygone, puis contraindre la génération de nombres aléatoires pour les valeurs X et Y dans ces étendues.

Processus de base: 1) Déterminer maxx, maxy, minx, miny des sommets du polygone, 2) Générer des points aléatoires en utilisant ces valeurs comme limites 3) Tester chaque point pour l'intersection avec votre polygone, 4) Arrêter de générer lorsque vous avez suffisamment de points satisfaisant l'intersection tester

Voici un algorithme (C #) pour le test d'intersection:

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

Il existe de bonnes bibliothèques qui font le gros du travail pour vous.

Exemple utilisant [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)

Ou utilisez .representative_point()pour obtenir un point dans l'objet (comme mentionné par Dain):

Renvoie un point calculé à moindre coût qui se trouve à l'intérieur de l'objet géométrique.

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

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

Si R est une option, voir ?spsampledans le sppackage. Les polygones peuvent être lus à partir de n'importe quel format pris en charge par GDAL intégré au package rgdal, puis spsamplefonctionnent directement sur l'objet importé avec une variété d'options d'échantillonnage.

Je voudrais proposer une solution qui nécessite très peu en termes d'analyse SIG. En particulier, il ne nécessite la triangulation d'aucun polygone.

L'algorithme suivant, donné en pseudocode, fait référence à quelques opérations simples en plus des capacités de gestion de liste de base (créer, trouver la longueur, ajouter, trier, extraire des sous-listes et concaténer) et la génération de flottants aléatoires dans l'intervalle [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.

Ils sont tous disponibles dans presque tous les environnements de programmation SIG ou graphiques (et faciles à coder sinon). Clipne doit pas renvoyer de polygones dégénérés (c'est-à-dire ceux avec une zone nulle).

La procédure SimpleRandomSampleobtient efficacement une liste de points répartis aléatoirement dans un polygone. C'est un emballage pour SRS, qui divise le polygone en petits morceaux jusqu'à ce que chaque morceau soit suffisamment compact pour être échantillonné efficacement. Pour ce faire, il utilise une liste précalculée de nombres aléatoires pour décider du nombre de points à allouer à chaque pièce.

Le SRS peut être "réglé" en changeant le paramètre t. Il s'agit du rapport maximum du cadre de délimitation: surface surfacique qui peut être toléré. Le rendre petit (mais supérieur à 1) provoquera la division de la plupart des polygones en plusieurs morceaux; l'agrandir peut entraîner le rejet de nombreux points d'essai pour certains polygones (sinueux, avec des éclats ou plein de trous). Cela garantit que la durée maximale d'échantillonnage du polygone d'origine est prévisible.

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

La procédure suivante s'appelle récursivement si nécessaire. L'expression mystérieuse t*N + 5*Sqrt(t*N)estime prudemment une limite supérieure du nombre de points nécessaires, ce qui tient compte de la variabilité aléatoire. La probabilité que cela échoue n'est que de 0,3 par million d'appels de procédure. Augmentez de 5 à 6 ou même 7 pour réduire cette probabilité si vous le souhaitez.

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 votre polygone est convexe et que vous connaissez tous les sommets, vous voudrez peut-être envisager de faire une pondération convexe "aléatoire" des sommets pour échantillonner un nouveau point qui est garanti de se trouver à l'intérieur de la coque convexe (polygone dans ce cas).

Par exemple, disons que vous avez un polygone convexe à N côtés avec des sommets

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

Générer ensuite aléatoirement N poids convexes

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

Le point échantillonné au hasard est alors donné par

Y= ∑ w_i*V_i

Il peut y avoir différentes manières d'échantillonner N poids convexes

• Choisissez N-1 nombres uniformément au hasard dans une plage (sans remplacement), triez-les et normalisez les N intervalles entre eux pour obtenir les poids.
• Vous pouvez également échantillonner à partir de la distribution de Dirichlet qui est souvent utilisée comme un conjugué avant pour la distribution multinomiale, qui est similaire aux poids convexes dans votre cas.

Lorsque votre polygone n'est pas très sévèrement non convexe, vous pouvez d'abord envisager de le convertir en coque convexe. Cela devrait au moins limiter le nombre de points se trouvant en dehors de votre polygone dans une large mesure.

La tâche est très facile à résoudre dans GRASS GIS (une commande) en utilisant v.random .

Ci-dessous un exemple sur la façon d'ajouter 3 points aléatoires dans les polygones sélectionnés (ici les zones de code postal de la ville de Raleigh, NC) à partir de la page de manuel. En modifiant l'instruction SQL "where", le ou les polygones peuvent être sélectionnés.

Lien de réponse

/gis//a/307204/103524

Trois algorithmes utilisant différentes approches.

Lien Git Repo

1. Voici une approche simple et optimale, utilisant la distance réelle des coordonnées des directions x et y. L'algorithme interne utilise le WGS 1984 (4326) et la transformation des résultats en SRID inséré.

Fonction ================================================= ==================

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;

Utilisez la fonction avec une requête simple, la géométrie doit être valide et polygone, multi-polygones ou enveloppe

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

Résultat ================================================= =====================

2. Deuxième fonction basée sur l' algorithme de Nicklas Avén . J'ai essayé de gérer n'importe quel SRID.

   J'ai appliqué les changements suivants dans l'algorithme.

   1. Variable séparée pour les directions x et y pour la taille des pixels,
   2. Nouvelle variable pour calculer la distance en sphéroïde ou ellipsoïde.
   3. Saisissez n'importe quel SRID, transformez la fonction Geom dans l'environnement de travail de Spheroid ou Ellipsoid Datum, puis appliquez la distance de chaque côté, obtenez le résultat et transformez-le en SRID d'entrée.

Fonction ================================================= ==================

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;

Utilisez-le avec une simple requête.

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

Résultat ================================================= ==================

3. Fonction basée sur un générateur en série.

Fonction ================================================= =================

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;

Utilisez-le avec une simple requête.

SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons; Résultat ================================================= =========================