Changeset 2102 in lmdz_wrf for trunk

Timestamp:

Aug 26, 2018, 1:30:36 AM (7 years ago)

Author:

lfita

Message:

Adding:

• `draw_2Dshad_map': plotting a shadow field with a background map

Location:

trunk/tools

Files:

• drawing.py (modified) (2 diffs)
• drawing_tools.py (modified) (2 diffs)

trunk/tools/drawing.py

@@ -81 +81 @@
 # draw_2D_shad_line_time: plotting a fields with shading and a line with time values
 # draw_2D_shad_time: plotting a fields with shading with time values
+# draw_2Dshad_map: plotting a shadow field with a background map
 # draw_2lines: Fucntion to plot two lines in different axes (x/x2 or y/y2)
 # draw_2lines_time: Function to plot two time-lines in different axes (x/x2 or y/y2)
@@ -10571 +10572 @@
 #draw_WRFeta_levels(fils, vals)
 
+def draw_2Dshad_map(ncfile, values, varn, axfig=None, fig=None):
+    """ plotting a shadow field with a background map
+    draw_2Dshad_map(ncfile, values, varn)
+      ncfile= file to use
+      values=[vnamemap];[vnamefs];[dimvals];[dimxvn];[dimyvn];[dimxyfmt];
+       [colorbarvals];[sminv],[smaxv];[figt];[kindfig];[mapkind];[mapv];[close]
+        [vnamemap]: Name of the variable to be used for map
+        [vnamefs]: Name in the figure of the variable to be shaded
+        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the 
+          variable a given value is required:
+            * [integer]: which value of the dimension
+            * -1: all along the dimension
+            * -9: last value of the dimension
+            * [beg]@[end]@[inc] slice from [beg] to [end] every [inc]
+            * NOTE, no dim name all the dimension size
+        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
+        [dimxyfmt]=[dxs],[dxf],[Ndx],[ordx],[dys],[dyf],[Ndy],[ordx]: format of the values at each axis (or 'auto')
+          [dxs]: style of x-axis ('auto' for 'pretty')
+            'Nfix', values computed at even 'Ndx'
+            'Vfix', values computed at even 'Ndx' increments
+            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
+          [dxf]: format of the labels at the x-axis ('auto' for '%5g')
+          [Ndx]: Number of ticks at the x-axis ('auto' for 5)
+          [ordx]: angle of orientation of ticks at the x-axis ('auto' for horizontal)
+          [dys]: style of y-axis ('auto' for 'pretty')
+          [dyf]: format of the labels at the y-axis ('auto' for '%5g')
+          [Ndy]: Number of ticks at the y-axis ('auto' for 5)
+          [ordy]: angle of orientation of ticks at the y-axis ('auto' for horizontal)
+        [colorbarvals]=[colbarn]#[fmtcolorbar]#[orientation]#[colorbarticks]
+          [colorbarn]: name of the color bar for the values to plot
+          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
+          [orientation]: orientation of the colorbar ('vertical' (default, by 'auto'), 'horizontal')
+          [colorbarticks]: kind of colorbarticks
+            'direct': direct values
+            'time'@[units]|[kind]|[tfmt]|[label] time labels characteristics
+              [units]: units string according to CF conventions ([tunits] since 
+                [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces), 
+                'auto' for minutes!since!1949-12-01!00:00:00
+              [kind]: kind of output
+                'Nval': according to a given number of values as 'Nval',[Nval]
+                'exct': according to an exact time unit as 'exct',[tunit]; 
+                  tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
+                   'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 
+                   'l': milisecond
+              [tfmt]: desired format (combination of 'C'-style values and LaTeX)
+              [label]: label at the graph ('!' for spaces, combination of 'C'-style 
+                values and LaTeX)
+        [smin/axv]: minimum and maximum value for the shading 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)
+        [figt]: title of the figure ('|' for spaces)
+        [kindfig]: kind of file for the figure output (ps, eps, png, pdf, ...)
+        [mapfmt] = [cbarmap]|[mapkind] values for the map 
+          cbarmap: name of the colorbar for the map
+          mapkind: kind of map to use in the plot
+            'direct': values are they come
+            'shadow',[pos]: pseudo-shadding from a given location of the sun
+              [pos]: 'N', 'NW'
+        [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
+        [close]: Whether figure should be finished or not
+      varn= [varsn] name of the variable to plot with shading
+    """
+
+    fname = 'draw_2Dshad_map'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_2Dshad_map.__doc__
+        quit()
+
+    availmapkind = ['direct', 'shadow']
+    availsunpos = ['N', 'NW']
+
+    expectargs = '[vnamemap];[vnamefs];[dimvals];[dimxvn];[dimyvn];[dimxyfmt];' +    \
+      '[colorbarvals];[sminv],[smaxv];[figt];[kindfig];[mapfmt];[mapv];[close]'
+
+    drw.check_arguments(fname,values,expectargs,';')
+
+    vnamemap = values.split(';')[0]
+    vnamesfig = values.split(';')[1]
+    dimvals= values.split(';')[2].replace('|',':')
+    vdimxn = values.split(';')[3]
+    vdimyn = values.split(';')[4]
+    dimxyf = values.split(';')[5]
+    colorbarvals = values.split(';')[6]
+    shadminmax = values.split(';')[7]
+    figtitle = values.split(';')[8].replace('|',' ')
+    figkind = values.split(';')[9]
+    mapfmt = values.split(';')[10].split('|')
+    mapvalue = values.split(';')[11]
+    close = gen.Str_Bool(values.split(';')[12])
+
+    ncfiles = ncfile
+    
+    if not os.path.isfile(ncfiles):
+        print errormsg
+        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
+        quit(-1)    
+
+    objsf = NetCDFFile(ncfiles, 'r')
+    
+    # Map values
+    if not objsf.variables.has_key(vnamemap):
+        print errormsg
+        print '  ' + fname + ": file does not have map variable '" + vnamemap + "' !!"
+        varns = list(objsf.variables.keys())
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+
+    omap = objsf.variables[vnamemap]
+    valsmap, dimsmap = drw.slice_variable(omap, dimvals.replace(',','|'))
+    if len(valsmap.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': map values have to be 2-rank and they are:',         \
+          len(valsmap.shape), '!!'
+        quit(-1)
+
+    dx = valsmap.shape[1]
+    dy = valsmap.shape[0]
+
+    mapkind = mapfmt[1]
+
+    if mapkind == 'direct':
+        mapv = valsmap[:]
+    elif mapkind[0:6] == 'shadow':
+        mapv = np.zeros((dy,dx), dtype=np.float)
+        sunpos = mapkind.split(',')[1]
+        if sunpos == 'N':
+            ## Linear shadow from N
+            mapv[0:dy-1,:] = valsmap[1:dy,:] - valsmap[0:dy-1,:]
+        elif sunpos == 'NW':
+            ## Diagonal shadow from NW
+            for i in range(dx-1):
+                for j in range(dy-1):
+                    if (i-1 >= 0) and (j+1 <= dy-1): 
+                        mapv[j,i] = valsmap[j+1,i-1] - valsmap[j,i]
+        else:
+            print errormsg
+            print '  ' + fname + ": sun location for map sahdding '" + sunpos +      \
+              "' not ready !!"
+            print '    available ones:', availsunpos
+            quit(-1)
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of map '" + mapkind + "' not ready !!"
+        print '    available ones:', availmapkind
+        quit(-1)
+
+    if  not objsf.variables.has_key(varn):
+        print errormsg
+        print '  ' + fname + ': shading file "' + ncfiles +                          \
+          '" does not have variable "' +  varn + '" !!'
+        varns = list(objsf.variables.keys())
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+
+    # Variables' values
+    objvars = objsf.variables[varn]
+    if type(objvars[:]) != type(gen.mamat):
+        print errormsg
+        print '  ' + fname + ": drawing variable is not masked !!"
+        quit(-1)
+
+    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
+
+    dimnamesv = [vdimxn, vdimyn]
+
+    if drw.searchInlist(objvars.ncattrs(),'units'):
+        varunits = objvars.getncattr('units')
+    else:
+        print warnmsg
+        print '  ' + fname + ": variable '" + varn + "' without units!!"
+        varunits = '-'
+
+    if  not objsf.variables.has_key(vdimxn):
+        print errormsg
+        print '  ' + fname + ': shading file "' + ncfiles +                          \
+          '" does not have dimension variable "' +  vdimxn + '" !!'
+        quit(-1)
+    if  not objsf.variables.has_key(vdimyn):
+        print errormsg
+        print '  ' + fname + ': shading file "' + ncfiles +                          \
+          '" does not have dimension variable "' +  vdimyn + '" !!'
+        quit(-1)
+
+    objdimx = objsf.variables[vdimxn]
+    objdimy = objsf.variables[vdimyn]
+    if drw.searchInlist(objdimx.ncattrs(),'units'):
+        odimxu = objdimx.getncattr('units')
+    else:
+        print warnmsg
+        print '  ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
+        odimxu = '-'
+
+    if drw.searchInlist(objdimy.ncattrs(),'units'):
+        odimyu = objdimy.getncattr('units')
+    else:
+        print warnmsg
+        print '  ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
+        odimyu = '-'
+
+    odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions,     \
+      objdimy.dimensions, dimvals.replace(':','|').split(','))
+
+    shading_nx = []
+    if shadminmax.split(',')[0][0:1] != 'S':
+            shading_nx.append(np.float(shadminmax.split(',')[0]))
+    else:
+        shading_nx.append(shadminmax.split(',')[0])
+
+    if shadminmax.split(',')[1][0:1] != 'S':
+        shading_nx.append(np.float(shadminmax.split(',')[1]))
+    else:
+        shading_nx.append(shadminmax.split(',')[1])
+
+    if mapvalue == 'None': mapvalue = None
+
+    colorbarv = colorbarvals.split('#')
+    colorbarvalsS = ','.join(colorbarv[0:3])
+
+    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvalsS, ',')
+    nbarv = np.min(valshad)
+    xbarv = np.max(valshad)
+    rbarv = xbarv - nbarv
+    dbarv = rbarv/5.
+
+    if colorbarv[3] == 'direct':
+        cbarpos = gen.pretty_int(nbarv,xbarv+dbarv,5)
+        cbarlab = []
+        for cpos in cbarpos:
+            cbarlab.append('{:6g}'.format(cpos))
+    elif colorbarv[3][0:4] == 'time':
+        print 'Lluis colorbarv', colorbarv
+        timevals = colorbarv[3].split('@')[1]
+        timeunit = timevals.split('|')[0].replace('!',' ')
+        timekind = timevals.split('|')[1]
+        timefmt = timevals.split('|')[2]
+        timelabel = timevals.split('|')[3].replace('!',' ')
+        cbarpos, cbarlab = drw.CFtimes_plot([nbarv, xbarv], timeunit, timekind, timefmt)
+        vnamesfig = timelabel
+        varunits = '-'
+        fmtcolbar = '%s'
+
+    colormapv = [vnamesfig, varunits, colbarn, fmtcolbar, colbaror, cbarpos, cbarlab]
+
+    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimxyf,',')
+    xaxis = [xstyl, xaxf, Nxax, xaxor]
+    yaxis = [ystyl, yaxf, Nyax, yaxor]
+
+    drw.plot_2Dshad_map(mapv, valshad, mapfmt[0], colormapv, xaxis, yaxis, odimxv,   \
+      odimyv, mapvalue, shading_nx, figtitle, figkind, close)
+
+    return
+
+#ncfile = '/home/lluis/sandbox/get/UBA_ERA-I_1a_2D.nc'
+##values='orog;$conv^{ini}$;lat|-1,lon|-1;lon;lat;auto;rainbow#auto#auto#direct;' +    \
+##  'Srange,Srange;Case1|1a;pdf;binary|shadow,NW;cyl,c;yes'
+#values='orog;$conv^{ini}$;lat|-1,lon|-1;lon;lat;auto;rainbow#auto#auto#' +     \
+#  'time@minutes!since!1949-12-01!00:00:00|exct,12,h|%d$^{%H}$|date!'+          \
+#  '([DD]$^{[HH]}$);Srange,Srange;Case1|1a;pdf;binary|shadow,NW;cyl,c;yes'
+#draw_2Dshad_map(ncfile, values, 'convini', axfig=None, fig=None)
+
 #quit()
 

trunk/tools/drawing_tools.py

@@ -162 +162 @@
 # plot_topo_geogrid_boxes:
 # plot_2D_shadow: Plot of a 2D field with shadow
+# plot_2Dshad_map: Function to plot a map of mask values on top of a map in shadow
 # plot_2D_shadow_time: Plotting a 2D field with one of the axes being time
 # plot_Neighbourghood_evol:Plotting neighbourghood evolution
@@ -13463 +13464 @@
 #fontsize=6,fontsize=6,fontsize=6,fontsize=6
 #prop={'size':legvs[1]}, ncol=legvs[2])prop={'size':8})prop={'size':8})prop={'size':8})
+
+def plot_2Dshad_map(mapvals, values, cbarmap, cbarv, xaxv, yaxv, dimxv, dimyv, mapv, \
+  vs, figtitle, figk, close):
+    """ Function to plot a map of mask values on top of a map in shadow
+      mapvals: 2D values to use for the map
+      values: 2D values to plot
+      cbarmap: cbar name for the behind map
+      cbarv = [varn, varu, cbar, orien, fmt, poscbar, labcbar] list of values of the 
+         colorbar for the values
+        varn: name of the variable in the plot
+        varu: units of the vraiable
+        orien: orientation ('vertical', 'horizontal')
+        fmt: format of ticks labels
+        cbar: cbar name for the values
+        poscbar: position of the ticks of the cbar
+        labcbar: labels of the ticks of the cbar
+      xaxv= list with the x-axis paramteres [style, format, number and orientation]
+      yaxv= list with the y-axis paramteres [style, format, number and orientation]
+      dim[x/y]v: values along x/y axes
+      mapv: values of the map
+      vs: range of values in the figure
+      figtitle: titile of the figure (LaTeX like)
+      figk: kind of output file of the figure (ps, eps, png, pdf, ...)
+      close: whether figure should be closed
+    """
+    fname = 'plot_2Dshad_map'
+
+    if len(mapvals.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': map values have to be 2-rank and they are:',         \
+          len(mapvals.shape), '!!'
+        quit(-1)
+
+    if len(values.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': values have to be 2-rank and they are:',             \
+          len(values.shape), '!!'
+        quit(-1)
+
+    dimxv0 = np.array(dimxv)
+    dimyv0 = np.array(dimyv)
+
+    dxn = dimxv0.min()
+    dxx = dimxv0.max()
+    dyn = dimyv0.min()
+    dyx = dimyv0.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 = 'x'
+    dimyT0 = 'y'
+
+    vsend = graphic_range(vs, values)
+
+    if not mapv is None:
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        lon0 = dimxv.copy()
+        lat0 = dimyv.copy()
+
+        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)
+
+    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))
+
+    plt.rc('text',usetex=True)
+
+    # 2D continuous shaded field
+    plt.pcolormesh(x, y, mapvals, cmap=plt.get_cmap(cbarmap))
+    plt.pcolormesh(x, y, values, cmap=plt.get_cmap(cbarv[2]), vmin=vsend[0], vmax=vsend[1])
+
+    print 'Lluis cbarv:', cbarv
+
+    cbar = plt.colorbar(format=cbarv[3],orientation=cbarv[4], ticks=cbarv[5])
+    if cbarv[4] == 'horizontal':
+        # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib
+        cbar.ax.set_xticklabels(cbarv[6],rotation=90)
+    else:
+        cbar.ax.set_xticklabels(cbarv[6])
+
+    if not mapv is None:
+        if cbarv[2] == '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])
+
+        plt.xlabel('W-E')
+        plt.ylabel('S-N')
+    else:
+        plt.xlabel(dimxT)
+        plt.ylabel(dimyT)
+
+    if mapv is None:
+        ldimxt = list(dimxt)
+        ldimyt = list(dimyt)
+        for val in ldimxt:
+            if val < np.min(limx): ldimxt.remove(val)
+            if val > np.max(limx): ldimxt.remove(val)
+        for val in ldimyt:
+            if val < np.min(limy): ldimyt.remove(val)
+            if val > np.max(limy): ldimyt.remove(val)
+        dimxt = np.array(ldimxt)
+        dimyt = np.array(ldimyt)
+        dxt = gen.interpolate_locs(dimxt,dimxv[0,:],'lin')
+        dyt = gen.interpolate_locs(dimyt,dimyv[:,0],'lin')
+        dimxt = dxt/varsv.shape[1]
+        dimyt = dyt/varsv.shape[0]
+
+        # Using output from transform
+        gxn = (gmaxmin[0]-xmin)/(xmax-xmin)
+        gxx = (gmaxmin[1]-xmin)/(xmax-xmin)
+        gyn = (gmaxmin[2]-ymin)/(ymax-ymin)
+        gyx = (gmaxmin[3]-ymin)/(ymax-ymin)
+
+        #gxn = (limx[0]-xmin)/(xmax-xmin)
+        #gxx = (limx[1]-xmin)/(xmax-xmin)
+        #gyn = (limy[0]-ymin)/(ymax-ymin)
+        #gyx = (limy[1]-ymin)/(ymax-ymin)
+
+        # Are axis flipped?
+        signxl = np.abs(limx[0]-limx[1])/(limx[0]-limx[1])
+        signxt = np.abs(ldimxt[0]-ldimxt[1])/(ldimxt[0]-ldimxt[1])
+        signyl = np.abs(limy[0]-limy[1])/(limy[0]-limy[1])
+        signyt = np.abs(ldimyt[0]-ldimyt[1])/(ldimyt[0]-ldimyt[1])
+
+        if signxl != signxt:
+            plt.xticks(1.-dimxt, dimxl, rotation=xaxv[3])
+        else:
+            plt.xticks(dimxt, dimxl, rotation=xaxv[3])
+        if signyl != signyt:
+            plt.yticks(1.-dimyt, dimyl, rotation=yaxv[3])
+        else:
+            plt.yticks(dimyt, dimyl, rotation=yaxv[3])
+
+        # Only if gmaxmin is used to compute limits of the figure !!
+        if signxl != signxt:
+            newgxn = gxx
+            newgxx = gxn
+        else:
+            newgxn = gxn
+            newgxx = gxx
+        if signyl != signyt:
+            newgyn = gyx
+            newgyx = gyn
+        else:
+            newgyn = gyn
+            newgyx = gyx
+
+        plt.axis([newgxn, newgxx, newgyn, newgyx])
+        # otherwise
+        #plt.axis([gxn, gxx, gyn, gyx])
+
+# units labels
+    cbar.set_label(gen.latex_text(cbarv[0]) + ' (' + units_lunits(cbarv[1]) + ')')
+    
+    plt.title(gen.latex_text(figtitle))
+
+    figname = '2Dshad_map'
+    output_kind(figk, figname, close)
+
+    return
+