Point in polygon - application sur des coordonnées GPS

Bonjour,
Alors voilà,
je dispose d'une fonction booléenne qui à partir d'un vecteurs de points V et d'un point P, détermine si le point est oui ou non à l'intérieur d'un polygone.
Les coordonnées ici sont dans un plan cartésien (X,Y).

Ma question :
Puis-je appliquer cette meme fonction mais en ayant des coordoonées GPS (longitude,Latitude) au lieu de coordonnées 2D (X,Y) ?
Mon but étant de délimiter une région sur une carte à partir d'un ensemble de points GPS (Polygone), puis d'appliquer cette fonction pour voir si un point est à l'intérieur de cette région ou non?

Merci pour toute réponse

Bien à vous.

Réda.

Bonjour, :D

Ce qui est délicat, c'est que ton outil se réfère à la géométrie du plan, dans lequel tout polygone est délimité par des segments de droite, tandis que les coordonnées géographiques font intervenir (au mieux) la géométrie sphérique; les contours "polygonaux" sont alors courbes (il peut s'agir d'arcs de grands cercles de la sphère), et leur projection sur un plan tangent ne conduira pas à des segments de droite.

Les réponses obtenues pourront apparaître contestables dans le cas de points "proches" de la frontière, et surtout lorsque la courbure de la surface n'est plus négligeable; le rapport du diamètre maximal (D_max) du polygone au rayon terrestre (R) n'étant alors plus petit devant l'unité: D_max/R >~ 0.1 .

De plus, les coordonnées angulaires poseront de grosse difficultés au niveau des pôles; il est heureusement peu probable que tu t'intéresses à ces régions.

Si tu te limites par contre à des régions d'étendue limitée (D_max <~ 100 km) et te contentes d'une précision raisonnable, tu pourras utiliser les formules de conversion entre systèmes de coordonnées (éventuellement sous une forme approchée)⁽¹⁾.

Cependant, cela me paraît assez difficile; ne serait-il pas plus simple d'envisager une réponse directe ? Cela suppose un minimum de familiarité dans le maniement des vecteurs.
Une autre solution pourrait être la recherche (éventuellement le téléchargement) de logiciels appropriés; d'autres intervenants pourraient te suggérer des liens sur ce sujet.

(1) Voir ce qui suit.

(1) Pour prolonger ce qui a été dit, avec les réserves déjà exprimées:

Au voisinage d'un point de longitude (a₀) et de latitude (b₀),
1) un écart de 1 degré en latitude correspond à la 360^me partie d'un méridien de rayon (R), donc à un arc de cercle de longueur
Dy = (2*Pi*R/360) ;
2) un écart de 1 degré en longitude correspond à la 360^me partie d'un parallèle de rayon (R*Cos(b₀)), donc à un arc de cercle de longueur Dx = (2*Pi*R*Cos(b₀)/360) .

