Changeset 539 in lmdz_wrf for trunk/tools

Timestamp:

Jun 29, 2015, 7:44:30 PM (9 years ago)

Author:

lfita

Message:

Adding a not tested versino of the barb plotting

Location:

trunk/tools

Files:

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

trunk/tools/drawing.py

@@ -30 +30 @@
 namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time',                 \
   'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line',                \
-  'draw_2D_shad_line_time', 'draw_timeSeries', 'draw_topo_geogrid',                  \
+  'draw_2D_shad_line_time', 'draw_barbs', 'draw_timeSeries', 'draw_topo_geogrid',    \
   'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
   'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', 'list_graphics',      \
@@ -1175 +1175 @@
     return
 
+def draw_barbs(ncfiles, values, varns):
+    """ Function to plot wind barbs
+      values= [X/Y/Z/T]|[dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:[mapvalues]:
+        [gtit]:[kindfig]:[figuren]
+        [dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
+          [dimname]: name of the dimension in the file
+          [vardimname]: name of the variable with the values for the dimension in the file
+          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
+          No value takes all the range of the dimension
+        [vecvals]= [frequency],[color],[length]
+          [freqv]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points; 
+            'auto', computed automatically to have 20 vectors along each axis)
+          [colorv]: color of the vectors ('None', for 'red')
+          [lengthv]: length of the wind barbs ('None', for 9)
+        [windlabs]= [windname],[windunits]
+          [windname]: name of the wind variable in the graph
+          [windunits]: units of the wind variable in the graph ('None', for the value in the file)
+        [mapvalues]= 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
+        gtit= title of the graph ('|', for spaces)
+        kindfig= kind of figure
+        figuren= name of the figure
+      ncfile= file to use
+      varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
+    """
+    fname = 'draw_barbs'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_barbs.__doc__
+        quit()
+
+    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:' +  \
+      '[mapvalues]:[gtit]:[kindfig]:[figuren]'
+ 
+    drw.check_arguments(fname,len(expectargs.split(':')),values,':',                 \
+      expectargs.split(':'))
+
+    dimvals = values.split(':')[0]
+    vecvals = values.split(':')[1]
+    windlabs = values.split(':')[2]
+    mapvalues = values.split(':')[3]
+    gtit = values.split(':')[4]
+    kindfig = values.split(':')[5]
+    figuren = values.split(':')[6]
+
+    of = NetCDFFile(ncfile,'r')
+
+    dims = {}
+    for dimv in dimvals.split(','):
+        dns = dimv.split('|')
+        dims[dns[0]] = [dns[1], dns[2]]
+
+    varNs = []
+    for dn in dims.keys():
+        if dn == 'X':
+            varNs.append(dims[dn][2])
+        elif dn == 'Y':
+            varNs.append(dims[dn][2])
+
+    ivar = 0
+    for wvar in varns.split(','):
+        if not drw.searchInlist(of.variables.keys(), wvar):
+            print errormsg
+            print '  ' + fname + ": file does not have variable '" + wvar + "' !!"
+            quit(-1)
+        if ivar == 0:
+            varNs.append(wvar)
+        else:
+            varNs.append(wvar)
+
+    ivar = 0
+    for varN in varNs:
+        varslice = []
+
+        ovarN = of.variables[varN]
+        vard = ovarN.dimensions.keys()
+        for vdn in vard:
+            for dd in dims.keys():
+                if dd == vdn:
+                    if dims[dd][1].find('@') != -1:
+                        rvals = dims[dd][1].replace('@')
+                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
+                    elif dims[dd][1] == '-1':
+                        varslice.append(slice(0,len(of.dimensions[dims[dd]])))
+                    break
+                else:
+                    varslice.append(slice(0,len(of.dimensions[dims[dd]])))
+
+        if ivar == 0:
+            lonvals0 = ovarN[tuple(varslice)]
+        elif ivar == 1:
+            latvals0 = ovarN[tuple(varslice)]
+        elif ivar == 2:
+            uwvals = ovarN[tuple(varslice)]
+        elif ivar == 3:
+            vwvals = ovarN[tuple(varslice)]
+
+    if len(lonvals0.shape) == 1:
+        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
+
+# Vecor values
+    freqv = vecvals.split(',')[0]
+    colorv = vecvals.split(',')[1]
+    lengthv = vecvals.split(',')[2]
+
+# Vector labels
+    windname = windlabels.split(',')[0]
+    windunits = windlabels.split(',')[1]
+
+    drw.plot_barbs(lonvals, latvals, uwvals, vwvals, freqv, colorv, lengthv, 
+      windname, windunits, mapvalues, gtit, kindfig, figuren)
+
+    return
+ 
 def draw_topo_geogrid(ncfile, values):
     """ plotting geo_em.d[nn].nc topography from WPS files
@@ -2518 +2642 @@
 elif oper == 'draw_2D_shad_line_time':
     draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
+elif oper == 'draw_barbs':
+    draw_barbs(opts.ncfile, opts.values, opts.varname)
 elif oper == 'draw_Neighbourghood_evol':
     draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -68 +68 @@
       if x == nameFind:
         return True
+        break
     return False
 
@@ -5170 +5171 @@
 
     return
+
+def plot_barbs(xvals,yvals,vecfreq,veccolor,veclength,windn,wuts,mapv,graphtit,kfig,figname):
+    """ Function to plot wind barbs
+      xvals= values for the 'x-axis'
+      yvals= values for the 'y-axis'
+      vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 
+        'auto', computed automatically to have 20 vectors along each axis)
+      veccolor= color of the vectors (None, for 'red')
+      veclength= length of the wind barbs (None, for 9)
+      windn= name of the wind variable in the graph
+      wuts= units of the wind variable in the graph
+      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
+      graphtit= title of the graph ('|', for spaces)
+      kfig= kind of figure
+      figname= name of the figure
+    """
+    fname = 'plot_barbs'
+ 
+    dx=xvals.shape[1]
+    dy=xvals.shape[0]
+
+# Frequency of vectors
+    if vecfreq is None:
+        xfreq = 1
+        yfreq = 1
+    elif vecfreq == 'auto':
+        xfreq = dx/20
+        yfreq = dy/20
+    else:
+        xfreq=int(vecfreq.split('@')[0])
+        yfreq=int(vecfreq.split('@')[1])
+
+    if veccolor is None:
+        vcolor = "red"
+    else:
+        vcolor = veccolor
+
+    if veclength is None:
+        vlength = 9
+    else:
+        vlength = veclength
+
+    plt.rc('text', usetex=True)
+
+    if not mapv is None:
+        lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:])
+        lat00 = yvals[:]
+
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        nlon = np.min(xvals[::yfreq,::xfreq])
+        xlon = np.max(xvals[::yfreq,::xfreq])
+        nlat = np.min(yvals[::yfreq,::xfreq])
+        xlat = np.max(yvals[::yfreq,::xfreq])
+
+        lon2 = xvals[dy/2,dx/2]
+        lat2 = yvals[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 + ": projection '" + map_proj + "' not ready!!"
+            print '    projections available: cyl, lcc'
+            quit(-1)
+
+        m.drawcoastlines()
+
+        meridians = pretty_int(nlon,xlon,5)
+        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")
+
+        parallels = pretty_int(nlat,xlat,5)
+        m.drawparallels(parallels,labels=[False,True,True,False],color="black")
+
+        plt.xlabel('W-E')
+        plt.ylabel('S-N')
+
+    plt.barbs(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq],                        \
+      uvals[timev,lev,::yfreq,::xfreq],vvals[timev,lev,::yfreq,::xfreq],color=vcolor,\
+      pivot='tip')
+
+    plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',           \
+      xy=(0.85,-0.10), xycoords='axes fraction', color=vcolor)
+
+    plt.title(graphtit.replace('|',' ').replace('&','\&'))
+
+## NOT WORKING ##
+
+# No legend so it is imposed
+##    windlabel=windn.replace('_','\_') +' (' + units_lunits(wuts[1]) + ')'
+##    vecpatch = mpatches.Patch(color=vcolor, label=windlabel)
+
+##    plt.legend(handles=[vecpatch])
+
+##    vecline = mlines.Line2D([], [], color=vcolor, marker='.', markersize=10, label=windlabel)
+##    plt.legend(handles=[vecline], loc=1)
+
+    output_kind(kfig, figname, True)
+
+    return