datum.py

gps pygarmin-0[1].6.tgz

"""

   Datums and other conversions.
   This is meant to be pretty and understandable rather than fast.

   Geodetic datums are models of the shape of the earth. The earth is
   rather irregular and attempts to simplify it by assuming that it is
   basically an ellipsoid, for example, tend to be reasonably accurate
   only for a local area. This is why so many datums are in use.

   The most common datum for international purposes is WGS84, and the
   parameters of the other ones are usually quoted as differences from
   WGS84.  In this module. a Datum object represents a datum and
   converts lat/lon points in that datum to and from WGS84.
   
   A reference datum consists of a particular shape of ellipsoid and
   an x, y, z offset.  The ellipsoid is defined by its semi-major-axis
   'a' (in this case, the equatorial radius) and its semi-minor-axis
   'b' (the polar radius).  In most cases, rather than using 'b'
   directly, we are interested in the difference between the two - the
   degree to which the earth is considered to be 'flattened' at the
   poles. So the flattening, f, is given by (a-b)/a.

   An example use:

     import datum
     osd = datum.DatumFromName ('Ordnance Survey of Great Britain 36')
     print osd.toWGS84deg(12.3,54.3) 

   (c) Quentin Stafford-Fraser <[e-mail]>
       AT&T Laboratories Cambridge 2001
   Partly derived from some python code by Joseph Newman
   which in turn came from C code by Alan Jones

"""

import math
import refdatum

# These are for convenience:

def DEG_TO_RAD(angle): return angle*0.0174532925199

def RAD_TO_DEG(angle): return angle/0.0174532925199

# First we define an ellipsoid

class Ellipsoid:
    def __init__(self, a, invf):
        """
        Create an ellipsoid.  Flattenings tend to be incoveniently
        small numbers (typically 0.0033), so they are often quoted
        as their reciprocal.  The invf parameter expected here is 1/f.
        The value es is the first eccentricity squared, which is
        useful later.
        """
        self.a = a
        self.f = 1.0/invf
        self.es = 2 * self.f  -  self.f * self.f


def EllipsoidFromName(name):
    """
    Create a standard ellipsoid from the parameters in
    refdatum.py
    """
    re =  refdatum.Ellipsoids[name]
    return Ellipsoid(re[0], re[1])



class Datum:
    def __init__(self, ellipsoid, dx, dy, dz):
        """
          All parameters give the WGS84 values relative to this datum.
          So dx is the WGS84 X value minus the local datum X value, etc.
          dx, dy, dz are the X,Y,Z offsets.
        """
        global WGS84Datum
        self.dx = dx
        self.dy = dy
        self.dz = dz
        # self.ell = ellipsoid
        self.a = ellipsoid.a
        self.f = ellipsoid.f
        self.wgse = EllipsoidFromName("WGS 84")
        self.da = self.wgse.a - self.a
        self.df = self.wgse.f - self.f

    def fromWGS84rad(self, lat, lon, alt=0.0):
        return self.molodensky(lat, lon, alt,
                               -self.dx, -self.dy, -self.dz,
                               self.wgse.a, self.wgse.f,
                               -self.da, -self.df)

    def fromWGS84deg(self, lat, lon, alt=0.0):
        latrad = DEG_TO_RAD(lat)
        lonrad = DEG_TO_RAD(lon)
        newlat, newlon, newalt = self.fromWGS84rad(latrad, lonrad, alt)
        return (RAD_TO_DEG(newlat), RAD_TO_DEG(newlon), newalt)

    def toWGS84rad(self, lat, lon, alt=0.0):
        return self.molodensky(lat, lon, alt,
                               self.dx, self.dy, self.dz,
                               self.a, self.f,
                               self.da, self.df)

    def toWGS84deg(self, lat, lon, alt=0.0):
        latrad = DEG_TO_RAD(lat)
        lonrad = DEG_TO_RAD(lon)
        newlat, newlon, newalt = self.toWGS84rad(latrad, lonrad, alt)
        return (RAD_TO_DEG(newlat), RAD_TO_DEG(newlon), newalt)

    # everything in radians here
    def molodensky(self, lat, lon,  alt, dx, dy, dz, from_a, from_f, da, df):
        sinlat = math.sin (lat)
        coslat = math.cos (lat)
        sinlon = math.sin (lon)
        coslon = math.cos (lon)
        sin1 = math.sin (math.pi / (3600.0 * 180.0))
        bdiva = 1.0 - from_f
        esq = (2.0 - from_f)*from_f
        exp = 1.0 - esq*sinlat*sinlat
        sqrtexp = math.sqrt (exp)
        rn = from_a / sqrtexp
        rm = from_a * (1.0 - esq) / (exp*sqrtexp)
        dlat = (- dx*sinlat*coslon - dy*sinlat*sinlon + dz*coslat +
                da*(rn*esq*sinlat*coslat)/from_a +
                df*(rm/bdiva + rn*bdiva)*sinlat*coslat) / (rm + alt)
        dlon = (-dx*sinlon +dy*coslon)/((rn+alt)*coslat)
        dh = (dx*coslat*coslon +
              dy*coslat*sinlon +
              dz*sinlat-
              da*from_a/rn +
              df*bdiva*rn*sinlat*sinlat)
        newlat = lat + dlat
        newlon = lon + dlon
        newalt = alt + dh
        return (newlat, newlon, newalt)


# Then we use the ellipsoids to define datums
def DatumFromName(name):
    rd = refdatum.Datums[name]
    e = EllipsoidFromName(rd[0])
    return Datum(e, rd[1], rd[2], rd[3])


# A simple test routine

def test():
     osd = DatumFromName ('Ordnance Survey of Great Britain 36')
     nad = DatumFromName ('North America 1927 mean')
     wlat, wlon, walt    = (12.3, 45.6, 78)
     print "lat, lon, alt =", wlat, wlon, walt
     olat, olon, oalt    = osd.fromWGS84deg(wlat,wlon,walt)
     wlat2, wlon2, walt2 = osd.toWGS84deg(olat, olon, oalt)
     print "Under OSGB datum =", olat, olon, oalt
     print "errors after conversion back = ",
     print wlat2-wlat, wlon2-wlon, walt2-walt

if __name__ == '__main__':
    test()