# Changeset 6404 for anuga_core/source/anuga/coordinate_transforms/redfearn.py

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

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

File:
1 edited

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

 r6149 return sign*dd, mm, ss def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None): def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None, central_meridian=None): """Compute UTM projection using Redfearn's formula 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 """ #GDA Specifications # GDA Specifications a = 6378137.0                       #Semi major axis inverse_flattening = 298.257222101  #1/f #Derived constants # Derived constants f = 1.0/inverse_flattening b = a*(1-f)       #Semi minor axis e6 = e2*e4 #Foot point latitude # Foot point latitude n = (a-b)/(a+b) #Same as e2 - why ? n2 = n*n t6 = t2*t4 #Radius of Curvature # Radius of Curvature rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5 nu = a/(1-e2*sinphi*sinphi)**0.5 psi4 = psi2*psi2 #Meridian distance # Meridian distance A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256 A2 = 3.0/8*(e2+e4/4+15*e6/128) m = term1 + term2 + term3 + term4 #OK #Zone 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 = zone*zone_width+longitude_of_central_meridian_zone0 # 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) omega8 = omega4*omega4 #Northing # Northing term1 = nu*sinphi*cosphi*omega2/2 term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24 northing = false_northing + K0*(m + term1 + term2 + term3 + term4) #Easting # Easting term1 = nu*omega*cosphi term2 = nu*cosphi3*(psi-t2)*omega3/6
Note: See TracChangeset for help on using the changeset viewer.