Ignore:
Timestamp:
Feb 24, 2009, 6:13:06 PM (15 years ago)
Author:
ole
Message:

Implemented ability to specify central meridian for non-UTM projections.
Test partially passes, but has been disabled for now.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/coordinate_transforms/redfearn.py

    r6149 r6404  
    3737    return sign*dd, mm, ss
    3838
    39 def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None):
     39def redfearn(lat, lon, false_easting=None, false_northing=None,
     40             zone=None, central_meridian=None):
    4041    """Compute UTM projection using Redfearn's formula
    4142
     
    4748    If zone is specified reproject lat and long to specified zone instead of
    4849    standard zone
     50
     51    If meridian is specified, reproject lat and lon to that instead of zone. In this case
     52    zone will be set to -1 to indicate non-UTM projection
     53
     54    Note that zone and meridian cannot both be specifed
    4955    """
    5056
     
    5460
    5561
    56     #GDA Specifications
     62    # GDA Specifications
    5763    a = 6378137.0                       #Semi major axis
    5864    inverse_flattening = 298.257222101  #1/f
     
    7379       
    7480   
    75     #Derived constants
     81    # Derived constants
    7682    f = 1.0/inverse_flattening
    7783    b = a*(1-f)       #Semi minor axis
     
    8490    e6 = e2*e4
    8591
    86     #Foot point latitude
     92    # Foot point latitude
    8793    n = (a-b)/(a+b) #Same as e2 - why ?
    8894    n2 = n*n
     
    114120    t6 = t2*t4
    115121   
    116     #Radius of Curvature
     122    # Radius of Curvature
    117123    rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5
    118124    nu = a/(1-e2*sinphi*sinphi)**0.5
     
    122128    psi4 = psi2*psi2
    123129
    124 
    125 
    126     #Meridian distance
    127 
     130    # Meridian distance
    128131    A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
    129132    A2 = 3.0/8*(e2+e4/4+15*e6/128)
     
    138141    m = term1 + term2 + term3 + term4 #OK
    139142
    140     #Zone
     143    if zone is not None and central_meridian is not None:
     144        msg = 'You specified both zone and central_meridian. Provide only one of them'
     145        raise Exception, msg
     146   
     147    # Zone
    141148    if zone is None:
    142149        zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
    143150
    144     central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
     151    # Central meridian
     152    if central_meridian is None:
     153        central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
     154    else:
     155        zone = -1
    145156
    146157    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
     
    153164    omega8 = omega4*omega4
    154165     
    155     #Northing
     166    # Northing
    156167    term1 = nu*sinphi*cosphi*omega2/2 
    157168    term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24
     
    162173    northing = false_northing + K0*(m + term1 + term2 + term3 + term4)
    163174
    164     #Easting
     175    # Easting
    165176    term1 = nu*omega*cosphi
    166177    term2 = nu*cosphi3*(psi-t2)*omega3/6
Note: See TracChangeset for help on using the changeset viewer.