Dans un plan tangent à la sphère au point M₀(a₀, b₀) et orienté par les vecteurs unitaires u_E (dirigé vers l'est) et u_N (dirigé vers le nord) la position d'un point M(a, b) est approximativement représentée par le vecteur
A₀M = Dx*(a - a₀).u_E + Dy*(b - b₀).u_N ,
ce qui revient à employer dans le repère local ainsi défini, les coordonnées cartésiennes:
X = Dx*(a - a₀) = (2*Pi*R*Cos(b₀))*((a - a₀)/360) et
Y = Dy*(b - b₀) = (2*Pi*R)*((b - b₀)/360) , les angles s'exprimant toujours en degrés.

Formule déconseillée pour le pistage des ours blancs d'Ellesmere Island, ou des manchots de Terre-Adélie :D .

Bonjour wiwaxia,
Avant toute chose, je vous remercie infiniment pour toutes ces explications bien détaillées et de votre effort !

Dans mon cas, le diamètre maximal du polygone est de l'ordre de quelques KMs seulement, je pense donc pouvoir me contenter des formules de conversion entre (longitude, latitude) vers (X,Y)


Si je comprends bien, le problème revient à calculer X et Y à partir d'un point de longitude a0 et latitude b0
soit déterminer les fonctions f1 et f2 tel que :
- X = f1(a0,b0)
- Y = f2(a0,b0).

En reprenant vos formules :

Citation:

---

X = Dx*(a - a0) = (2*Pi*R*Cos(b0))*((a - a0)/360) et
Y = Dy*(b - b0) = (2*Pi*R)*((b - b0)/360) , les angles s'exprimant toujours en degrés

---

je remarque que f1 et f2 font intervenir a0, b0, a et b
que représentent alors a et b ?
Merci encore

Réda

Citation:

---

Envoyé par b_reda31

... Dans mon cas, le diamètre maximal du polygone est de l'ordre de quelques KMs seulement, je pense donc pouvoir me contenter des formules de conversion entre (longitude, latitude) vers (X,Y) ...

---

Le rayon terrestre (R) admettant pour valeur moyenne 6371 km , les erreurs relatives affectant les distances et dues à la courbure de la surface n'atteignent pas un millionème; on a en effet:
e ~ (D_max/R)² ~ (4 / 6371)² = 4E^-7 .

Citation:

---

Envoyé par b_reda31

... Si je comprends bien, le problème revient à calculer X et Y à partir d'un point de longitude a0 et latitude b0
soit déterminer les fonctions f1 et f2 tel que :
- X = f1(a₀, b₀, a)
- Y = f2(b₀, b) ...

---

En effet, et à prendre le point (A) pour origine du repère local puisque l'on a par définition: X_A = Y_A = 0 .
On calculera de même les coordonnées des autres sommets du polygone: (X_B, Y_B), (X_C, Y_C), ... etc.

Citation:

---

Envoyé par b_reda31

En reprenant vos formules :
X = Dx*(a - a0) = (2*Pi*R*Cos(b0))*((a - a0)/360) et
Y = Dy*(b - b0) = (2*Pi*R)*((b - b0)/360) , les angles s'exprimant toujours en degrés
je remarque que f1 et f2 font intervenir a0, b0, a et b
que représentent alors a et b ?

---

La longitude et la latitude du point considéré (M) que l'on cherche à situer par rapport aux sommets (A, B, C, D ...) du polygone.
Les coordonnées angulaires (a_i, b_i)des autres points s'obtiennent par les mêmes relations; par exemple:
- X_B = X₁ = f1(a₀, b₀, a₁)
- Y_B = Y₁ = f2(b₀, b₁)
en affectant les sommets successifs du polygone (A, B, C ...) des valeurs (0, 1, 2 ... etc); l'indice (i) joue le rôle d'un compteur.

PS: J'ai corrigé une faute par omission, tardivement repérée: l'erreur relative est inférieure à la première estimation donnée.

Tu n'as pas besoin de faire de conversion. Les distances que tu traites sont petites, tout va bien...
Donc tu connais un algorithme qui permet de savoir si un point (X,Y) est à l'intérieur d'un polygône (Xi,Yi).
Tu peux parfaitement utiliser l'algorithme tel quel, en disant X=longitude, et Y=latitude.

Tant que les distances sont pas trop grandes (quelques kilomètres), et tant qu'on est assez loin des pôles (disons latitude entre -75° et +75°), ça ne pose aucun problème. Si tu te rapproches des pôles, il faut traiter des polygônes plus petits pour que ça reste valide. Et en particulier, si le polygône contient l'un des 2 pôles, l'approximation ne marche plus du tout.

Sur des latitudes 'classiques', en disant X=longitude, et Y=latitude, on fait une approximation. Mais c'est une approximation totalement négligeable. En gros, à 1 mètre près, les résultats sont corrects. Donc normalement, ça ne pose aucun problème, vu que l'erreur due au traitement est en fait inférieure à la précision du GPS lui-même.

Merci à vous deux.
j'ai appliqué l'algorithme directement sans passer par la conversion pour des petites surfaces, les résultats sont remarquables !
En voici une capture où les marqueurs rouges désignent les sommet du polygone , les marqueur vert, ceux détectés à l’intérieur du poly et les gris à l’extérieur
https://preview.ibb.co/czwSaS/Polygone_Carte.png

je vais tenter d'augmenter au fur et à mesure le diametre du polygone et voir à quel moment ca va clocher !

Merci encore à vous deux!

Le diamètre du polygone n'est pas forcément critique. Si tu as un polygone avec un diamètre 10 fois plus grand, mais avec 10 fois plus de points, la précision restera la même. Les problèmes peuvent apparaître si tu as des segments très longs, et d'autant plus si tu as des segments très longs qui sont orientés '''à 45%'''
Si tu as des segments très longs, qui sont parfaitement orientés nord/sud, aucun problème. Idem ou presque si tu as des segments très longs orientés est/ouest. Reste le cas où tu as des segments longs, orientés plus ou moins à 45°.
Mais pour un segment de 1 km, ce qui semble énorme par rapport aux données que tu traites, tu auras une imprécision de 1 mètre au maximum. Comme en plus tu traites des contours, qui se basent sur le centre de la route, et que tes marqueurs, sont des 'habitations', donc décalées de 3 ou 4 mètres au moins par rapport au centre de la route, tu es tranquille.
Et comme en plus, tu travailles sur une carte de l'Algérie, donc plus proche de l'équateur par rapport à la France, ça aussi, ça diminue les risques d'erreur.

Merci à vous deux. Je prends bien note de toutes vos remarques très instructives.

Bien à vous.
Réda

# L'injection directe des valeurs des coordonnées angulaires est effectivement le moyen le plus rapide de vérifier, à l'aide de la fonction disponible, la position d'un point par rapport à un polygone; la validité du procédé vient de ce l'on passe d'un type de coordonnées (a, b) à l'autre (x, y) par une transformation affine, qui ne modifie pas la disposition relative des points dans une figure.

Comme dans l'exemple ci-dessous, où les 2 groupes de points se déduisent l'un de l'autre par une affinité horizontale de rapport (2):
Pièce jointe 366142

# Le diamètre maximal (D_max) du polygone permet d'évaluer l'erreur relative affectant les distances dans le cas de la configuration la plus défavorable; il s'agit donc d'un seuil de sécurité, qui n'exclut pas des écarts beaucoup plus faibles, dans des cas très particuliers.
Le remplacement de ce paramètre global (D_max) par la plus longue des arêtes (A_max) supprimerait bien sûr tout désaccord, mais je doute que ce soit pertinent dans le cas d'un très grand nombre de sommets (comme celui de la figure proposée (#7, 30 sommets), où le second seuil est beaucoup plus petit (A_max < 0.1 * D_max); il faudrait connaître l'algorithme de la fonction conduisant au résultat.⁽¹⁾

# On a bien évoqué la nécessité de ne pas s'approcher des pôles, mais en restant sur ce point beaucoup trop évasifs. Une justification de cette contrainte peut être établie, en reconsidérant l'erreur relative résultant de l'assimilation d'un arc de cercle à une segment de droite:
a) pour un déplacement Nord-Sud suivant un cercle méridien, l'erreur correspondante est: e₁ ~ (d₁/R_mer)² = (d₁/R)²
et son bornage conduit à celui des distances: e₁ <= e_max entraînant: d₁ <~ R*e_max^1/2 = D_max1 ;
b) pour un déplacement Est-Ouest sur un cercle parallèle à l'équateur, l'erreur devient: e₂ ~ (d₂/R_par)² = (d₂/(R*Cos(b))²
et l'on obtient par les mêmes considérations: e₂ <= e_max , et d₂ <~ R*Cos(b)*e_max^1/2 = D_max2 = D_max1*Cos(b) .

Le diamètre maximal du polygone ne peut dépasser le plus petit de ces seuils, et doit donc vérifier:
D_max <~ Min(D_max1, D_max2) = D_max2 = R*Cos(b)*e_max^1/2 ;
la contrainte imposée est d'autant plus étroite que l'on se rapproche plus des pôles.
Ainsi pour une erreur tolérée (e_max) de l'ordre de un dix-millième (0.01 %), on obtient à l'équateur (b = 0): D_max <~ 64 km ,
et à une latitude de 75° : D_max <~ 16 km .

(1) Au fait, quel en est le principe ? Je crois deviner un procédé très simple, mais qui ne change en rien le problème présent.