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

Timestamp:

Feb 24, 2009, 6:13:06 PM (16 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:

• anuga_core/source/anuga/coordinate_transforms/redfearn.py (modified) (10 diffs)

anuga_core/source/anuga/coordinate_transforms/redfearn.py

@@ -37 +37 @@
     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
 
@@ -47 +48 @@
     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
     """
 
@@ -54 +60 @@
 
 
-    #GDA Specifications
+    # GDA Specifications
     a = 6378137.0                       #Semi major axis
     inverse_flattening = 298.257222101  #1/f
@@ -73 +79 @@
         
     
-    #Derived constants
+    # Derived constants
     f = 1.0/inverse_flattening
     b = a*(1-f)       #Semi minor axis
@@ -84 +90 @@
     e6 = e2*e4
 
-    #Foot point latitude
+    # Foot point latitude
     n = (a-b)/(a+b) #Same as e2 - why ?
     n2 = n*n
@@ -114 +120 @@
     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
@@ -122 +128 @@
     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)
@@ -138 +141 @@
     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)
@@ -153 +164 @@
     omega8 = omega4*omega4
      
-    #Northing
+    # Northing
     term1 = nu*sinphi*cosphi*omega2/2  
     term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24
@@ -162 +173 @@
     northing = false_northing + K0*(m + term1 + term2 + term3 + term4)
 
-    #Easting
+    # Easting
     term1 = nu*omega*cosphi
     term2 = nu*cosphi3*(psi-t2)*omega3/6