source: issm/oecreview/Archive/25834-26739/ISSM-26400-26401.diff

../trunk-jpl/src/m/coordsystems/ll2xy.py

@@ -1 +1 @@
 import numpy as np
 
 
-def ll2xy(lat, lon, sgn=-1, central_meridian=0, standard_parallel=71):
+def ll2xy(lat, lon, sgn, *args):
     """LL2XY - converts lat lon to polar stereographic
 
     Converts from geodetic latitude and longitude to Polar
@@ -15 +15 @@
     - sgn = Sign of latitude    1 : north latitude (default is mer = 45 lat = 70)
                                                    -1 : south latitude (default is mer = 0  lat = 71)
     """
-
     assert sgn == 1 or sgn == -1, 'error: sgn should be either 1 or -1'
-
-    #Get central_meridian and standard_parallel depending on hemisphere
-    if sgn == 1:
-        delta = 45
-        slat = 70
-        print('        ll2xy: creating coordinates in north polar stereographic (Std Latitude: 70N Meridian: 45)')
+         
+         # Get central_meridian and standard_parallel depending on hemisphere
+    if len(args) == 2:
+       delta = args[0]
+       slat = args[1]
+    elif len(args) == 0:
+       if sgn == 1:
+          delta = 45.
+          slat = 70.
+          print('        ll2xy: creating coordinates in north polar stereographic (Std Latitude: 70degN Meridian: 45deg)')
+       elif sgn == -1:
+          delta = 0.
+          slat = 71.
+          print('        ll2xy: creating coordinates in south polar stereographic (Std Latitude: 71degS Meridian: 0deg)')
+       else:
+          raise ValueError('sgn should be either 1 or -1')
     else:
-        delta = central_meridian
-        slat = standard_parallel
-        print('        ll2xy: creating coordinates in south polar stereographic (Std Latitude: 71S Meridian: 0)')
+       raise Exception('bad usage: type "help(ll2xy)" for details')
 
+    # if lat, lon passed as lists, convert to np.arrays
+    if type(lat) != "np.ndarray":
+       lat = np.array(lat)
+    if type(lon) != "np.ndarray":
+       lon = np.array(lon)
+
     # Conversion constant from degrees to radians
     #cde = 57.29577951
     # Radius of the earth in meters
@@ -57 +70 @@
     cnt1 = np.nonzero(latitude >= np.pi / 2.)[0]
 
     if len(cnt1)>0:
-        x[cnt1, 0] = 0.0
-        y[cnt1, 0] = 0.0
+        x[cnt1] = 0.0
+        y[cnt1] = 0.0
+
     return x, y