Changeset 1872 in lmdz_wrf for trunk/tools/drawing_tools.py


Ignore:
Timestamp:
Mar 27, 2018, 5:50:52 PM (7 years ago)
Author:
lfita
Message:

Adding:

  • `draw_multi_2D_shad': plotting multiple 2D fields with same projection with shading and sharing colorbar
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/tools/drawing_tools.py

    r1802 r1872  
    1084310843
    1084410844    return
     10845
     10846def multi_plot_2D_shadow(lvarsv, dimxv, dimyv, dimxu, dimyu, xaxv, yaxv, dimn,       \
     10847  cbarv, vs, uts, kfig, reva, mapv, labs, Ncol, Nrow, ftit, ifclose):
     10848    """ plot of multiple 2D shadowing variables sharing colorbar
     10849      lvarsv= list of 2D values to plot with shading
     10850      dim[x/y]v = list of values at the axes of x and y
     10851      dim[x/y]u = units at the axes of x and y
     10852      xaxv= list with the x-axis paramteres [style, format, number and orientation]
     10853      yaxv= list with the y-axis paramteres [style, format, number and orientation]
     10854      dimn= dimension names to plot
     10855      cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror]
     10856        colorbar: name of the color bar to use
     10857        cbarfmt: format of the numbers in the colorbar
     10858        cbaror: orientation of the colorbar
     10859      vs= minmum and maximum values to plot in shadow or:
     10860        'Srange': for full range
     10861        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
     10862        'Saroundminmax@val': for min*val,max*val
     10863        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
     10864          percentile_(100-val)-median)
     10865        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
     10866        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
     10867        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
     10868           percentile_(100-val)-median)
     10869      uts= units of the variable to shadow
     10870      vtit= title of the variable
     10871      kfig= kind of figure (jpg, pdf, png)
     10872      reva= ('|' for combination)
     10873        * 'transpose': reverse the axes (x-->y, y-->x)
     10874        * 'flip'@[x/y]: flip the axis x or y
     10875      mapv= map characteristics: [proj],[res]
     10876        see full documentation: http://matplotlib.org/basemap/
     10877        [proj]: projection
     10878          * 'cyl', cilindric
     10879          * 'lcc', lambert conformal
     10880        [res]: resolution:
     10881          * 'c', crude
     10882          * 'l', low
     10883          * 'i', intermediate
     10884          * 'h', high
     10885          * 'f', full
     10886      labs: list of consecutive labels for each figure
     10887      ftit: global title of the figure
     10888      Ncol: number of column panels
     10889      Nrow: number of rows of panels
     10890      ifclose= boolean value whether figure should be close (finish) or not
     10891    """
     10892##    import matplotlib as mpl
     10893##    mpl.use('Agg')
     10894##    import matplotlib.pyplot as plt
     10895    fname = 'multi_plot_2D_shadow'
     10896
     10897#    print dimyv[73,21]
     10898#    dimyv[73,21] = -dimyv[73,21]
     10899#    print 'Lluis dimsv: ',np.min(dimxv), np.max(dimxv), ':', np.min(dimyv), np.max(dimyv)
     10900
     10901    if type(lvarsv) == type('S') and varsv == 'h':
     10902        print fname + '_____________________________________________________________'
     10903        print multi_plot_2D_shadow.__doc__
     10904        quit()
     10905 
     10906    # Getting the right lon values for plotting
     10907    if mapv is not None:
     10908        dimxv = np.where(dimxv > 180., dimxv-360., dimxv)
     10909
     10910    # Axis ticks
     10911    # Usually axis > x must be the lon, thus...
     10912    dimxv0 = dimxv.copy()
     10913    dimyv0 = dimyv.copy()
     10914
     10915    dxn = dimxv.min()
     10916    dxx = dimxv.max()
     10917    dyn = dimyv.min()
     10918    dyx = dimyv.max()
     10919
     10920    if xaxv[0] == 'pretty':
     10921        dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2]))
     10922    elif xaxv[0] == 'Nfix':
     10923        dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2]))
     10924    elif xaxv[0] == 'Vfix':
     10925        dimxt0 = np.arange(0,dxx,xaxv[2])
     10926    if yaxv[0] == 'pretty':
     10927        dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2]))
     10928    elif yaxv[0] == 'Nfix':
     10929        dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2]))
     10930    elif yaxv[0] == 'Vfix':
     10931        dimyt0 = np.arange(0,dyx,yaxv[2])
     10932
     10933    dimxl0 = []
     10934    for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], style=xaxv[1]))
     10935    dimyl0 = []
     10936    for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1]))
     10937
     10938    dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')'
     10939    dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')'
     10940
     10941    if mapv is not None:
     10942        pixkind = 'data'
     10943    else:
     10944        # No following data values
     10945        dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0))
     10946        dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0))
     10947        pixkind = 'fixpixel'
     10948
     10949    if reva is not None:
     10950        varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT =              \
     10951          transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0,     \
     10952          dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0)
     10953    else:
     10954        dimxv = dimxv0
     10955        dimyv = dimyv0
     10956        dimxt = dimxt0
     10957        dimyt = dimyt0
     10958        dimxl = dimxl0
     10959        dimyl = dimyl0
     10960        dimxT = dimxT0
     10961        dimyT = dimyT0
     10962
     10963    if len(dimxv[:].shape) == 3:
     10964        xdims = '1,2'
     10965    elif len(dimxv[:].shape) == 2:
     10966        xdims = '0,1'
     10967    elif len(dimxv[:].shape) == 1:
     10968        xdims = '0'
     10969
     10970    if len(dimyv[:].shape) == 3:
     10971        ydims = '1,2'
     10972    elif len(dimyv[:].shape) == 2:
     10973        ydims = '0,1'
     10974    elif len(dimyv[:].shape) == 1:
     10975        ydims = '0'
     10976
     10977    lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)
     10978
     10979    if not mapv is None:
     10980        map_proj=mapv.split(',')[0]
     10981        map_res=mapv.split(',')[1]
     10982
     10983        dx = lon0.shape[1]
     10984        dy = lat0.shape[0]
     10985
     10986        nlon = np.min(lon0)
     10987        xlon = np.max(lon0)
     10988        nlat = np.min(lat0)
     10989        xlat = np.max(lat0)
     10990
     10991        lon2 = lon0[dy/2,dx/2]
     10992        lat2 = lat0[dy/2,dx/2]
     10993
     10994        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
     10995          xlon, ',', xlat
     10996
     10997        if map_proj == 'cyl':
     10998            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
     10999              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
     11000        elif map_proj == 'lcc':
     11001            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
     11002              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
     11003        else:
     11004            print errormsg
     11005            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
     11006            print '    available: cyl, lcc'
     11007            quit(-1)
     11008 
     11009        x,y = m(lon0,lat0)
     11010        if cbarv[0] == 'gist_gray':
     11011            m.drawcoastlines(color="red")
     11012        else:
     11013            m.drawcoastlines()
     11014        meridians = gen.pretty_int(nlon,xlon,xaxv[2])
     11015        m.drawmeridians(meridians,labels=[True,False,False,True])
     11016        parallels = gen.pretty_int(nlat,xlat,yaxv[2])
     11017        m.drawparallels(parallels,labels=[False,True,True,False])
     11018
     11019    else:
     11020        # No following data values
     11021        x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv))
     11022        y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv))
     11023
     11024# Changing limits of the colors
     11025    vsend = graphic_range(vs,lvarsv)
     11026
     11027    # Number of panels
     11028    Npanels = len(lvarsv)
     11029    if Ncol*Nrow != Npanels:
     11030        print errormsg
     11031        print '  ' + fname + ': no coincident number of panels:', Npanels,           \
     11032          'with provied distribution of them:', Ncol, Nrow
     11033
     11034    plt.rc('text', usetex=True)
     11035    fig, axes = plt.subplots(Nrow, Ncol)
     11036
     11037    irow=0
     11038    icol=0
     11039    for iv in range(Npanels):
     11040        varsv = lvarsv[iv]
     11041
     11042        if len(varsv.shape) != 2:
     11043            print errormsg
     11044            print '  ' + fname + ': wrong variable shadow rank:', varsv.shape,       \
     11045              'is has to be 2D!!'
     11046            quit(-1)
     11047        axes[irow,icol] = plt.subplot(Nrow ,Ncol, iv+1)
     11048        im = plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], \
     11049          vmax=vsend[1])
     11050   
     11051        if not mapv is None:
     11052            m.drawcoastlines()
     11053            m.drawmeridians(meridians,labels=[True,False,False,True])
     11054            m.drawparallels(parallels,labels=[False,True,True,False])
     11055            if irow == Nrow -1: plt.xlabel('W-E')
     11056            if icol == 0: plt.ylabel('S-N')
     11057        else:
     11058            plt.xlabel(dimxT)
     11059            plt.ylabel(dimyT)
     11060            plt.xticks(dimxt, dimxl, rotation=xaxv[3])
     11061            plt.yticks(dimyt, dimyl, rotation=yaxv[3])
     11062
     11063        plt.axis([x.min(), x.max(), y.min(), y.max()])
     11064
     11065        graphtit = gen.latex_text(labs[iv].replace('!', ' '))
     11066        plt.title(gen.latex_text(graphtit))
     11067        icol = icol + 1
     11068        if icol > Ncol-1:
     11069            irow = irow + 1
     11070            icol = 0
     11071
     11072    fig.subplots_adjust(right=0.8, hspace=0.4)
     11073    cbar_ax = fig.add_axes([0.86, 0.25, 0.03, 0.50])
     11074
     11075    if cbarv[2] == 'horizontal':
     11076        cbar = fig.colorbar(im, cax=cbar_ax, format=cbarv[1],orientation=cbarv[2])
     11077        # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib
     11078        ticklabels= cbar.ax.get_xticklabels()
     11079        Nticks = len(ticklabels)
     11080        ticklabs = []
     11081        for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text())
     11082        cbar.ax.set_xticklabels(ticklabs,rotation=90)
     11083    else:
     11084        cbar = fig.colorbar(im, cax=cbar_ax, format=cbarv[1],orientation=cbarv[2])
     11085
     11086    #fig.colorbar(im, cax=cbar_ax)
     11087# units labels
     11088    ulab = units_lunits(uts)
     11089    cbar_ax.set_label(ulab)
     11090    plt.annotate(ulab, xy=(0.96,0.5), xycoords='figure fraction', rotation=90)
     11091    cbar_ax.tick_params(labelsize=8)
     11092
     11093    plt.suptitle(ftit)
     11094    figname = 'multi_2Dfields_shadow'
     11095
     11096    output_kind(kfig, figname, ifclose)
     11097
     11098    return
     11099
Note: See TracChangeset for help on using the changeset viewer.