source: anuga_core/source/anuga/coordinate_transforms/redfearn.py @ 3613

"""Implementation of Redfearn's formula to compute UTM projections from latitude and longitude

"""
from anuga.coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE
from Numeric import array

def degminsec2decimal_degrees(dd,mm,ss):
    assert abs(mm) == mm
    assert abs(ss) == ss    

    if dd < 0:
        sign = -1
    else:
        sign = 1
    
    return sign * (abs(dd) + mm/60. + ss/3600.)      
    
def decimal_degrees2degminsec(dec):

    if dec < 0:
        sign = -1
    else:
        sign = 1

    dec = abs(dec)    
    dd = int(dec)
    f = dec-dd

    mm = int(f*60)
    ss = (f*60-mm)*60
    
    return sign*dd, mm, ss

def redfearn(lat, lon, false_easting=None, false_northing=None):
    """Compute UTM projection using Redfearn's formula

    lat, lon is latitude and longitude in decimal degrees

    If false easting and northing are specified they will override
    the standard
    """


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


    #GDA Specifications
    a = 6378137.0                       #Semi major axis
    inverse_flattening = 298.257222101  #1/f
    K0 = 0.9996                         #Central scale factor    
    zone_width = 6                      #Degrees

    longitude_of_central_meridian_zone0 = -183    
    longitude_of_western_edge_zone0 = -186

    if false_easting is None:
        false_easting = 500000

    if false_northing is None:
        if lat < 0:
            false_northing = 10000000  #Southern hemisphere
        else:
            false_northing = 0         #Northern hemisphere)
        
    
    #Derived constants
    f = 1.0/inverse_flattening
    b = a*(1-f)       #Semi minor axis

    e2 = 2*f - f*f#    = f*(2-f) = (a^2-b^2/a^2   #Eccentricity
    e = sqrt(e2)
    e2_ = e2/(1-e2)   # = (a^2-b^2)/b^2 #Second eccentricity
    e_ = sqrt(e2_)
    e4 = e2*e2
    e6 = e2*e4

    #Foot point latitude
    n = (a-b)/(a+b) #Same as e2 - why ?
    n2 = n*n
    n3 = n*n2
    n4 = n2*n2

    G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180


    phi = lat*pi/180     #Convert latitude to radians

    sinphi = sin(phi)   
    sin2phi = sin(2*phi)
    sin4phi = sin(4*phi)
    sin6phi = sin(6*phi)

    cosphi = cos(phi)
    cosphi2 = cosphi*cosphi
    cosphi3 = cosphi*cosphi2
    cosphi4 = cosphi2*cosphi2
    cosphi5 = cosphi*cosphi4    
    cosphi6 = cosphi2*cosphi4
    cosphi7 = cosphi*cosphi6
    cosphi8 = cosphi4*cosphi4        

    t = tan(phi)
    t2 = t*t
    t4 = t2*t2
    t6 = t2*t4
    
    #Radius of Curvature
    rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5
    nu = a/(1-e2*sinphi*sinphi)**0.5
    psi = nu/rho
    psi2 = psi*psi
    psi3 = psi*psi2
    psi4 = psi2*psi2



    #Meridian distance

    A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
    A2 = 3.0/8*(e2+e4/4+15*e6/128)
    A4 = 15.0/256*(e4+3*e6/4)
    A6 = 35*e6/3072
    
    term1 = a*A0*phi
    term2 = -a*A2*sin2phi
    term3 = a*A4*sin4phi
    term4 = -a*A6*sin6phi

    m = term1 + term2 + term3 + term4 #OK

    #Zone
    zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
    central_meridian = zone*zone_width+longitude_of_central_meridian_zone0

    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
    omega2 = omega*omega
    omega3 = omega*omega2
    omega4 = omega2*omega2
    omega5 = omega*omega4
    omega6 = omega3*omega3
    omega7 = omega*omega6
    omega8 = omega4*omega4
    
    
    #Northing
    term1 = nu*sinphi*cosphi*omega2/2  
    term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24
    term3 = nu*sinphi*cosphi5*\
            (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\
             psi2*(1-32*t2)-psi*2*t2+t4-t2)*omega6/720
    term4 = nu*sinphi*cosphi7*(1385-3111*t2+543*t4-t6)*omega8/40320
    northing = false_northing + K0*(m + term1 + term2 + term3 + term4)

    #Easting
    term1 = nu*omega*cosphi
    term2 = nu*cosphi3*(psi-t2)*omega3/6
    term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120
    term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040
    easting = false_easting + K0*(term1 + term2 + term3 + term4)
    
    return zone, easting, northing


# FIXME (Ole) 
def convert_points_from_latlon_to_utm(latitudes, longitudes,
                                      false_easting=None,
                                      false_northing=None):
    """Wrapper
    """
    
    return convert_lats_longs(latitudes,
                              longitudes,
                              false_easting,
                              false_northing)

    

def convert_lats_longs(latitudes,
                       longitudes,
                       false_easting=None,
                       false_northing=None):
    """
    Redfearn a list of latitudes and longitudes.

    Assume the false_easting and false_northing are the same for each list.

    Note, if a zone turns out to be different, an ANUGAerror is thrown.
    """
    # Using geo_ref so there is only one point in the code where zones
    # are checked.
    old_geo = Geo_reference()
    points = []
    for lat, long in map(None, latitudes, longitudes,):
        zone, easting, northing = redfearn(float(lat),
                                           float(long),
                                           false_easting=false_easting,
                                           false_northing=false_northing)
        new_geo = Geo_reference(zone)
        old_geo.reconcile_zones(new_geo)
        #print "old_geo.get_zone()", old_geo.get_zone()
        points.append([easting, northing])
        
    return old_geo.get_zone(), points