source: branches/numpy/anuga/coordinate_transforms/redfearn.py @ 6883

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

Based in part on spreadsheet
www.icsm.gov.au/gda/gdatm/redfearn.xls
downloaded from INTERGOVERNMENTAL COMMITTEE ON SURVEYING & MAPPING (ICSM)
http://www.icsm.gov.au/icsm/

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


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,
             zone=None, central_meridian=None, scale_factor=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

    If zone is specified reproject lat and long to specified zone instead of
    standard zone
    If meridian is specified, reproject lat and lon to that instead of zone. In this case
    zone will be set to -1 to indicate non-UTM projection

    Note that zone and meridian cannot both be specifed
    """


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


    #GDA Specifications
    a = 6378137.0                       #Semi major axis
    inverse_flattening = 298.257222101  #1/f
    if scale_factor is None:
        K0 = 0.9996                         #Central scale factor    
    else:
        K0 = scale_factor
    #print 'scale', K0
    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

    if zone is not None and central_meridian is not None:
        msg = 'You specified both zone and central_meridian. Provide only one of them'
        raise Exception, msg
    
    # Zone
    if zone is None:
        zone = int((lon - longitude_of_western_edge_zone0)/zone_width)

    # Central meridian
    if central_meridian is None:
        central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
    else:
        zone = -1

    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



def convert_from_latlon_to_utm(points=None,
                               latitudes=None,
                               longitudes=None,
                               false_easting=None,
                               false_northing=None):
    """Convert latitude and longitude data to UTM as a list of coordinates.


    Input

    points: list of points given in decimal degrees (latitude, longitude) or
    latitudes: list of latitudes   and
    longitudes: list of longitudes 
    false_easting (optional)
    false_northing (optional)

    Output

    points: List of converted points
    zone:   Common UTM zone for converted points


    Notes

    Assume the false_easting and false_northing are the same for each list.
    If points end up in different UTM zones, an ANUGAerror is thrown.    
    """

    old_geo = Geo_reference()    
    utm_points = []
    if points == None:
        assert len(latitudes) == len(longitudes)
        points =  map(None, latitudes, longitudes)
        
    for point in points:
        
        zone, easting, northing = redfearn(float(point[0]),
                                           float(point[1]),
                                           false_easting=false_easting,
                                           false_northing=false_northing)
        new_geo = Geo_reference(zone)
        old_geo.reconcile_zones(new_geo)        
        utm_points.append([easting, northing])

    return utm_points, old_geo.get_zone()