Changeset 1930 in lmdz_wrf for trunk

Timestamp:

Jul 13, 2018, 11:17:30 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `draw_SkewT': creation of a SkewT-logP diagram using matplotlib's API example

Location:

trunk/tools

Files:

• drawing.py (modified) (5 diffs)
• drawing_tools.py (modified) (2 diffs)
• skewt.py (added)

trunk/tools/drawing.py

@@ -57 +57 @@
 ## e.g. # drawing.py -o draw_multi_2D_shad -f '../PY/wrfout_d01_1995-01-01_00:00:00@T2@west_east|-1,south_north|-1,Time|0;../PY/wrfout_d01_1995-01-01_00:00:00@T2@west_east|-1,south_north|-1,Time|1;../PY/wrfout_d01_1995-01-01_00:00:00@T2@west_east|-1,south_north|-1,Time|2;../PY/wrfout_d01_1995-01-01_00:00:00@T2@west_east|-1,south_north|-1,Time|3' -S 'tas:XLONG:XLAT:auto:rainbow,auto,auto:Srange,Srange:png:None:cyl,l:0!UTC,1!UTC,2!UTC,3!UTC:2:2:tas!at!2001-11-11:True'
 # movie_2D_shad -f '../PY/wrfout_d01_1995-01-01_00:00:00' -S 'tas:west_east|-1,south_north|-1,Time|-1:XLONG:XLAT:auto:rainbow,auto,auto:Srange,Srange:png:None:cyl,l:Time,Times:WRFdate,$%d^{%H}$:15:mp4' -v T2
+## e.g. # drawing.py -o draw_SkewT -f UWyoming_snd_87576.nc -S 'time|0:auto:auto:Sounding!at!Ezeiza!airport!on!3rd!July!1998:png:yes' -v ta,tda,pres
 
 #######
@@ -85 +86 @@
 # draw_points_lonlat: Function to plot a series of lon/lat points
 # draw_ptZvals: Function to plot a given list of points by their Z value and a colorbar
+# draw_SkewT: creation of a SkewT-logP diagram using matplotlib's API example  
 # draw_time_lag: Function to plot a time-lag figure with multiple sources (x, previous values; y, future values)
 # draw_timeSeries: Function to draw a time-series
@@ -116 +118 @@
   'draw_Neighbourghood_evol',                                                        \
   'draw_points', 'draw_points_lonlat',                                               \
-  'draw_ptZvals', 'draw_river_desc', 'draw_subbasin', 'draw_Taylor',                 \
+  'draw_ptZvals', 'draw_river_desc', 'draw_SkewT', 'draw_subbasin', 'draw_Taylor',   \
   'draw_time_lag', 'draw_timeSeries', 'draw_topo_geogrid',                           \
   'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
@@ -8871 +8873 @@
 #python drawing.py -o movie_2D_shad -f '../PY/wrfout_d01_1995-01-01_00:00:00' -S 'tas:west_east|-1,south_north|-1,Time|-1:XLONG:XLAT:auto:rainbow,auto,auto:Srange,Srange:png:None:cyl,l:Time,Times:WRFdate,$%d^{%H}$:15:mp4' -v T2
 
+def draw_SkewT(ncfile, values, varname):
+    """ creation of a SkewT-logP diagram using matplotlib's API example  
+      https://matplotlib.org/examples/api/skewt.html (skewt.py, included as external python)
+    draw_SkewT(ncfile, values, varn)
+      ncfile= file to use
+      values = [dimvals]:[taminv],[tamaxv]:[presminv],[presmaxv]:[figt]:[kindfig]:
+        [close]
+        [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
+        [taminv],[tamaxv]: minimum and maximum of temperatures to plot [C] 
+          ('auto' for -50,50)
+        [presminv],[presmaxv]: minimum and maximum of pressures to plot [hPa]
+          ('auto' for 100.,1020.)
+        [figt]: title of figure ('!' for spaces)
+        [kindfig]: kind of graphical output (ps, png, pdf)
+        [close]: wether figure should be closed or not
+      varn= [tavarn],[tdavarn],[presvarn]
+    """
+    fname = 'draw_SkewT'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_SkewT.__doc__
+        quit()
+
+    expectargs = '[dimvals]:[taminv],[tamaxv]:[presminv],[presmaxv]:[figt]:' +       \
+      '[kindfig]:[close]'
+ 
+    drw.check_arguments(fname,values,expectargs,':')
+
+    # Variables
+    varns = varname.split(',')
+    if len(varns) != 3:
+        print errormsg
+        print '  ' +fname+ ": requires 3 name of variables as '[tavarn],[tdavarn]" + \
+          ",[tavarn]' !!"
+        print "    variables provided: '" + varname + "' "
+        quit(-1)
+
+    dimvals= values.split(':')[0].replace('|',':')
+    if values.split(':')[1] != 'auto':
+        tanx = gen.str_list_k(values.split(':')[1], ',', 'np.float')
+    else:
+        tanx = [-50.,50.]
+    if values.split(':')[2] != 'auto':
+        presnx = gen.str_list_k(values.split(':')[2], ',', 'np.float')
+    else:
+        presnx = [100.,1020.]
+    figt = values.split(':')[3].replace('!',' ')
+    kindfig = values.split(':')[4]
+    close = gen.Str_Bool(values.split(':')[5])
+
+    onc = NetCDFFile(ncfile, 'r')
+    for varn in varns:
+        if not onc.variables.has_key(varn):
+            print errormsg
+            print ' ' + fname + ": file '" + ncfile + "' does not have variable '" + \
+              varn + "' !!"
+            print '    available ones:', onc.variables-keys()
+            quit(-1)
+
+    ota = onc.variables[varns[0]]
+    otda = onc.variables[varns[1]]
+    opres = onc.variables[varns[2]]
+
+    # Getting values
+    tavals, tadims = drw.slice_variable(ota, dimvals.replace(',','|'))
+    tdavals, tdadims = drw.slice_variable(otda, dimvals.replace(',','|'))
+    presvals, presdims = drw.slice_variable(opres, dimvals.replace(',','|'))
+
+    # Expected masked values
+    if type(tavals) != type(gen.mamat):
+        tavals = ma.masked_array(tavals)
+    if type(tdavals) != type(gen.mamat):
+        tdavals = ma.masked_array(tdavals)
+    if type(tavals) != type(gen.mamat):
+        presvals = ma.masked_array(presvals)
+
+    drw.plot_SkewT(tavals, tdavals, presvals, tanx, presnx, figt, kindfig, close)
+
+    onc.close()
+
+    return
+
+#filen='/home/lluis/estudios/ChemGBsAs/tests/199807/obs/snd/UWyoming_snd_87576.nc'
+#values='time|0:auto:auto:Sounding!at!Ezeiza!airport!on!3rd!July!1998:png:yes'
+#varn='ta,tda,pres'
+#draw_SkewT(filen, values, varn)
+
 #quit()
 
@@ -8994 +9090 @@
     elif oper == 'draw_subbasin':
         draw_subbasin(opts.ncfile, opts.values)
+    elif oper == 'draw_SkewT':
+        draw_SkewT(opts.ncfile, opts.values, opts.varname)
     elif oper == 'draw_Taylor':
         draw_Taylor(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -143 +143 @@
 # plot_cycle: Function to plot a variale with a circular cycle
 # plot_lines: Function to plot a collection of lines
+# plot_SkewT: Function to plot a SkewT-logP diagram using matplotlib's API example 
 # plot_Taylor: Function to draw a Taylor diagram (Taylor 2001)
 # plot_TimeEnsembles: Function to plot a of an Ensemble of values with a time-axis
@@ -11369 +11370 @@
     return
 
-
+def SaturationPressure(temp):
+    """ Provides de saturation pressure at a given temperature following The August-Roche-Magnus equation
+      temp: temperature [K]
+    """
+    fname = 'SaturationPressure'
+
+    ARM1 = 6.1094 # hPa
+    ARM2 = 17.625 # C
+    ARM3 = 243.04 # C
+    
+    tempC = temp - 273.15  
+    satpres = 100.*ARM1*(ARM2*tempC)/(tempC+ARM3)
+
+    return satpres
+
+#print 'satpres:', SaturationPressure(293.15)
+
+def SatMixingRatio(ta,pres):
+    """ Provides de saturation mixing ratio
+      ta: temperature [K]
+      pres: pressure [Pa]
+    """
+    fname = 'SatMixingRatio'
+
+    es = SaturationPressure(ta)
+    rs = (0.6222*es)/(pres-es)
+
+    return rs
+
+#print 'saturation mixing ratio:', SatMixingRatio(293.,101300.)
+
+def ethetapotct(tsfc, dz=50, pmin=5000.):
+    """ Provides the curve for a given constant equivalent potential temperature starting at surface
+     After: 
+      http://www.atmo.arizona.edu/students/courselinks/fall14/atmo551a/Site/ATMO_451a_551a_files/WaterVapor2.pdf
+      https://www.meteo.physik.uni-muenchen.de/~roger/Lectures/Moist%20Convection/Conv_L03_2005.pdf
+      
+      tsfc= temperature at the surface (C)
+      dz= number of points for the curve
+      R =  287.058 (Jkg-1k-1) # 8.3144598 (Jmol-1K-1)
+      Cp =  1.0035 (Jg-1K-1) # 29.07 (Jmol-1K-1)
+      p0 = 100000 # (Pa)
+      psfc = 101325 # (Pa)
+      L = 2264.705 (kJkg-1)
+    """
+    fname = 'ethetapotct'
+    R = 287.058 # (Jkg1-K-1)
+    Cp = 1003.5 # (Jkg1-K-1)
+    psfc = 101325. # Pa
+    p0 = 100000. # Pa
+    L = 2264.705 # kJkg-1
+
+    ethetapotv = np.zeros((dz,2), dtype=np.float)
+
+    dp = (psfc-pmin)/(dz-1)
+    
+    for iz in range(dz):
+        ethetapotv[iz,0] = psfc-dp*iz
+        thetapotv = (tsfc+273.15)*((psfc-dp*iz)/p0)**(R/Cp) 
+#        print iz,psfc-dp*iz,thetapotv,SatMixingRatio(tsfc+273.15,psfc),np.exp(L*1000.*SatMixingRatio(tsfc+273.15,psfc)/Cp*(tsfc+273.15))
+        ethetapotv[iz,1] = thetapotv*np.exp(L*1000.*SatMixingRatio(tsfc+273.15,psfc)/(Cp*(tsfc+273.15))) - 273.15
+
+#    quit()
+    return ethetapotv
+
+
+def thetapotct(tsfc, dz=50, pmin=5000.):
+    """ Provides the curve for a given constant potential temperature starting at surface
+      tsfc= temperature at the surface (C)
+      dz= number of points for the curve
+      R =  287.058 (Jkg-1k-1) # 8.3144598 (Jmol-1K-1)
+      Cp =  1.0035 (Jg-1K-1) # 29.07 (Jmol-1K-1)
+      p0 = 100000 # (Pa)
+      psfc = 101325 # (Pa)
+    """
+    fname = 'thetapotct'
+    R = 287.058 # (Jkg1-K-1)
+    Cp = 1003.5 # (Jkg1-K-1)
+    psfc = 101325. # Pa
+    p0 = 100000. # Pa
+
+    thetapotv = np.zeros((dz,2), dtype=np.float)
+    dp = (psfc-pmin)/(dz-1)
+    
+    for iz in range(dz):
+        thetapotv[iz,0] = psfc-dp*iz
+        thetapotv[iz,1] = (tsfc+273.15)*((psfc-dp*iz)/p0)**(R/Cp) - 273.15
+
+    return thetapotv
+
+def plot_SkewT(ta, tda, pres, taxtrm, presxtrm, figtitle, kfig, close):
+    """ Function to plot a SkewT-logP diagram using matplotlib's API example  
+      ta: air temperature [K]
+      tda: dew point air temperature [K]
+      pres: pessures for ta and tda [Pa]
+      figtitle: title of figure
+      kfig: kind of graphical output (jpg, pdf, png)
+      close: whether figure should be closed or not
+    """
+    import skewt
+
+    fname = 'plot_SkewT'
+    # Checking for variable consistency
+    if len(ta.shape) != 1:
+        print errormsg
+        print '  ' + fname + ": values of 'ta' must have 1D-rank!!"
+        print '    passed shape:', ta.shape
+        quit(-1)
+    if len(tda.shape) != 1:
+        print errormsg
+        print '  ' + fname + ": values of 'tda' must have 1D-rank!!"
+        print '    passed shape:', tda.shape
+        quit(-1)
+    if len(pres.shape) != 1:
+        print errormsg
+        print '  ' + fname + ": values of 'pres' must have 1D-rank!!"
+        print '    passed shape:', pres.shape
+        quit(-1)
+
+    # Create a new figure. The dimensions here give a good aspect ratio
+    fig = plt.figure(figsize=(6.5875, 6.2125))
+    ax = fig.add_subplot(111, projection='skewx')
+
+    plt.grid(True)
+
+    # Quantity of dryadiabats
+    xmax = taxtrm[1]
+    xmin = taxtrm[0]
+    print 'len:', np.arange(xmin,xmax,10.)
+    if len(np.arange(xmin,xmax,10.)) < 10: xfrqtk = 10.
+    else: xfrqtk = 20.
+
+    # Including dry adiabats
+    ddt = int((xmax+xfrqtk-xmin)/xfrqtk)
+    for itt in range(ddt):
+        dryt = thetapotct(xmin+itt*xfrqtk, pmin=10000.)
+        ax.semilogy(dryt[:,1], dryt[:,0]/100., color='#FFCCCC', linewidth=0.75)
+
+    # Including moist- adiabats
+    #for itt in range(ddt):
+    #    moist = ethetapotct(xmin+itt*xfrqtk, pmin=10000.)
+    #    if xmin+itt*xfrqtk == 10.: print moist[:,1], moist[:,0]/100.
+    #    ax.semilogy(moist[:,1], moist[:,0]/100., color='#CCFFCC', linewidth=0.75)   
+
+    # Plot the data using normal plotting functions, in this case using
+    # log scaling in Y, as dictated by the typical meteorological plot
+    tamask = ta.mask
+    tdamask = tda.mask
+
+    ax.semilogy(ta.compressed(), pres[~tamask], color='C3')
+    ax.semilogy(tda.compressed(), pres[~tdamask], color='C2')
+
+    # An example of a slanted line at constant X
+    l = ax.axvline(0, color='C0')
+
+    # Disables the log-formatting that comes with semilogy
+    ax.yaxis.set_major_formatter(ScalarFormatter())
+    ax.yaxis.set_minor_formatter(NullFormatter())
+    ax.set_yticks(np.linspace(100, 1000, 10))
+    ax.set_ylim(presxtrm[1], presxtrm[0])
+
+    ax.xaxis.set_major_locator(MultipleLocator(10))
+    ax.set_xlim(taxtrm[0], taxtrm[1])
+
+    # Indices
+    #plt.annotate('CAPE: ' + str(stvals['CAPE']), xy=(0.98,0.97),                         \
+    #  xycoords='axes fraction', ha="right", backgroundcolor='#DDDDDD')
+    #plt.annotate('CI: '+ str(stvals['CI']), xy=(0.98,0.93), xycoords='axes fraction',    \
+    #  ha='right', backgroundcolor='#DDDDDD')
+
+    ax.set_title(gen.latex_text(figtitle))
+
+    figname = 'SkewT'
+    output_kind(kfig, figname, close)
+
+    return