Calcul de distance entre coordonnées (lat/lon) sur terre

[PhM]

Bonjour,

Aujourd'hui, juste pour s'amuser, je me suis penché sur les méthodes de calculs de distance sur terre à partir des coordonnées lat/lon de deux lieux différents.

Il y a essentiellement 3 méthodes disponibles depuis longtemps (http://villemin.gerard.free.fr/aGeograp/Distance.htm) :

• Méthode simple de Pythagore

• Méthode avec loi des sinus

• Méthode avec la formule de Haversine (la plus précise)

Je me suis amusé à programmer suivant la formule de Haversine et voici le résultat :

Code : Tout sélectionner

; -------------------------------------------------------------------------------------------------------
;                       Distance entre deux coordonnées géographiques sur terre
;                                 suivant la formule de Haversine
;                       Philippe Mijon - décembre 2023 PureBasic 6.03 LTS (x64)
;
; -------------------------------------------------------------------------------------------------------

Declare haversine(lat1.f, lon1.f, lat2.f, lon2.f)

Global PI.f = 3.141592653589793       ; Nombre PI avec 15 décimales
Global ra.f = 6371.0087714            ; Rayon moyen de la terre en km
Global DIST.f

      ; Coordonnées de 2 lieux en valeurs décimales ex : KMSP,Minneapolis,MN,USA/KMKE,Milwaukee,WI,USA
      
      lat1.f = 44.8780191
      lon1.f = -93.2209281
      
      lat2.f = 42.946680
      lon2.f = -87.896750
      
      haversine (lat1.f, lon1.f, lat2.f, lon2.f)
     
      Debug StrD(DIST.f,2) + " Km"
      Debug StrD(DIST.f/1.852,2) + " NM"
      
      ; Contrôle des calculs sur https://fr.distance.To
      RunProgram ("https://fr.distance.To/" + lat1.f + "," + lon1.f + "/" + lat2.f + "," + lon2.f)

End

Procedure haversine(lat1.f, lon1.f, lat2.f, lon2.f)

        ; différences entre latitudes et longitudes
        dLat.f = (lat2.f - lat1.f) * (PI / 180)
        dLon.f = (lon2.f - lon1.f) * (PI / 180)
        
        ; conversion en radians
        lat1.f = (lat1.f) * (PI / 180)
        lat2.f = (lat2.f) * (PI / 180)
        
        ; formule de Haversine
        a.f = Pow(Sin(dLat.f / 2), 2) + Pow(Sin(dLon.f / 2), 2) * Cos(lat1.f) * Cos(lat2.f)
        c.f = 2 * ASin(Sqr(a.f))
        DIST.f = ra.f * c.f

EndProcedure

[PhM]

Variante de la conversion en radians :
; conversion en radians
lat1.f = Radian(lat1.f)
lat2.f = Radian(lat2.f)

[PhM]

Avec la conversion en radians dans la formule :

Code : Tout sélectionner

Procedure haversine(lat1.f, lon1.f, lat2.f, lon2.f)

        ; différences entre latitudes et longitudes
        dLat.f = (lat2.f - lat1.f) * (PI / 180)
        dLon.f = (lon2.f - lon1.f) * (PI / 180)
        
        ; formule de Haversine
        a.f = Pow(Sin(dLat.f / 2), 2) + Pow(Sin(dLon.f / 2), 2) * Cos(Radian(lat1.f)) * Cos(Radian(lat2.f))
        c.f = 2 * ASin(Sqr(a.f))
        DIST.f = ra.f * c.f

EndProcedure

[PhM]

La version, moins lisible, mais très condensée du même programme :

Code : Tout sélectionner

; -------------------------------------------------------------------------------------------------------
;                       Distance entre deux coordonnées géographiques sur terre
;                                 suivant la formule de Haversine
;                       Philippe Mijon - décembre 2023 PureBasic 6.03 LTS (x64)
;      RunProgram ("https://fr.distance.To/" + lat1.f + "," + lon1.f + "/" + lat2.f + "," + lon2.f)
;
; -------------------------------------------------------------------------------------------------------

; Coordonnées de 2 lieux en valeurs décimales ex : KMSP,Minneapolis,MN,USA/KMKE,Milwaukee,WI,USA

lat1.f = 44.8780191
lon1.f = -93.2209281

lat2.f = 42.946680
lon2.f = -87.896750

; formule de Haversine
DIST.f = 6371.0087714 * (2 * ASin(Sqr(Pow(Sin(Radian(lat2.f - lat1.f) / 2), 2) + Pow(Sin(Radian(lon2.f - lon1.f) / 2), 2) * Cos(Radian(lat1.f)) * Cos(Radian(lat2.f)))))

Debug StrD(DIST.f,2) + " Km (" + StrD(DIST.f/1.852,2) + " NM)"

End