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

"""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
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, zone=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
    """


    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
    if zone is None:
        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



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()