Changeset 6553 for branches/numpy/anuga/coordinate_transforms/redfearn.py

Timestamp:

Mar 19, 2009, 1:43:34 PM (16 years ago)

Author:

rwilson

Message:

Merged trunk into numpy, all tests and validations work.

File:

• branches/numpy/anuga/coordinate_transforms/redfearn.py (modified) (4 diffs)

branches/numpy/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, scale_factor=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
     """
 
@@ -57 +62 @@
     a = 6378137.0                       #Semi major axis
     inverse_flattening = 298.257222101  #1/f
-    K0 = 0.9996                         #Central scale factor    
+    if scale_factor is None:
+        K0 = 0.9996                         #Central scale factor    
+    else:
+        K0 = scale_factor
+    #print 'scale', K0
     zone_width = 6                      #Degrees
 
@@ -138 +147 @@
     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)