Vive la géométrie sphérique !
Le dimanche 20 mars 2011 à 19:40 - Lien permanent
Soit le problème simple de calculer la longueur d'une balade, donnée par exemple sous forme d'une trace GPX, c'est-à-dire une suite de points dont on connaît les coordonnées dans le système WGS 84.
Ma première idée a été de projeter les points sur le plan (en utilisant par exemple une projection transverse de Mercator dans le bon fuseau UTM), puis de calculer la distance entre deux points successifs en utilisant bêtement le théorème de Pythagore. Cela fonctionne très bien, à condition que la projection soit valide pour toute la trace GPX. Aucun problème dans le cas d'une balade, mais pour un trajet réalisé en avion, qui s'étend sur plusieurs fuseaux UTM, cela n'est plus le cas. De plus, les formules de projection sont souvent assez compliquées.
En réalité, il est beaucoup plus simple, et plus universel, de raisonner directement en géométrie sphérique et de calculer les longueurs des arcs de grands cercles entre deux points successifs (les grands cercles sont les géodésiques d'une sphère, c'est-à-dire que les arcs de grands cercles sont les plus courts chemins entre deux points). Ceci est indépendant de toute projection.
Il suffit de calculer l'angle de l'arc de grand cercle entre deux points :
$$\Delta\theta = 2\arcsin\sqrt{\sin^2{\left(\frac{\varphi' - \varphi}{2}\right)} + \cos{\varphi} \cdot \cos{\varphi'} \cdot \sin^2{\left(\frac{\lambda' - \lambda}{2}\right)}\ }$$\( \varphi\) et \( \varphi'\) sont les latitudes des deux points ; \( \lambda\) et \( \lambda'\) leurs longitudes.
La longueur est alors donnée par \( R \Delta\theta \). J'ai pris la même valeur du rayon que celle choisie dans l'API Google Maps[1], la valeur du demi-grand axe de l'ellipsoïde du système WGS 84, soit \(R = 6378,137 \mathrm{~km}\).
En bonus, la fonction Python ci-dessous effectue le calcul, à partir de coordonnées exprimées en radians :
import math def distance(fromLat, fromLon, toLat, toLon): radius = 6378.137 deltaLat = toLat - fromLat deltaLon = toLon - fromLon angle = 2 * math.asin( math.sqrt( \ math.pow(math.sin(deltaLat/2), 2) + \ math.cos(fromLat) * math.cos(toLat) * \ math.pow(math.sin(deltaLon/2), 2) ) ) return radius * angle
Au final, c'est beaucoup plus simple que d'effectuer une projection, et cela fonctionne dans tous les cas !
Notes
[1] Le calcul est déjà implémenté en Javascript dans le package de géométrie sphérique de l'API Google Maps