Package killerbee :: Package zbwardrive :: Package gps :: Module misc
[hide private]
[frames] | no frames]

Source Code for Module killerbee.zbwardrive.gps.misc

  1  # misc.py - miscellaneous geodesy and time functions 
  2  # 
  3  # This file is Copyright (c) 2010 by the GPSD project 
  4  # BSD terms apply: see the file COPYING in the distribution root for details. 
  5   
  6  import time, calendar, math 
  7   
  8  # some multipliers for interpreting GPS output 
  9  METERS_TO_FEET  = 3.2808399     # Meters to U.S./British feet 
 10  METERS_TO_MILES = 0.00062137119 # Meters to miles 
 11  KNOTS_TO_MPH    = 1.1507794     # Knots to miles per hour 
 12  KNOTS_TO_KPH    = 1.852         # Knots to kilometers per hour 
 13  KNOTS_TO_MPS    = 0.51444444    # Knots to meters per second 
 14  MPS_TO_KPH      = 3.6           # Meters per second to klicks/hr 
 15  MPS_TO_MPH      = 2.2369363     # Meters/second to miles per hour 
 16  MPS_TO_KNOTS    = 1.9438445     # Meters per second to knots 
 17   
 18  # EarthDistance code swiped from Kismet and corrected 
 19   
20 -def Deg2Rad(x):
21 "Degrees to radians." 22 return x * (math.pi/180)
23
24 -def Rad2Deg(x):
25 "Radians to degress." 26 return x * (180/math.pi)
27
28 -def CalcRad(lat):
29 "Radius of curvature in meters at specified latitude." 30 a = 6378.137 31 e2 = 0.081082 * 0.081082 32 # the radius of curvature of an ellipsoidal Earth in the plane of a 33 # meridian of latitude is given by 34 # 35 # R' = a * (1 - e^2) / (1 - e^2 * (sin(lat))^2)^(3/2) 36 # 37 # where a is the equatorial radius, 38 # b is the polar radius, and 39 # e is the eccentricity of the ellipsoid = sqrt(1 - b^2/a^2) 40 # 41 # a = 6378 km (3963 mi) Equatorial radius (surface to center distance) 42 # b = 6356.752 km (3950 mi) Polar radius (surface to center distance) 43 # e = 0.081082 Eccentricity 44 sc = math.sin(Deg2Rad(lat)) 45 x = a * (1.0 - e2) 46 z = 1.0 - e2 * sc * sc 47 y = pow(z, 1.5) 48 r = x / y 49 50 r = r * 1000.0 # Convert to meters 51 return r
52
53 -def EarthDistance((lat1, lon1), (lat2, lon2)):
54 "Distance in meters between two points specified in degrees." 55 x1 = CalcRad(lat1) * math.cos(Deg2Rad(lon1)) * math.sin(Deg2Rad(90-lat1)) 56 x2 = CalcRad(lat2) * math.cos(Deg2Rad(lon2)) * math.sin(Deg2Rad(90-lat2)) 57 y1 = CalcRad(lat1) * math.sin(Deg2Rad(lon1)) * math.sin(Deg2Rad(90-lat1)) 58 y2 = CalcRad(lat2) * math.sin(Deg2Rad(lon2)) * math.sin(Deg2Rad(90-lat2)) 59 z1 = CalcRad(lat1) * math.cos(Deg2Rad(90-lat1)) 60 z2 = CalcRad(lat2) * math.cos(Deg2Rad(90-lat2)) 61 a = (x1*x2 + y1*y2 + z1*z2)/pow(CalcRad((lat1+lat2)/2), 2) 62 # a should be in [1, -1] but can sometimes fall outside it by 63 # a very small amount due to rounding errors in the preceding 64 # calculations (this is prone to happen when the argument points 65 # are very close together). Thus we constrain it here. 66 if abs(a) > 1: a = 1 67 elif a < -1: a = -1 68 return CalcRad((lat1+lat2) / 2) * math.acos(a)
69
70 -def MeterOffset((lat1, lon1), (lat2, lon2)):
71 "Return offset in meters of second arg from first." 72 dx = EarthDistance((lat1, lon1), (lat1, lon2)) 73 dy = EarthDistance((lat1, lon1), (lat2, lon1)) 74 if lat1 < lat2: dy *= -1 75 if lon1 < lon2: dx *= -1 76 return (dx, dy)
77
78 -def isotime(s):
79 "Convert timestamps in ISO8661 format to and from Unix time." 80 if type(s) == type(1): 81 return time.strftime("%Y-%m-%dT%H:%M:%S", time.gmtime(s)) 82 elif type(s) == type(1.0): 83 date = int(s) 84 msec = s - date 85 date = time.strftime("%Y-%m-%dT%H:%M:%S", time.gmtime(s)) 86 return date + "." + `msec`[2:] 87 elif type(s) == type(""): 88 if s[-1] == "Z": 89 s = s[:-1] 90 if "." in s: 91 (date, msec) = s.split(".") 92 else: 93 date = s 94 msec = "0" 95 # Note: no leap-second correction! 96 return calendar.timegm(time.strptime(date, "%Y-%m-%dT%H:%M:%S")) + float("0." + msec) 97 else: 98 raise TypeError
99 100 # End 101