gps/misc.py

   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