Update to Mapping

Author: antonio-navarra

zapata/__init__.py

@@ -27,6 +27,7 @@
 **koopman** :
     Routines for Koopman decomposition and eigenvalues handling
 
+    Version 2.1
 """
 import warnings
 warnings.simplefilter('ignore')

zapata/colormap.py

@@ -29,7 +29,7 @@
 from lxml import etree
 
 homedir = os.path.expanduser("~")
-COLPATHDIR = homedir + '/Dropbox (CMCC)/GitHub/Zapata/zapata/SciVis_colormaps'
+COLPATHDIR = homedir + '/Dropbox (CMCC)/ZapataLibrary/Zapata/zapata/SciVis_colormaps'
 
 def make_cmap(xml,colpath=COLPATHDIR):
     '''

zapata/computation.py

@@ -196,6 +196,8 @@ def anomaly(var,option='anom',freq='month'):
     anom :  xarray
 
     """
+    # Correction for very small standard deviations
+    TINY = 1.0e-14
 
     frequency = 'time.' + freq
     if option == 'deviation':
@@ -211,7 +213,7 @@ def anomaly(var,option='anom',freq='month'):
         anom = xr.apply_ufunc(lambda x, m, s: (x - m) / s,
                 var.groupby(frequency),
                 clim,
-                climstd )
+                climstd + TINY)
     else:
         print(' Wrong option in `anomaly` {}'.format(option))
         raise SystemExit
@@ -275,8 +277,7 @@ def __init__(
         self.A = X.stack(z=dims).transpose()
         self._ntime = len(X.time.data)
         self._npoints = len(X.stack(z=dims).z.data)
-        print(' Created mathematical matrix A, \n \
-                stacked along dimensions {} '.format(dims))
+        print('Created data Matrix X, stacked along dimensions {} '.format(dims))
         if self._opt == 'DropNaN':
             self.A = self.A.dropna(dim='z')
             print(' Creating Matrix with Drop NaN values')
@@ -289,7 +290,7 @@ def __call__(self, v ):
 
     def __repr__(self):
         '''  Printing Information '''
-        print(' \n Math Data Matrix \n {} \n'.format(self.A))
+        print(f' \n Math Data Matrix \n {self.A}\n')
         print(f' Shape of A numpy array {self.A.shape}')
         return  '\n'
 

zapata/mapping.py

@@ -24,7 +24,7 @@
 '''
 import math as mp
 from operator import is_
-import sys
+from collections import namedtuple
 import cartopy.crs as car
 import cartopy.util as utl
 
@@ -151,27 +151,31 @@ def choose_projection(proview):
         Projection object   
     '''
 
-#Defaults for Satellite
+    #Defaults for Satellite
     default_satellite = {'centlon': 0.0, 'centlat': 0.0, 'z': 35785831}
 
     # This fixes the projection for the mapping: 
     #      Data Projection is then fixed in mapping routine
     # Select Projection parameters
     if isinstance(proview,str):
-        if   proview == 'Pacific':
-            projection = car.PlateCarree(central_longitude=180.)
-        elif proview == 'Atlantic':
-            projection = car.PlateCarree(central_longitude=0.)
-        elif proview == 'NHStereoEurope':
-            projection = car.NorthPolarStereo(central_longitude=0.)
-        elif proview == 'NHStereoAmerica':
-            projection = car.NorthPolarStereo(central_longitude=-90.)
-        elif proview == 'SHStereoAfrica':
-            projection = car.SouthPolarStereo(central_longitude=90.)
-        elif proview == 'RobAtlantic':
-            projection = car.Robinson(central_longitude=0.)
-        elif proview == 'RobPacific':
-            projection = car.Robinson(central_longitude=180.)
+        match proview:
+            case 'Pacific':
+                projection = car.PlateCarree(central_longitude=180.)
+            case 'Atlantic':
+                projection = car.PlateCarree(central_longitude=0.)
+            case 'NHStereoEurope':
+                projection = car.NorthPolarStereo(central_longitude=0.)
+            case 'NHStereoAmerica':
+                projection = car.NorthPolarStereo(central_longitude=-90.)
+            case 'SHStereoAfrica':
+                projection = car.SouthPolarStereo(central_longitude=90.)
+            case 'RobAtlantic':
+                projection = car.Robinson(central_longitude=0.)
+            case 'RobPacific':
+                projection = car.Robinson(central_longitude=180.)
+            case _:
+                print(f' Error in init_figure projection {proview}')
+                print(f'It should be one of the following: Pacific, Atlantic, NHStereoEurope, NHStereoAmerica, SHStereoAfrica, RobAtlantic, RobPacific')
     elif isinstance(proview,dict):
         temp = { **default_satellite,**proview}
         print(temp['z'])
@@ -202,7 +206,8 @@ def xmap(data, cont, prol, ax=None, fill=True,contour=True, \
                     cyclic=False,\
                     colorbar=False,cmap='coolwarm',\
                     coasts=True,color_land='lightgray',\
-                    quiet=True,\
+                    quiet=True,
+                    custom=False,\
                     label_style={}):
     """
     Mapping function based on cartopy and matplotlib.
@@ -273,6 +278,8 @@ def xmap(data, cont, prol, ax=None, fill=True,contour=True, \
         The projection used is the geographical projection `pro`
     quiet:
         (False)  Suppress all output
+    custom:
+        (False)  Return a named tuple with full detail on the plot, including handles, gl etc...
     
     
     Returns
@@ -360,7 +367,11 @@ def xmap(data, cont, prol, ax=None, fill=True,contour=True, \
             xlim=xlimit
 
         if ylimit is  None:
-            ylim = [np.amin(data.lat.values),np.amax(data.lat.values)]
+            match this:
+                case 'NorthPolarStereo':
+                    ylim = [20,90]
+                case 'SouthPolarStereo':
+                    ylim = [-90,-20]
         else:
             ylim=ylimit
         print(' Plotting with x limits {}  '.format(xlim)  ) 
@@ -386,12 +397,12 @@ def xmap(data, cont, prol, ax=None, fill=True,contour=True, \
 
     cc = choose_contour(vmin,vmax,cont)
     if not quiet:
+        print(title)
         print(f'Contouring from {vmin.data} to {vmax.data} with {cc} contours')
     if cc is None:
         print(f'******** Constant field ***********')
         return None
 
-
     handles = dict()
     if fill:
         handles['filled'] = ax.contourf(data.lon, data.lat, data, levels=cc, cmap=cmap, transform=data_proj)
@@ -434,7 +445,12 @@ def xmap(data, cont, prol, ax=None, fill=True,contour=True, \
         cax = divider.append_axes('right',size="2.5%", pad=0.2, axes_class=plt.Axes)
         ax.get_figure().colorbar(handles['filled'], cax=cax,orientation='vertical')
 
-    return handles  
+    # Choose output
+    Output = namedtuple('Output', 'handles gl' )
+    if custom:
+        return Output(handles,gl)
+    else:
+        return handles
 
 @d.get_sections(base='zonal_plot')
 @d.dedent
@@ -665,7 +681,7 @@ def choose_contour(vmin,vmax,cont):
     """
     if len(cont) == 3:
         print("Setting Fixed Contours")
-        cc=np.arange(cont[0]-cont[2],cont[1]+2*cont[2],cont[2])
+        cc=np.arange(cont[0],cont[1]+cont[2],cont[2])
         print(' Contouring from ', cont[0], '  to', cont[1],' with interval  ',cont[2])  
     elif len(cont)> 3:
         cc=cont