1 #! /usr/bin/env python 2 # -*- coding: Latin-1 -*- 3 4 __author__ = "Jean-Philippe Durand (jpdurand@free.fr)" 5 __version__ = "$Revision: 1.0 $" 6 __date__ = "$Date: jeu mai 4 14:01:02 CEST 2006 $" 7 __copyright__ = "Copyright (c) 2006 J.-Ph. Durand" 8 __license__ = "Python" 9 10 import os 11 import sys 12 #import math 13 from math import cos 14 from math import sin 15 from math import acos 16 from math import asin 17 from math import atan2 18 from math import degrees 19 from math import radians 20 from math import sqrt 21 22 rt = 6371 # rayon terrestre moyen en km 23 24 class Point: 25 def __init__(self,lat,lon): 26 self.lat=radians(lat) 27 self.lon=radians(lon) 28 def afficher(self): 29 print "lat=",degrees(self.lat) 30 print "long=",degrees(self.lon) 31 32 # 33 # Calcul de la distance (en km) entre 2 points specifies par leurs 34 # latitude/longitude avec la formule de Haversine 35 # 36 # de : Haversine formula - R. W. Sinnott, "Virtues of the Haversine", 37 # Sky and Telescope, vol 68, no 2, 1984 38 # http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1 39 # 40 41 42 43 def distHaversine(p1, p2): 44 dLat = p2.lat - p1.lat 45 dLong = p2.lon - p1.lon 46 a = sin(dLat/2) * sin(dLat/2) + cos(p1.lat) * cos(p2.lat) * sin(dLong/2) * sin(dLong/2) 47 c = 2 * atan2(sqrt(a), sqrt(1-a)) 48 d = rt * c 49 return d 50 51 # 52 # Calcul de la distance (en km) entre 2 points specifies par leurs 53 # latitude/longitude en utilisant les fonctions trigonometriques 54 # 55 def distCosineLaw(p1, p2): 56 d = acos(sin(p1.lat)*sin(p2.lat) + cos(p1.lat)*cos(p2.lat)*cos(p2.lon-p1.lon)) * rt 57 return d 58 59 60 # 61 # calcul du cap (en degres ) entre 2 points 62 # 63 # de : Ed Williams' Aviation Formulary, http://williams.best.vwh.net/avform.htm#Crs 64 # 65 def cap(p1, p2): 66 y = sin(p2.lon-p1.lon) * cos(p2.lat) 67 x = cos(p1.lat)*sin(p2.lat) - sin(p1.lat)*cos(p2.lat)*cos(p2.lon-p1.lon) 68 return degrees(atan2(y, x)) 69 70 # 71 # calcul le point de destination a parit du point de depart, du cap et de la distance 72 # voir http://williams.best.vwh.net/avform.htm#LL 73 # 74 def destPoint(p1, capinit, dist): 75 if (dist !=0): 76 p2=Point(0,0) 77 d = float(dist)/rt # d = distance angulaire parcourue sur la surface terrestre 78 capinit = radians(capinit) 79 p2.lat = asin( sin(p1.lat)*cos(d) + cos(p1.lat)*sin(d)*cos(capinit) ) 80 p2.lon = p1.lon + atan2(sin(capinit)*sin(d)*cos(p1.lat), cos(d)-sin(p1.lat)*sin(p2.lat)) 81 else: 82 p2=p1 83 return p2 84 85 86 87 p1=Point(30,60) 88 p2=Point(38,61) 89 p1.afficher() 90 p2.afficher() 91 dist = distHaversine(p1, p2) 92 print dist 93 dist = distCosineLaw(p1, p2) 94 print dist 95 capinit=cap(p1,p2) 96 print capinit 97 p3=destPoint(p1, 35.2643896828, 6671.69559867) 98 p3.afficher() 99 print degrees(p3.lat) 100 print degrees(p3.lon)