Bonjour,
On se place dans l'espace euclidien usuel.
On considère l'éllipsoide de paramètres et (demis-grands axes) définie par :
, pour et .
L'éllipsoide représente en fait un model approché de la Terre ( longitude ; latitude).
Et nous nous donnons 2 points de l'éllipsoide : et . Ces points sont exprimés en coordonnées éllipsoidiques dans un repère dont l'origine est le centre de notre l'éllipsoide ().
Après calculs trigonométriques, je trouve l'angle au centre entre et .
Mais comment calcul-t-on la distance (orthodromique si on en revient au modèle terrestre) exacte ? Dans le cas sphèrique , il suffirait d'effectuer ...
Bonsoir,
ton paramétrage me semble mal adapté.
La terre est à peu près un ellipsoïde de révolution.
L'axe de révolution étant l'axe des pôles.
Il me semblerai judicieux de prendre cet axe dans la base utilisée.
Traditionnellement c'est le troisième axe.
Et on a
Ce qui me semble plus simple.
Bonjour verdurin,
Tout d'abords merci de votre réponse !
Effectivement, j'ai inversé l'axe de révolution ... Une fois cette correction faite, le calcul de l'angle se simplifie un peu.
Mais ma question persiste toujours pour le calcul de la distance exacte de l'arc sur l'éllipsoide... En précisant que pour une éllipsoide approximant la Terre, les paramètres et étant relativement proche, une moyenne de ces 2 réels multipliée à donne le résultat avec une erreur correcte de l'ordre de . Mais dans le cas général, comment procédé avec ?
Si vous auriez une idée de la façon dont procédé, je suis preneur.
Effectivement, il est expliqué comment calculer la longueur d'un arc elliptique mais celui ci doit (si j'ai bien compris) être contenu dans un quadrant de l'ellipse...
Et en regardant les intégrales, je commence à réaliser la complexité des calculs ...
De plus, le calcul de l'équation cartésienne de l'ellipse intersection de l'ellipsoïde avec le plan vectoriel (centre de l'ellipsoïde considéré comme origine) passant par et devient vite compliquée dans le cas général..
Enfin bref, merci pour les indications !
Bonjour,
Merci sylvainc2, ces liens sont excellents !
Dans la section {Inverse problem} du lien Wikipedia, j'ai codé les formules puis tracé un graphique de l'erreur faite en fonction du nombre d'itérations. On y voit aussi la convergence de .
L'erreur tends vers 6400 km en gros ... Mais je ne comprends pas mon erreur ...
Ci-dessous le code (en python, mais si nécéssaire je peux l'expliquer en français parlé) pour calculer une distance orthodromique suivi du graphique évoqué :
def DISTANCE_ORTHODROMIQUE(Point1, Point2, Model = 'WGS84') :
# On importe les bon paramètres suivant le modèle choisi.
[a, b, f] = PARAMETRES_MODEL_ELLIPSOIDE_TERRE(Model)
# On convertit en mètres (comme sur Wikipedia).
a = 1000 * a
b = 1000 * b
# On importe les coordonnées Longitude / Latitude en radian de 2 points.
[Long1, Lat1] = CONVERT_COORDONNEES_DEGRES_TO_RADIANS(Point1)
[Long2, Lat2] = CONVERT_COORDONNEES_DEGRES_TO_RADIANS(Point2)
# Chaque variable calculée porte le même nom que sur Wikipedia pour rester clair.
U1 = math.atan((1 - f) * math.tan(Lat1))
U2 = math.atan((1 - f) * math.tan(Lat2))
L = Long1 - Long2
# Nombre d'itérations qu'on effectuera : 1000
for iter in range(0, 1000) :
sin_sigma = math.sqrt(pow(math.cos(U2)*math.sin(L), 2)+pow(math.cos(U1)*math.sin(U2) \
-math.sin(U1)*math.cos(U2)*math.cos(L), 2))
cos_sigma = math.sin(U1)*math.sin(U2) + math.cos(U1)*math.cos(U2)*math.cos(L)
sigma = math.atan2(sin_sigma, cos_sigma)
sin_alpha = (math.cos(U1)*math.cos(U2)*math.sin(L))/sin_sigma
cos_2sigma_m = cos_sigma - (2*math.sin(U1)*math.sin(U2))/(1-pow(sin_alpha, 2))
C = (f/16)*(1-pow(sin_alpha, 2))*(4 + f*(4 - 3*(1-pow(sin_alpha, 2))))
lamb = L + (1-C)*f*sin_alpha*(sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma*(-1+2*pow(cos_2sigma_m, 2))))
L = lamb
u_carre = (1-pow(sin_alpha, 2))*((a*a-b*b)/(b*b))
A = 1 + (u_carre/16384)*(4096 + u_carre*(-768 + u_carre*(320 - 175*u_carre)))
B = (u_carre/1024)*(256 + u_carre*(-128 + u_carre*(74 - 47*u_carre)))
Delta_sigma = B * sin_sigma * (cos_2sigma_m + 0.25*B * (cos_sigma * (-1 + 2*pow(cos_2sigma_m, 2)) \
- B/6 * cos_2sigma_m * (-3 + 4*pow(sin_sigma, 2)) * (-3 + 4*pow(cos_2sigma_m, 2))))
s = b*A*(sigma - Delta_sigma)
return s
La variable "L" ne doit pas changer à l'intérieur de la boucle, c'est la variable que tu nommes "lamb" qui change. Regarde bien dans la page wiki, il y a une distinction entre L et .
Aussi tu ne devrais pas coder le nombre d'itérations avec un "for", tu devrais faire comme ceci:
while True:
...
avant la ligne lamb = L + ...etc... tu mets:
D = lamb
puis pour fermer la boucle:
if (not(math.fabs(lamb-D)>EPSILON)):
break
où EPSILON est une constante que tu définis du genre EPSILON = 0.0000000001
Ca devrait fonctionner maintenant.
Et puis pour être clair à l'intérieur de la boucle c'est la variable soit "lamb" qui doit être utilisée dans les sin et cos, pas L sauf à la ligne
lamb = L + (1-C)*f ... etc...
Là c'est bien L, qui ne doit pas changer je le répète.
Donc tu dois initialiser lamb = L avant d'entrer dans la boucle pour que ca fonctionne correctement.
Bonjour,
Merci pour les corrections. J'ai modifié le code en conséquence :
def DISTANCE_ORTHODROMIQUE_3(Point1, Point2, Model = 'WGS84') :
# On importe les bon paramètres suivant le modèle choisi.
[a, b, f] = PARAMETRES_MODEL_ELLIPSOIDE_TERRE(Model)
# On convertit en mètres (comme sur Wikipedia).
a = 1000 * a
b = 1000 * b
# On importe les coordonnées Longitude / Latitude en radian de 2 points.
[Long1, Lat1] = CONVERT_COORDONNEES_DEGRES_TO_RADIANS(Point1)
[Long2, Lat2] = CONVERT_COORDONNEES_DEGRES_TO_RADIANS(Point2)
# Chaque variable calculée porte le même nom que sur Wikipedia pour rester clair.
U1 = math.atan((1 - f) * math.tan(Lat1))
U2 = math.atan((1 - f) * math.tan(Lat2))
L = Long1 - Long2
tolerance = 0.000000000001
lamb = L
while True :
D = lamb
sin_sigma = math.sqrt(pow(math.cos(U2)*math.sin(lamb), 2)+pow(math.cos(U1)* \
math.sin(U2)-math.sin(U1)*math.cos(U2)*math.cos(lamb), 2))
cos_sigma = math.sin(U1)*math.sin(U2) + math.cos(U1)*math.cos(U2)*math.cos(lamb)
sigma = math.atan2(sin_sigma, cos_sigma)
sin_alpha = (math.cos(U1)*math.cos(U2)*math.sin(lamb))/sin_sigma
cos_2sigma_m = cos_sigma - (2*math.sin(U1)*math.sin(U2))/(1-pow(sin_alpha, 2))
C = (f/16)*(1-pow(sin_alpha, 2))*(4 + f*(4 - 3*(1-pow(sin_alpha, 2))))
lamb = L + (1-C)*f*sin_alpha*(sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma* \
(-1+2*pow(cos_2sigma_m, 2))))
if abs(lamb - D) < tolerance :
break
u_carre = (1-pow(sin_alpha, 2))*((a*a-b*b)/(b*b))
A = 1 + (u_carre/16384)*(4096 + u_carre*(-768 + u_carre*(320 - 175*u_carre)))
B = (u_carre/1024)*(256 + u_carre*(-128 + u_carre*(74 - 47*u_carre)))
Delta_sigma = B * sin_sigma * (cos_2sigma_m + 0.25*B * (cos_sigma * (-1 \
+ 2*pow(cos_2sigma_m, 2)) - B/6 * cos_2sigma_m * (-3 + 4*pow(sin_sigma, 2)) \
* (-3 + 4*pow(cos_2sigma_m, 2))))
s = b*A*(sigma - Delta_sigma)
return s/1000
Je comprends pas, quand tu parles d'erreur, c'est par rapport à quoi? Parce que j'ai ressorti mon vieux code en langage C et je l'ai comparé avec le tien en python avec quelques jeux de test et les deux me donnent le même résultat au mm près pour la distance. D'après moi, ton programme est correct.
Aussi il existe des packages en python déjà faits, par exemple geopy et une autre référence donnée en bas de page ici: https://movable-type.co.uk/scripts/latlong-vincenty.html
Bonjour,
Par exemple on se donne les coordonnées :
Point1 = [138.363531, 84.946308]
Point2 = [-22.095115, 70.148730]
>>> print(DISTANCE_ORTHODROMIQUE_3(Point1, Point2), "km - Modèle WGS84")
2754.54660537193 km - Modèle WGS-84
Je viens de reprendre la fonction .distance.vincenty du module geopy évoqué :
def measure(Point1, Point2):
[lng1, lat1] = CONVERT_COORDONNEES_DEGRES_TO_RADIANS(Point1)
[lng2, lat2] = CONVERT_COORDONNEES_DEGRES_TO_RADIANS(Point2)
[major, minor, f] = PARAMETRES_MODEL_ELLIPSOIDE_TERRE('WGS-84')
delta_lng = lng2 - lng1
reduced_lat1 = math.atan((1 - f) * math.tan(lat1))
reduced_lat2 = math.atan((1 - f) * math.tan(lat2))
sin_reduced1, cos_reduced1 = math.sin(reduced_lat1), math.cos(reduced_lat1)
sin_reduced2, cos_reduced2 = math.sin(reduced_lat2), math.cos(reduced_lat2)
lambda_lng = delta_lng
lambda_prime = 2 * pi
iter_limit = 100 # Je crois que dans l'original c'était 20 ...
i = 0
while (i == 0 or (abs(lambda_lng - lambda_prime) > 10e-12 and i <= iter_limit)):
i += 1
sin_lambda_lng, cos_lambda_lng = math.sin(lambda_lng), math.cos(lambda_lng)
sin_sigma = math.sqrt( (cos_reduced2 * sin_lambda_lng) ** 2 + \
(cos_reduced1 * sin_reduced2 - sin_reduced1 * cos_reduced2 * cos_lambda_lng) ** 2 )
if sin_sigma == 0:
return 0 # Coincident points
cos_sigma = ( sin_reduced1 * sin_reduced2 + cos_reduced1 * cos_reduced2 * cos_lambda_lng )
sigma = math.atan2(sin_sigma, cos_sigma)
sin_alpha = ( cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma )
cos_sq_alpha = 1 - sin_alpha ** 2
if cos_sq_alpha != 0:
cos2_sigma_m = cos_sigma - 2 * ( sin_reduced1 * sin_reduced2 / cos_sq_alpha )
else:
cos2_sigma_m = 0.0 # Equatorial line
C = f / 16. * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))
lambda_prime = lambda_lng
lambda_lng = ( delta_lng + (1 - C) * f * sin_alpha * ( sigma + C * sin_sigma * ( cos2_sigma_m + C \
* cos_sigma * ( -1 + 2 * cos2_sigma_m ** 2 ) ) ) )
if i > iter_limit:
raise ValueError("Vincenty formula failed to converge!")
u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2
A = 1 + u_sq / 16384. * ( 4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)) )
B = u_sq / 1024. * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
delta_sigma = ( B * sin_sigma * ( cos2_sigma_m + B / 4. * ( cos_sigma * ( -1 + 2 * cos2_sigma_m ** 2 \
) - B / 6. * cos2_sigma_m * ( -3 + 4 * sin_sigma ** 2 ) * ( -3 + 4 * cos2_sigma_m ** 2 ) ) ) )
s = minor * A * (sigma - delta_sigma)
return s
>>> print(GPS.DISTANCE_ORTHODROMIQUE_3(A[3], B[3]), "km.")
print(GPS.measure(A[6], B[6]), "km - Fonction reprise de geopy.")
2754.54660537193 km.
2754.5466053719297 km - Fonction reprise de geopy.
Google maps semble calculer la distance sur grand cercle, avec un rayon terrestre de 6371km qui est la valeur moyenne acceptée du rayon. Avec ton exemple ci-haut j'obtiens aussi 2743122.395855m sur grand cercle. Donc il ne faut pas comparer la formule de Vincenty avec google maps, c'est pas le même calcul.
Remarque qu'il y a une version de Vincenty qui tient compte de l'altitude des points au-dessus ou sous l'ellipsoide de référence, je pensais que c'était ca que google maps calculait, mais non finalement.
Effectivement tu as raison. Je suis tout même surpris que Google Map se contente d'une telle approximation ... Enfin, on peut leur pardonner puisque la version Globe vient récemment de remplacer la version Mercator ...
Enfin bref, merci pour tout le temps que tu as pu consacrer à cette conversation.
Vous devez être membre accéder à ce service...
Pas encore inscrit ?
1 compte par personne, multi-compte interdit !
Ou identifiez-vous :