Index: trunk/tools/drawing_tools.py
===================================================================
--- trunk/tools/drawing_tools.py	(revision 1802)
+++ trunk/tools/drawing_tools.py	(revision 1872)
@@ -10843,2 +10843,257 @@
 
     return
+
+def multi_plot_2D_shadow(lvarsv, dimxv, dimyv, dimxu, dimyu, xaxv, yaxv, dimn,       \
+  cbarv, vs, uts, kfig, reva, mapv, labs, Ncol, Nrow, ftit, ifclose):
+    """ plot of multiple 2D shadowing variables sharing colorbar
+      lvarsv= list of 2D values to plot with shading
+      dim[x/y]v = list of values at the axes of x and y
+      dim[x/y]u = units at the axes of x and y
+      xaxv= list with the x-axis paramteres [style, format, number and orientation]
+      yaxv= list with the y-axis paramteres [style, format, number and orientation]
+      dimn= dimension names to plot
+      cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror]
+        colorbar: name of the color bar to use
+        cbarfmt: format of the numbers in the colorbar
+        cbaror: orientation of the colorbar
+      vs= minmum and maximum values to plot in shadow or:
+        'Srange': for full range
+        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
+        'Saroundminmax@val': for min*val,max*val
+        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
+          percentile_(100-val)-median)
+        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
+        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
+        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
+           percentile_(100-val)-median)
+      uts= units of the variable to shadow
+      vtit= title of the variable
+      kfig= kind of figure (jpg, pdf, png)
+      reva= ('|' for combination)
+        * 'transpose': reverse the axes (x-->y, y-->x)
+        * 'flip'@[x/y]: flip the axis x or y
+      mapv= map characteristics: [proj],[res]
+        see full documentation: http://matplotlib.org/basemap/
+        [proj]: projection
+          * 'cyl', cilindric
+          * 'lcc', lambert conformal
+        [res]: resolution:
+          * 'c', crude
+          * 'l', low
+          * 'i', intermediate
+          * 'h', high
+          * 'f', full
+      labs: list of consecutive labels for each figure
+      ftit: global title of the figure
+      Ncol: number of column panels
+      Nrow: number of rows of panels
+      ifclose= boolean value whether figure should be close (finish) or not
+    """
+##    import matplotlib as mpl
+##    mpl.use('Agg')
+##    import matplotlib.pyplot as plt
+    fname = 'multi_plot_2D_shadow'
+
+#    print dimyv[73,21]
+#    dimyv[73,21] = -dimyv[73,21]
+#    print 'Lluis dimsv: ',np.min(dimxv), np.max(dimxv), ':', np.min(dimyv), np.max(dimyv)
+
+    if type(lvarsv) == type('S') and varsv == 'h':
+        print fname + '_____________________________________________________________'
+        print multi_plot_2D_shadow.__doc__
+        quit()
+ 
+    # Getting the right lon values for plotting
+    if mapv is not None:
+        dimxv = np.where(dimxv > 180., dimxv-360., dimxv)
+
+    # Axis ticks
+    # Usually axis > x must be the lon, thus...
+    dimxv0 = dimxv.copy()
+    dimyv0 = dimyv.copy()
+
+    dxn = dimxv.min()
+    dxx = dimxv.max()
+    dyn = dimyv.min()
+    dyx = dimyv.max()
+
+    if xaxv[0] == 'pretty':
+        dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2]))
+    elif xaxv[0] == 'Nfix':
+        dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2]))
+    elif xaxv[0] == 'Vfix':
+        dimxt0 = np.arange(0,dxx,xaxv[2])
+    if yaxv[0] == 'pretty':
+        dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2]))
+    elif yaxv[0] == 'Nfix':
+        dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2]))
+    elif yaxv[0] == 'Vfix':
+        dimyt0 = np.arange(0,dyx,yaxv[2])
+
+    dimxl0 = []
+    for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], style=xaxv[1]))
+    dimyl0 = []
+    for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1]))
+
+    dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')'
+    dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')'
+
+    if mapv is not None:
+        pixkind = 'data'
+    else:
+        # No following data values
+        dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0))
+        dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0))
+        pixkind = 'fixpixel'
+
+    if reva is not None:
+        varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT =              \
+          transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0,     \
+          dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0)
+    else:
+        dimxv = dimxv0
+        dimyv = dimyv0
+        dimxt = dimxt0
+        dimyt = dimyt0
+        dimxl = dimxl0
+        dimyl = dimyl0
+        dimxT = dimxT0
+        dimyT = dimyT0
+
+    if len(dimxv[:].shape) == 3:
+        xdims = '1,2'
+    elif len(dimxv[:].shape) == 2:
+        xdims = '0,1'
+    elif len(dimxv[:].shape) == 1:
+        xdims = '0'
+
+    if len(dimyv[:].shape) == 3:
+        ydims = '1,2'
+    elif len(dimyv[:].shape) == 2:
+        ydims = '0,1'
+    elif len(dimyv[:].shape) == 1:
+        ydims = '0'
+
+    lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)
+
+    if not mapv is None:
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        dx = lon0.shape[1]
+        dy = lat0.shape[0]
+
+        nlon = np.min(lon0)
+        xlon = np.max(lon0)
+        nlat = np.min(lat0)
+        xlat = np.max(lat0)
+
+        lon2 = lon0[dy/2,dx/2]
+        lat2 = lat0[dy/2,dx/2]
+
+        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
+          xlon, ',', xlat
+
+        if map_proj == 'cyl':
+            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
+              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        elif map_proj == 'lcc':
+            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
+              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        else:
+            print errormsg
+            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
+            print '    available: cyl, lcc'
+            quit(-1)
+ 
+        x,y = m(lon0,lat0)
+        if cbarv[0] == 'gist_gray':
+            m.drawcoastlines(color="red")
+        else:
+            m.drawcoastlines()
+        meridians = gen.pretty_int(nlon,xlon,xaxv[2])
+        m.drawmeridians(meridians,labels=[True,False,False,True])
+        parallels = gen.pretty_int(nlat,xlat,yaxv[2])
+        m.drawparallels(parallels,labels=[False,True,True,False])
+
+    else:
+        # No following data values
+        x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv))
+        y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv))
+
+# Changing limits of the colors
+    vsend = graphic_range(vs,lvarsv)
+
+    # Number of panels
+    Npanels = len(lvarsv)
+    if Ncol*Nrow != Npanels:
+        print errormsg
+        print '  ' + fname + ': no coincident number of panels:', Npanels,           \
+          'with provied distribution of them:', Ncol, Nrow
+
+    plt.rc('text', usetex=True)
+    fig, axes = plt.subplots(Nrow, Ncol)
+
+    irow=0
+    icol=0
+    for iv in range(Npanels):
+        varsv = lvarsv[iv]
+
+        if len(varsv.shape) != 2:
+            print errormsg
+            print '  ' + fname + ': wrong variable shadow rank:', varsv.shape,       \
+              'is has to be 2D!!'
+            quit(-1)
+        axes[irow,icol] = plt.subplot(Nrow ,Ncol, iv+1)
+        im = plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], \
+          vmax=vsend[1])
+    
+        if not mapv is None:
+            m.drawcoastlines()
+            m.drawmeridians(meridians,labels=[True,False,False,True])
+            m.drawparallels(parallels,labels=[False,True,True,False])
+            if irow == Nrow -1: plt.xlabel('W-E')
+            if icol == 0: plt.ylabel('S-N')
+        else:
+            plt.xlabel(dimxT)
+            plt.ylabel(dimyT)
+            plt.xticks(dimxt, dimxl, rotation=xaxv[3])
+            plt.yticks(dimyt, dimyl, rotation=yaxv[3])
+
+        plt.axis([x.min(), x.max(), y.min(), y.max()])
+
+        graphtit = gen.latex_text(labs[iv].replace('!', ' '))
+        plt.title(gen.latex_text(graphtit))
+        icol = icol + 1
+        if icol > Ncol-1:
+            irow = irow + 1
+            icol = 0
+
+    fig.subplots_adjust(right=0.8, hspace=0.4)
+    cbar_ax = fig.add_axes([0.86, 0.25, 0.03, 0.50])
+
+    if cbarv[2] == 'horizontal':
+        cbar = fig.colorbar(im, cax=cbar_ax, format=cbarv[1],orientation=cbarv[2])
+        # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib
+        ticklabels= cbar.ax.get_xticklabels()
+        Nticks = len(ticklabels)
+        ticklabs = []
+        for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text())
+        cbar.ax.set_xticklabels(ticklabs,rotation=90)
+    else:
+        cbar = fig.colorbar(im, cax=cbar_ax, format=cbarv[1],orientation=cbarv[2])
+
+    #fig.colorbar(im, cax=cbar_ax)
+# units labels
+    ulab = units_lunits(uts)
+    cbar_ax.set_label(ulab)
+    plt.annotate(ulab, xy=(0.96,0.5), xycoords='figure fraction', rotation=90)
+    cbar_ax.tick_params(labelsize=8)
+
+    plt.suptitle(ftit)
+    figname = 'multi_2Dfields_shadow'
+
+    output_kind(kfig, figname, ifclose)
+
+    return
+
