source: branches/numpy/anuga/coordinate_transforms/lat_long_UTM_conversion.py @ 6442

#!/usr/bin/env python

# Lat Long - UTM, UTM - Lat Long conversions
#
# see http://www.pygps.org
#

from math import pi, sin, cos, tan, sqrt

#LatLong- UTM conversion..h
#definitions for lat/long to UTM and UTM to lat/lng conversions
#include <string.h>

_deg2rad = pi / 180.0
_rad2deg = 180.0 / pi

_EquatorialRadius = 2
_eccentricitySquared = 3

_ellipsoid = [
#  id, Ellipsoid name, Equatorial Radius, square of eccentricity
# first once is a placeholder only, To allow array indices to match id numbers
        [ -1, "Placeholder", 0, 0],
        [ 1, "Airy", 6377563, 0.00667054],
        [ 2, "Australian National", 6378160, 0.006694542],
        [ 3, "Bessel 1841", 6377397, 0.006674372],
        [ 4, "Bessel 1841 (Nambia] ", 6377484, 0.006674372],
        [ 5, "Clarke 1866", 6378206, 0.006768658],
        [ 6, "Clarke 1880", 6378249, 0.006803511],
        [ 7, "Everest", 6377276, 0.006637847],
        [ 8, "Fischer 1960 (Mercury] ", 6378166, 0.006693422],
        [ 9, "Fischer 1968", 6378150, 0.006693422],
        [ 10, "GRS 1967", 6378160, 0.006694605],
        [ 11, "GRS 1980", 6378137, 0.00669438],
        [ 12, "Helmert 1906", 6378200, 0.006693422],
        [ 13, "Hough", 6378270, 0.00672267],
        [ 14, "International", 6378388, 0.00672267],
        [ 15, "Krassovsky", 6378245, 0.006693422],
        [ 16, "Modified Airy", 6377340, 0.00667054],
        [ 17, "Modified Everest", 6377304, 0.006637847],
        [ 18, "Modified Fischer 1960", 6378155, 0.006693422],
        [ 19, "South American 1969", 6378160, 0.006694542],
        [ 20, "WGS 60", 6378165, 0.006693422],
        [ 21, "WGS 66", 6378145, 0.006694542],
        [ 22, "WGS-72", 6378135, 0.006694318],
        [ 23, "WGS-84", 6378137, 0.00669438]
]

#Reference ellipsoids derived from Peter H. Dana's website-
#http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html
#Department of Geography, University of Texas at Austin
#Internet: pdana@mail.utexas.edu
#3/22/95

#Source
#Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
#1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency

#def LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long,
#                        double &UTMNorthing, double &UTMEasting, char* UTMZone)

def LLtoUTM( Lat, Long, ReferenceEllipsoid=23):
    """
    converts lat/long to UTM coords.  Equations from USGS Bulletin 1532
    East Longitudes are positive, West longitudes are negative.
    North latitudes are positive, South latitudes are negative
    Lat and Long are in decimal degrees
    Written by Chuck Gantz- chuck.gantz@globalstar.com
    """
    a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
    eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
    k0 = 0.9996

#Make sure the longitude is between -180.00 .. 179.9
    LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9

    LatRad = Lat*_deg2rad
    LongRad = LongTemp*_deg2rad

    ZoneNumber = int((LongTemp + 180)/6) + 1

    if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0:
        ZoneNumber = 32

    # Special zones for Svalbard
    if Lat >= 72.0 and Lat < 84.0:
        if  LongTemp >= 0.0  and LongTemp <  9.0:ZoneNumber = 31
        elif LongTemp >= 9.0  and LongTemp < 21.0: ZoneNumber = 33
        elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35
        elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37

    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 #+3 puts origin in middle of zone
    LongOriginRad = LongOrigin * _deg2rad

    #compute the UTM Zone from the latitude and longitude
    UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat))

    eccPrimeSquared = (eccSquared)/(1-eccSquared)
    N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad))
    T = tan(LatRad)*tan(LatRad)
    C = eccPrimeSquared*cos(LatRad)*cos(LatRad)
    A = cos(LatRad)*(LongRad-LongOriginRad)

    M = a*((1
            - eccSquared/4
            - 3*eccSquared*eccSquared/64
            - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
           - (3*eccSquared/8
              + 3*eccSquared*eccSquared/32
              + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
           + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
           - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad))

    UTMEasting = (k0*N*(A+(1-T+C)*A*A*A/6
                        + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
                  + 500000.0)

    UTMNorthing = (k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
                                        + (61
                                           -58*T
                                           +T*T
                                           +600*C
                                           -330*eccPrimeSquared)*A*A*A*A*A*A/720)))

    if Lat < 0:
        UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere
    #UTMZone was originally returned here.  I don't know what the
    #letter at the end was for.
    #print "UTMZone", UTMZone 
    return (ZoneNumber, UTMEasting, UTMNorthing)


def _UTMLetterDesignator(Lat):
#This routine determines the correct UTM letter designator for the given latitude
#returns 'Z' if latitude is outside the UTM limits of 84N to 80S
#Written by Chuck Gantz- chuck.gantz@globalstar.com

    if 84 >= Lat >= 72: return 'X'
    elif 72 > Lat >= 64: return 'W'
    elif 64 > Lat >= 56: return 'V'
    elif 56 > Lat >= 48: return 'U'
    elif 48 > Lat >= 40: return 'T'
    elif 40 > Lat >= 32: return 'S'
    elif 32 > Lat >= 24: return 'R'
    elif 24 > Lat >= 16: return 'Q'
    elif 16 > Lat >= 8: return 'P'
    elif  8 > Lat >= 0: return 'N'
    elif  0 > Lat >= -8: return 'M'
    elif -8> Lat >= -16: return 'L'
    elif -16 > Lat >= -24: return 'K'
    elif -24 > Lat >= -32: return 'J'
    elif -32 > Lat >= -40: return 'H'
    elif -40 > Lat >= -48: return 'G'
    elif -48 > Lat >= -56: return 'F'
    elif -56 > Lat >= -64: return 'E'
    elif -64 > Lat >= -72: return 'D'
    elif -72 > Lat >= -80: return 'C'
    else: return 'Z'    # if the Latitude is outside the UTM limits

#void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
#                         double& Lat,  double& Long )

def UTMtoLL(northing, easting, zone, isSouthernHemisphere=True,
            ReferenceEllipsoid=23):
    """
    converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
    East Longitudes are positive, West longitudes are negative.
    North latitudes are positive, South latitudes are negative
    Lat and Long are in decimal degrees.
    Written by Chuck Gantz- chuck.gantz@globalstar.com
    Converted to Python by Russ Nelson <nelson@crynwr.com>

    FIXME: This is set up to work for the Southern Hemisphere.

Using
http://www.ga.gov.au/geodesy/datums/redfearn_geo_to_grid.jsp

    Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid) 
Zone:   36    
Easting:  511669.521  Northing: 19328195.112 
Latitude:   84  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 '' 
Grid Convergence:  0  -59 ' 40.28 ''  Point Scale: 0.99960166

____________
Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid) 
Zone:   36    
Easting:  519384.803  Northing: 1118247.585 
Latitude:   -80  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 '' 
Grid Convergence:  0  59 ' 5.32 ''  Point Scale: 0.99960459 
___________
Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid) 
Zone:   36    
Easting:  611263.812  Northing: 10110547.106 
Latitude:   1  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 '' 
Grid Convergence:  0  -1 ' 2.84 ''  Point Scale: 0.99975325 
______________
Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid) 
Zone:   36    
Easting:  611263.812  Northing: 9889452.894 
Latitude:   -1  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 '' 
Grid Convergence:  0  1 ' 2.84 ''  Point Scale: 0.99975325 

So this uses a false northing of 10000000 in the both hemispheres.
ArcGIS used a false northing of 0 in the northern hem though.
Therefore it is difficult to actually know what hemisphere you are in.
    """
    k0 = 0.9996
    a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
    eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
    e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared))

    x = easting - 500000.0 #remove 500,000 meter offset for longitude
    y = northing

    ZoneNumber = int(zone)
    if isSouthernHemisphere:
        y -= 10000000.0         # remove 10,000,000 meter offset used
                                # for southern hemisphere

    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3  # +3 puts origin in middle of zone

    eccPrimeSquared = (eccSquared)/(1-eccSquared)

    M = y / k0
    mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256))

    phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
               + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
               +(151*e1*e1*e1/96)*sin(6*mu))
    phi1 = phi1Rad*_rad2deg;

    N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad))
    T1 = tan(phi1Rad)*tan(phi1Rad)
    C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad)
    R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5)
    D = x/(N1*k0)

    Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
                                          +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720)
    Lat = Lat * _rad2deg

    Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
            *D*D*D*D*D/120)/cos(phi1Rad)
    Long = LongOrigin + Long * _rad2deg
    return (Lat, Long)

if __name__ == '__main__':
    (z, e, n) = LLtoUTM(45.00, -75.00, 23)
    print z, e, n
    (lat, lon) = UTMtoLL(n, e, z, 23)
    print lat, lon