Changeset 1947 in lmdz_wrf for trunk


Ignore:
Timestamp:
Jul 19, 2018, 4:58:53 AM (6 years ago)
Author:
lfita
Message:

Adding:

  • `draw_2Dshad_contdisc': 2D shadding field with a discrete one
  • `draw_2Dshad_contdisc_time': 2D shadding field with a discrete one wtih a time-axis
  • Fixing 'interpolation'
  • Adding more variables
Location:
trunk/tools
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/drawing.py

    r1946 r1947  
    6060## 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
    6161## e.g. # drawing.py -o draw_multi_SkewT -f '/home/lluis/estudios/ChemGBsAs/tests/199807/obs/snd/UWyoming_snd_87576.nc:pres|-1,time|0:ta,pres;/home/lluis/estudios/ChemGBsAs/tests/199807/obs/snd/UWyoming_snd_87576.nc:pres|-1,time|0:tda,pres' -S 'auto:auto:multilines!ta,tda!#6464FF,#FF6464!-!,!2:0,auto:Sounding!at!Ezeiza!airport!on!3rd!July!1998:png:yes'
     62## e.g. # drawing.py -o draw_2D_shad_contdisc -f '/home/lluis/estudios/ChemGBsAs/tests/199501/simout_snddiags.nc;ta;time;pres;time|-1,pres|-1@/home/lluis/estudios/ChemGBsAs/tests/199501/obs/snd/UWyoming_snd_87576.nc;ta;time;pres;time|-1,pres|-1' -S 'ta:time,bottom_top:Vfix,auto,3600.,auto,Vfix,auto,50.,auto:auto:Srange,Srange:auto:obs!&!sim!Ezeiza!airport!sounding:pdf:flip@y:None:yes'
    6263
    6364#######
    … …  
    6566# draw_2D_shad: plotting a fields with shading
    6667# draw_2D_shad_cont: plotting two fields, one with shading and the other with contour lines
     68# draw_2D_shad_contdisc: plotting one continuous fields with shading and another discrete one with points
    6769# draw_2D_shad_2cont: plotting three fields, one with shading and the other two with contour lines
    6870# draw_2D_shad_cont_time: plotting two fields, one with shading and the other with contour lines being
    … …  
    113115
    114116namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time',                 \
    115   'draw_2D_shad_cont', 'draw_2D_shad_2cont', 'draw_2D_shad_cont_time',               \
     117  'draw_2D_shad_cont', 'draw_2D_shad_contdisc', 'draw_2D_shad_2cont',                \
     118  'draw_2D_shad_cont_time',                                                          \
    116119  'draw_2D_shad_line',                                                               \
    117120  'draw_2D_shad_line_time', 'draw_bar', 'draw_bar_line', 'draw_bar_line_time',       \
    … …  
    91269129#draw_multi_SkewT(filen, values)
    91279130
     9131def draw_2D_shad_contdisc(ncfiles, values, axfig=None, fig=None):
     9132    """ plotting one continuous fields with shading and another discrete one with points
     9133    draw_2D_shad_contdisc(ncfile, values)
     9134      ncfiles= [contfilen];[contvarn];[dimxvarn];[dimyvarn];[dimvals]@[discfilen];
     9135          [discvarn];[dimxvarn];[dimyvarn];[dimvals] files and variables to use
     9136        [contfilen]: name of the file with the continuous varible
     9137        [contvarn]: name of the continuos variable
     9138        [dimxvarn]: name of the variable with values for the x-dimension
     9139        [dimyvarn]: name of the variable with values for the y-dimension
     9140        [dimvals]: ',' list of [dimname]|[value] values to slice the variable:
     9141            * [integer]: which value of the dimension
     9142            * -1: all along the dimension
     9143            * -9: last value of the dimension
     9144            * [beg]@[end]@[inc] slice from [beg] to [end] every [inc]
     9145            * NOTE, no dim name all the dimension size
     9146        [discfilen]: name of the file with the discrete varible
     9147        [discvarn]: name of the discrete variable
     9148        * NOTE: limits of the graph will be computed from the continuous variable
     9149      values=[vnamefs]:[dvarxn],[dvaryn]:[dimxyfmt]:[colorbarvals]:[sminv],[smaxv]:
     9150       [discvals]:[figt]:[kindfig]:[reverse]:[mapv]:[close]
     9151        [vnamefs]: Name in the figure of the variable to be shaded
     9152        [dvarxn],[dvaryn]: name of the dimensions for the final x-axis and y-axis at
     9153          the figure (from contfilen)
     9154        [dimxyfmt]=[dxs],[dxf],[Ndx],[ordx],[dys],[dyf],[Ndy],[ordx]: format of the
     9155            values at each axis (or 'auto')
     9156          [dxs]: style of x-axis ('auto' for 'pretty')
     9157            'Nfix', values computed at even 'Ndx'
     9158            'Vfix', values computed at even 'Ndx' increments
     9159            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
     9160          [dxf]: format of the labels at the x-axis ('auto' for '%5g')
     9161          [Ndx]: Number of ticks at the x-axis ('auto' for 5)
     9162          [ordx]: angle of orientation of ticks at the x-axis ('auto' for horizontal)
     9163          [dys]: style of y-axis ('auto' for 'pretty')
     9164          [dyf]: format of the labels at the y-axis ('auto' for '%5g')
     9165          [Ndy]: Number of ticks at the y-axis ('auto' for 5)
     9166          [ordy]: angle of orientation of ticks at the y-axis ('auto' for horizontal)
     9167        [colorbarvals]=[colbarn],[fmtcolorbar],[orientation]
     9168          [colorbarn]: name of the color bar
     9169          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
     9170          [orientation]: orientation of the colorbar ('vertical' ['auto'], 'horizontal')
     9171          * NOTE: single 'auto' for 'rainbow,%6g,vertical'
     9172        [smin/axv]: minimum and maximum values for the shading or string for each for:
     9173          'Srange': for full range
     9174          'Saroundmean@val': for mean-xtrm,mean+xtrm where
     9175            xtrm = np.min(mean-min@val,max@val-mean)
     9176          'Saroundminmax@val': for min*val,max*val
     9177          'Saroundpercentile@val': for median-xtrm,median+xtrm where
     9178            xtrm = np.min(median-percentile_(val),percentile_(100-val)-median)
     9179          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
     9180          'Smedian@val': for -xtrm,xtrm where
     9181             xtrm = np.min(median-min@val,max@val-median)
     9182          'Spercentile@val': for -xtrm,xtrm where
     9183             xtrm = np.min(median-percentile_(val),percentile_(100-val)-median)
     9184        [discvals]= [type],[size],[lwidth],[lcol] characteristics of the points for the discrete field
     9185          [type]: type of point. Any marker from matoplib must be filled !
     9186          [size]: size of point
     9187          [lwidth]: width of the line around the point
     9188          [lcol]: color of the line around the point
     9189          'auto': for [type]='o', [size]=5, [lwdith]=0.25, [lcol]='#000000'
     9190        [figt]: title of the figure ('!' for spaces)
     9191        [kindfig]: kind of figure output (ps, png, pdf)
     9192        [reverse]: Transformation of the values
     9193          * 'transpose': reverse the axes (x-->y, y-->x)
     9194          * 'flip'@[x/y]: flip the axis x or y
     9195        [mapv]: map characteristics: [proj],[res]
     9196          see full documentation: http://matplotlib.org/basemap/
     9197          [proj]: projection
     9198            * 'cyl', cilindric
     9199            * 'lcc', lambert conformal
     9200          [res]: resolution:
     9201            * 'c', crude
     9202            * 'l', low
     9203            * 'i', intermediate
     9204            * 'h', high
     9205            * 'f', full
     9206        [close]: Whether figure should be finished or not
     9207    """
     9208
     9209    fname = 'draw_2D_shad_contdisc'
     9210    if values == 'h':
     9211        print fname + '_____________________________________________________________'
     9212        print draw_2D_shad_contdisc.__doc__
     9213        quit()
     9214
     9215    expectargs = '[vnamefs]:[dvarxn],[dvaryn]:[dimxyfmt]:[colorbarvals]:[sminv],' +  \
     9216      '[smaxv]:[discvals]:[figt]:[kindfig]:[reverse]:[mapv]:[close]'
     9217 
     9218    drw.check_arguments(fname,values,expectargs,':')
     9219
     9220    vnamesfig = values.split(':')[0]
     9221    [dvarxn, dvaryn] = values.split(':')[1].split(',')
     9222    dimxyf = values.split(':')[2]
     9223    colorbarvals = values.split(':')[3]
     9224    shadminmax = values.split(':')[4]
     9225    discvals = values.split(':')[5]
     9226    figt = values.split(':')[6].replace('!',' ')
     9227    kindfig = values.split(':')[7]
     9228    revals = values.split(':')[8]
     9229    mapvalue = values.split(':')[9]
     9230    close = gen.Str_Bool(values.split(':')[10])
     9231
     9232    ncs = ncfiles.split('@')
     9233    filesinf = {}
     9234    ifile = 0
     9235    for nc in ncs:
     9236        if ifile == 0:
     9237            hfn = 'cont'
     9238        else:
     9239            hfn = 'disc'
     9240
     9241        filesinf[hfn + 'filen'] = nc.split(';')[0]
     9242        filesinf[hfn + 'varn'] = nc.split(';')[1]
     9243        filesinf[hfn + 'dimxvn'] = nc.split(';')[2]
     9244        filesinf[hfn + 'dimyvn'] = nc.split(';')[3]
     9245        filesinf[hfn + 'dimvals'] = nc.split(';')[4].replace('|',':')
     9246        ifile = ifile + 1
     9247
     9248    for hfn in ['cont', 'disc']:
     9249        ncfilen = filesinf[hfn + 'filen']
     9250        varn = filesinf[hfn + 'varn']
     9251        vdimxn = filesinf[hfn + 'dimxvn']
     9252        vdimyn = filesinf[hfn + 'dimyvn']
     9253        dimvals = filesinf[hfn + 'dimvals']
     9254
     9255        if not os.path.isfile(ncfilen):
     9256            print errormsg
     9257            print '  ' + fname + ": '" + hfn + "' file '" + ncfile + "' does " +     \
     9258              "not exist !!"
     9259            quit(-1)
     9260
     9261        onc = NetCDFFile(ncfilen, 'r')
     9262
     9263        if  not onc.variables.has_key(varn):
     9264            print errormsg
     9265            print '  ' + fname + ": '" + hfn + "' file '" + ncfile +                 \
     9266              "' does not have variable '" + varn + "' !!"
     9267            varns = sorted(onc.variables.keys())
     9268            print '    available ones:', varns
     9269            quit(-1)
     9270
     9271        # Variables' values
     9272        ovar = onc.variables[varn]
     9273        vals, dims = drw.slice_variable(ovar, dimvals.replace(',','|'))
     9274
     9275        # Values have to be 2D (might not be true for disc?)
     9276        if len(vals.shape) != 2:
     9277            print errormsg
     9278            print '  ' + fname + ": values have to be 2D !!"
     9279            print '    provided shape:', vals.shape, 'for slice:', dimvals
     9280            quit(-1)
     9281
     9282        if drw.searchInlist(ovar.ncattrs(),'units'):
     9283            varunits = ovar.getncattr('units')
     9284        else:
     9285            print warnmsg
     9286            print '  ' + fname + ": variable '" + varn + "' without units!!"
     9287            varunits = '-'
     9288
     9289        # dimensions
     9290        dimnamesv = [vdimxn, vdimyn]
     9291        if  not onc.variables.has_key(vdimxn):
     9292            print errormsg
     9293            print '  ' + fname + ": '" + hfn + "' file '" + ncfile +                 \
     9294              "' does not have x-dimension variable '" +  vdimxn + "' !!"
     9295            varns = sorted(onc.variables.keys())
     9296            print '    available ones:', varns
     9297            quit(-1)
     9298        if  not onc.variables.has_key(vdimyn):
     9299            print errormsg
     9300            print '  ' + fname + ": '" + hfn + "' file '" + ncfile +                 \
     9301              "' does not have y-dimension variable '" +  vdimyn + "' !!"
     9302            varns = sorted(onc.variables.keys())
     9303            print '    available ones:', varns
     9304            quit(-1)
     9305
     9306        objdimx = onc.variables[vdimxn]
     9307        objdimy = onc.variables[vdimyn]
     9308        if drw.searchInlist(objdimx.ncattrs(),'units'):
     9309            odimxu = objdimx.getncattr('units')
     9310        else:
     9311            print warnmsg
     9312            print '  ' +fname+": variable dimension '" + vdimxn + "' without units!!"
     9313            odimxu = '-'
     9314
     9315        if drw.searchInlist(objdimy.ncattrs(),'units'):
     9316            odimyu = objdimy.getncattr('units')
     9317        else:
     9318            print warnmsg
     9319            print '  ' +fname+": variable dimension '" + vdimyn + "' without units!!"
     9320            odimyu = '-'
     9321
     9322        odimxv0, odimyv0 = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:],                \
     9323          objdimx.dimensions, objdimy.dimensions, dimvals.replace(':','|').split(','),\
     9324          False)
     9325
     9326        # Some statistics
     9327        vn = np.min(vals)
     9328        vx = np.max(vals)
     9329        dxn = np.min(odimxv0)
     9330        dxx = np.max(odimxv0)
     9331        dyn = np.min(odimyv0)
     9332        dyx = np.max(odimyv0)
     9333
     9334        print gen.infmsg
     9335        print '  ' + fname + ": for '" + hfn + "' ________"
     9336        print '    values min:', vn, 'max:', vx
     9337        print '    dimx min:', dxn, 'max:', dxx
     9338        print '    dimy min:', dyn, 'max:', dyx
     9339
     9340        if hfn == 'cont':
     9341            if not gen.searchInlist(onc.dimensions, dvarxn):
     9342                print errormsg
     9343                print '  ' + fname + ": file '" + ncfilen + "' does not have " +     \
     9344                  "dimension '" + dvarxn + "' for final x-axis in figure !!"
     9345                print '     available ones: ', onc.dimensions.keys()
     9346                quit(-1)
     9347            if not gen.searchInlist(onc.dimensions, dvaryn):
     9348                print errormsg
     9349                print '  ' + fname + ": file '" + ncfilen + "' does not have " +     \
     9350                  "dimension '" + dvaryn + "' for final y-axis in figure !!"
     9351                print '     available ones: ', onc.dimensions.keys()
     9352                quit(-1)
     9353
     9354            # Getting the right shape of dimensions
     9355            gdimx = len(onc.dimensions[dvarxn])
     9356            gdimy = len(onc.dimensions[dvaryn])
     9357
     9358            if vals.shape[1] == gdimx: contv = vals[:]
     9359            else: contv = vals[:].transpose()
     9360
     9361            odxv = np.zeros((gdimy,gdimx), dtype=np.float)
     9362            odyv = np.zeros((gdimy,gdimx), dtype=np.float)
     9363
     9364            if len(odimxv0.shape) == 2:
     9365                if odimxv0.shape[1] == gdimx: odxv = odimxv0[:]
     9366                else: odxv = odimxv0.transpose()
     9367            else:
     9368                for k in range(gdimy):
     9369                    odxv[k,:] = odimxv0[:]
     9370
     9371            if len(odimyv0.shape) == 2:
     9372                if odimyv0.shape[0] == gdimy: odyv = odimyv0[:]
     9373                else: odyv = odimyv0.transpose()
     9374            else:
     9375                for k in range(gdimx):
     9376                    odyv[:,k] = odimyv0[:]
     9377
     9378            diminfo = {}
     9379            diminfo['units'] = [odimxu, odimyu]
     9380            diminfo['names'] = [gen.variables_values(vdimxn)[0],                     \
     9381              gen.variables_values(vdimyn)[1]]
     9382            varinfo = {}
     9383            varinfo['units'] = varunits
     9384            varinfo['name'] = vnamesfig
     9385
     9386            # Absolute xtremes for the plot
     9387            absxn = dxn
     9388            absxx = dxx
     9389            absyn = dyn
     9390            absyx = dyx
     9391
     9392        else:
     9393            discvarv = []
     9394            discx = []
     9395            discy = []
     9396            # Assuming that discrete values are masked
     9397            if type(vals) != type(gen.mamat):
     9398                madiscvarv = ma.masked_array(vals)
     9399            else:
     9400                madiscvarv = vals
     9401
     9402            # Getting 2D dimensions
     9403            [d2Dx, d2Dy] = np.meshgrid(odimxv0, odimyv0)
     9404
     9405            # Checking for flipping axis?
     9406            if d2Dx.shape[0] == vals.shape[1] and d2Dx.shape[1] == vals.shape[0]:
     9407                print '  rotating values !'
     9408                #madiscvarv = madiscvarv.transpose()
     9409                d2Dx = d2Dx.transpose()
     9410                d2Dy = d2Dy.transpose()
     9411
     9412            # vals have to be 2D otherwise... ?
     9413            dx = madiscvarv.shape[1]
     9414            dy = madiscvarv.shape[0]
     9415            for j in range(dy):
     9416                for i in range(dx):
     9417                    if not madiscvarv.mask[j,i]:
     9418                        discvarv.append(madiscvarv[j,i])
     9419                        discx.append(d2Dx[j,i])
     9420                        discy.append(d2Dy[j,i])
     9421
     9422            Ndiscvals = len(discvarv)
     9423            print '    getting:', Ndiscvals, 'discrete values'
     9424            discv = np.zeros((Ndiscvals,3), dtype=np.float)
     9425            discv[:,0] = discx[:]
     9426            discv[:,1] = discy[:]
     9427            discv[:,2] = discvarv[:]
     9428
     9429            # Absolute xtremes for the plot
     9430            absxn = np.max([absxn,dxn])
     9431            absxx = np.min([absxx,dxx])
     9432            absyn = np.max([absyn,dyn])
     9433            absyx = np.min([absyx,dyx])
     9434
     9435        onc.close()
     9436
     9437    graphnx = [absxn, absxx, absyn, absyx]
     9438    print 'limits of the graphic:', graphnx
     9439
     9440    shading_nx = []
     9441    if shadminmax.split(',')[0][0:1] != 'S':
     9442            shading_nx.append(np.float(shadminmax.split(',')[0]))
     9443    else:
     9444        shading_nx.append(shadminmax.split(',')[0])
     9445
     9446    if shadminmax.split(',')[1][0:1] != 'S':
     9447        shading_nx.append(np.float(shadminmax.split(',')[1]))
     9448    else:
     9449        shading_nx.append(shadminmax.split(',')[1])
     9450
     9451    if discvals == 'auto':
     9452        discinfo = ['o', 5., 0.25, '#000000']
     9453    else:
     9454        discinfo = [discvals.split(',')[0], np.float(discvals.split(',')[1]),        \
     9455          np.float(discvals.split(',')[2]), discvals.split(',')[3]]
     9456
     9457    if mapvalue == 'None': mapvalue = None
     9458
     9459    if colorbarvals == 'auto': colorbarvals = 'rainbow,auto,auto'
     9460    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvals,',')
     9461    colormapv = [colbarn, fmtcolbar, colbaror]
     9462
     9463    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimxyf,',')
     9464    xaxis = [xstyl, xaxf, Nxax, xaxor]
     9465    yaxis = [ystyl, yaxf, Nyax, yaxor]
     9466
     9467    if revals == 'None':
     9468        revals = None
     9469
     9470    drw.plot_2Dshad_obssim(discv, contv, odxv, odyv, revals, diminfo, xaxis, yaxis,  \
     9471      colormapv, shading_nx, varinfo, discinfo, mapvalue, graphnx, figt, kindfig,    \
     9472      close)
     9473
     9474    return
     9475
     9476#contfilen = '/home/lluis/estudios/ChemGBsAs/tests/199501/simout_snddiags.nc;ta;time;pres;time|-1,pres|-1'
     9477#discfilen = '/home/lluis/estudios/ChemGBsAs/tests/199501/obs/snd/UWyoming_snd_87576.nc;ta;time;pres;time|-1,pres|-1'
     9478
     9479#values = 'ta:time,bottom_top:Vfix,auto,3600.,auto,Vfix,auto,50.,auto:auto:Srange,Srange:o,5.:obs!&!sim!Ezeiza!airport!sounding:pdf:flip@y:None:yes'
     9480
     9481#draw_2D_shad_contdisc(contfilen+'@'+discfilen, values, axfig=None, fig=None)
     9482   
     9483def draw_2D_shad_contdisc_time(ncfiles, values, axfig=None, fig=None):
     9484    """ plotting one continuous fields with shading and another discrete one with points
     9485    draw_2D_shad_contdisc(ncfile, values)
     9486      ncfiles= [contfilen];[contvarn];[dimtvarn];[dimvvarn];[dimvals]@[discfilen];
     9487          [discvarn];[dimtvarn];[dimvvarn];[dimvals] files and variables to use
     9488        [contfilen]: name of the file with the continuous varible
     9489        [contvarn]: name of the continuos variable
     9490        [dimtvarn]: name of the variable with values for the time-dimension
     9491        [dimvvarn]: name of the variable with values for the variable-dimension
     9492        [dimvals]: ',' list of [dimname]|[value] values to slice the variable:
     9493            * [integer]: which value of the dimension
     9494            * -1: all along the dimension
     9495            * -9: last value of the dimension
     9496            * [beg]@[end]@[inc] slice from [beg] to [end] every [inc]
     9497            * NOTE, no dim name all the dimension size
     9498        [discfilen]: name of the file with the discrete varible
     9499        [discvarn]: name of the discrete variable
     9500        * NOTE: limits of the graph will be computed from the continuous variable
     9501      values=[vnamefs];[varaxis];[timevals];[dimxyfmt];[colorbarvals];[sminv],[smaxv];
     9502       [discvals];[figt];[kindfig];[reverse];[mapv];[close]
     9503        [vnamefs]: Name in the figure of the variable to be shaded
     9504        [varaxis]: axis in the figure for the variable-axis ('x' or 'y')
     9505        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
     9506          [timen]: name of the time variable ('WRFtime' for WRF times)
     9507          [units]: units string according to CF conventions ([tunits] since
     9508            [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
     9509          [kind]: kind of output
     9510            'Nval': according to a given number of values as 'Nval',[Nval]
     9511            'exct': according to an exact time unit as 'exct',[tunit];
     9512              tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
     9513               'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
     9514               'l': milisecond
     9515          [tfmt]: desired format
     9516          [label]: label at the graph ('!' for spaces)
     9517        [dimvfmt]=[[dvs],[dvf],[Ndv],[ordv]: format of the values at variable-axis
     9518          (or 'auto')
     9519          [dvs]: style of variable-axis ('auto' for 'pretty')
     9520            'Nfix', values computed at even 'Ndx'
     9521            'Vfix', values computed at even 'Ndx' increments
     9522            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
     9523          [dvf]: format of the labels at the var-axis ('auto' for '%5g')
     9524          [Ndv]: Number of ticks at the var-axis ('auto' for 5)
     9525          [ordv]: angle of orientation of ticks at the var-axis ('auto' for horizontal)
     9526        [colorbarvals]=[colbarn],[fmtcolorbar],[orientation]
     9527          [colorbarn]: name of the color bar
     9528          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
     9529          [orientation]: orientation of the colorbar ('vertical' ['auto'], 'horizontal')
     9530          * NOTE: single 'auto' for 'rainbow,%6g,vertical'
     9531        [smin/axv]: minimum and maximum values for the shading or string for each for:
     9532          'Srange': for full range
     9533          'Saroundmean@val': for mean-xtrm,mean+xtrm where
     9534            xtrm = np.min(mean-min@val,max@val-mean)
     9535          'Saroundminmax@val': for min*val,max*val
     9536          'Saroundpercentile@val': for median-xtrm,median+xtrm where
     9537            xtrm = np.min(median-percentile_(val),percentile_(100-val)-median)
     9538          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
     9539          'Smedian@val': for -xtrm,xtrm where
     9540             xtrm = np.min(median-min@val,max@val-median)
     9541          'Spercentile@val': for -xtrm,xtrm where
     9542             xtrm = np.min(median-percentile_(val),percentile_(100-val)-median)
     9543        [discvals]= [type],[size],[lwidth],[lcol] characteristics of the points for the discrete field
     9544          [type]: type of point. Any marker from matoplib must be filled !
     9545          [size]: size of point
     9546          [lwidth]: width of the line around the point
     9547          [lcol]: color of the line around the point
     9548          'auto': for [type]='o', [size]=5, [lwdith]=0.25, [lcol]='#000000'
     9549        [figt]: title of the figure ('!' for spaces)
     9550        [kindfig]: kind of figure output (ps, png, pdf)
     9551        [reverse]: Transformation of the values
     9552          * 'transpose': reverse the axes (x-->y, y-->x)
     9553          * 'flip'@[x/y]: flip the axis x or y
     9554        [mapv]: map characteristics: [proj],[res]
     9555          see full documentation: http://matplotlib.org/basemap/
     9556          [proj]: projection
     9557            * 'cyl', cilindric
     9558            * 'lcc', lambert conformal
     9559          [res]: resolution:
     9560            * 'c', crude
     9561            * 'l', low
     9562            * 'i', intermediate
     9563            * 'h', high
     9564            * 'f', full
     9565        [close]: Whether figure should be finished or not
     9566    """
     9567
     9568    fname = 'draw_2D_shad_contdisc_time'
     9569    if values == 'h':
     9570        print fname + '_____________________________________________________________'
     9571        print draw_2D_shad_contdisc_time.__doc__
     9572        quit()
     9573
     9574    expectargs = '[vnamefs];[varaxis];[timevals];[dimxyfmt];[colorbarvals];' +       \
     9575       '[sminv],[smaxv];[discvals];[figt];[kindfig];[reverse];[mapv];[close]';;;
     9576
     9577    drw.check_arguments(fname,values,expectargs,';')
     9578
     9579    vnamesfig = values.split(';')[0]
     9580    varaxis = values.split(';')[1]
     9581    timevals = values.split(';')[2].split(',')
     9582    colorbarvals = values.split(';')[3]
     9583    shadminmax = values.split(';')[4]
     9584    discvals = values.split(';')[5]
     9585    figt = values.split(';')[6].replace('!',' ')
     9586    kindfig = values.split(';')[7]
     9587    revals = values.split(';')[8]
     9588    mapvalue = values.split(';')[9]
     9589    close = gen.Str_Bool(values.split(';')[10])
     9590
     9591    ncs = ncfiles.split('@')
     9592    filesinf = {}
     9593    ifile = 0
     9594    for nc in ncs:
     9595        if ifile == 0:
     9596            hfn = 'cont'
     9597        else:
     9598            hfn = 'disc'
     9599
     9600        filesinf[hfn + 'filen'] = nc.split(';')[0]
     9601        filesinf[hfn + 'varn'] = nc.split(';')[1]
     9602        filesinf[hfn + 'dimtvn'] = nc.split(';')[2]
     9603        filesinf[hfn + 'dimvvn'] = nc.split(';')[3]
     9604        filesinf[hfn + 'dimvals'] = nc.split(';')[4].replace('|',':')
     9605        ifile = ifile + 1
     9606
     9607    for hfn in ['cont', 'disc']:
     9608        ncfilen = filesinf[hfn + 'filen']
     9609        varn = filesinf[hfn + 'varn']
     9610        vdimtn = filesinf[hfn + 'dimtvn']
     9611        vdimvn = filesinf[hfn + 'dimvvn']
     9612        dimvals = filesinf[hfn + 'dimvals']
     9613
     9614        if not os.path.isfile(ncfilen):
     9615            print errormsg
     9616            print '  ' + fname + ": '" + hfn + "' file '" + ncfile + "' does " +     \
     9617              "not exist !!"
     9618            quit(-1)
     9619
     9620        onc = NetCDFFile(ncfilen, 'r')
     9621
     9622        if  not onc.variables.has_key(varn):
     9623            print errormsg
     9624            print '  ' + fname + ": '" + hfn + "' file '" + ncfile +                 \
     9625              "' does not have variable '" + varn + "' !!"
     9626            varns = sorted(onc.variables.keys())
     9627            print '    available ones:', varns
     9628            quit(-1)
     9629
     9630        # Variables' values
     9631        ovar = onc.variables[varn]
     9632        vals, dims = drw.slice_variable(ovar, dimvals.replace(',','|'))
     9633
     9634        # Values have to be 2D (might not be true for disc?)
     9635        if len(vals.shape) != 2:
     9636            print errormsg
     9637            print '  ' + fname + ": values have to be 2D !!"
     9638            print '    provided shape:', vals.shape, 'for slice:', dimvals
     9639            quit(-1)
     9640
     9641        if drw.searchInlist(ovar.ncattrs(),'units'):
     9642            varunits = ovar.getncattr('units')
     9643        else:
     9644            print warnmsg
     9645            print '  ' + fname + ": variable '" + varn + "' without units!!"
     9646            varunits = '-'
     9647
     9648        # dimensions
     9649        dimnamesv = [vdimtn, vdimvn]
     9650        if  not onc.variables.has_key(vdimtn):
     9651            print errormsg
     9652            print '  ' + fname + ": '" + hfn + "' file '" + ncfile +                 \
     9653              "' does not have time-dimension variable '" +  vdimtn + "' !!"
     9654            varns = sorted(onc.variables.keys())
     9655            print '    available ones:', varns
     9656            quit(-1)
     9657        if  not onc.variables.has_key(vdimvn):
     9658            print errormsg
     9659            print '  ' + fname + ": '" + hfn + "' file '" + ncfile +                 \
     9660              "' does not have var-dimension variable '" +  vdimvn + "' !!"
     9661            varns = sorted(onc.variables.keys())
     9662            print '    available ones:', varns
     9663            quit(-1)
     9664
     9665        if vdimtn != 'WRFtime':
     9666            objdimt = onc.variables[vdimtn]
     9667            if drw.searchInlist(objdimt.ncattrs(),'units'):
     9668                odimxu = objdimt.getncattr('units')
     9669            else:
     9670                print erormsg
     9671                print '  ' +fname+": variable dimension '" + vdimxn + "' without units!!"
     9672                quit(-1)
     9673        else:
     9674            objdimt = onc.variables['Times']
     9675            timev, odimxv = ncvar.compute_WRFtime(objdimt[:])
     9676            odimxv0 = ncvar.nericNCvariable_Dict({'time':timev.shape[0]}, ['time'],  \
     9677              'time', 'time', 'Time', tunits, timev)
     9678
     9679        objdimv = onc.variables[vdimvn]
     9680        if drw.searchInlist(objdimv.ncattrs(),'units'):
     9681            odimvu = objdimv.getncattr('units')
     9682        else:
     9683            print warnmsg
     9684            print '  ' +fname+": variable dimension '" + vdimvn + "' without units!!"
     9685            odimvu = '-'
     9686
     9687        odimxv0, odimyv0 = drw.dxdy_lonlatDIMS(objdimt[:], objdimv[:],                \
     9688          objdimt.dimensions, objdimv.dimensions, dimvals.replace(':','|').split(','),\
     9689          False)
     9690
     9691        # Some statistics
     9692        vn = np.min(vals)
     9693        vx = np.max(vals)
     9694        dxn = np.min(odimxv0)
     9695        dxx = np.max(odimxv0)
     9696        dyn = np.min(odimyv0)
     9697        dyx = np.max(odimyv0)
     9698
     9699        print gen.infmsg
     9700        print '  ' + fname + ": for '" + hfn + "' ________"
     9701        print '    values min:', vn, 'max:', vx
     9702        print '    dimx min:', dxn, 'max:', dxx
     9703        print '    dimy min:', dyn, 'max:', dyx
     9704
     9705        if hfn == 'cont':
     9706            timeinfo = []
     9707
     9708            if varaxis == 'x':
     9709                # Getting the right shape of dimensions
     9710                if vals.shape[0] != odimxv0.shape[0]: gdimx = vals.shape[0]
     9711                else: gdimx = vals.shape[1]
     9712                gdimy = odimxv0.shape[0]
     9713            else:
     9714                if vals.shape[0] != odimxv0.shape[0]: gdimy = vals.shape[0]
     9715                else: gdimy = vals.shape[1]
     9716                gdimx = odimxv0.shape[0]
     9717
     9718            if vals.shape[1] == gdimx: contv = vals[:]
     9719            else: contv = vals[:].transpose()
     9720
     9721            odxv = np.zeros((gdimy,gdimx), dtype=np.float)
     9722            odyv = np.zeros((gdimy,gdimx), dtype=np.float)
     9723
     9724            if len(odimxv0.shape) == 2:
     9725                if odimxv0.shape[1] == gdimx: odxv = odimxv0[:]
     9726                else: odxv = odimxv0.transpose()
     9727            else:
     9728                for k in range(gdimy):
     9729                    odxv[k,:] = odimxv0[:]
     9730
     9731            if len(odimyv0.shape) == 2:
     9732                if odimyv0.shape[0] == gdimy: odyv = odimyv0[:]
     9733                else: odyv = odimyv0.transpose()
     9734            else:
     9735                for k in range(gdimx):
     9736                    odyv[:,k] = odimyv0[:]
     9737
     9738            diminfo = {}
     9739            diminfo['units'] = [odimxu, odimyu]
     9740            diminfo['names'] = [gen.variables_values(vdimxn)[0],                     \
     9741              gen.variables_values(vdimyn)[1]]
     9742            varinfo = {}
     9743            varinfo['units'] = varunits
     9744            varinfo['name'] = vnamesfig
     9745
     9746            # Absolute xtremes for the plot
     9747            absxn = dxn
     9748            absxx = dxx
     9749            absyn = dyn
     9750            absyx = dyx
     9751
     9752        else:
     9753            discvarv = []
     9754            discx0 = []
     9755            discy0 = []
     9756            # Assuming that discrete values are masked
     9757            if type(vals) != type(gen.mamat):
     9758                madiscvarv = ma.masked_array(vals)
     9759            else:
     9760                madiscvarv = vals
     9761
     9762            # Getting 2D dimensions
     9763            [d2Dx, d2Dy] = np.meshgrid(odimxv0, odimyv0)
     9764
     9765            # Checking for flipping axis?
     9766            if d2Dx.shape[0] == vals.shape[1] and d2Dx.shape[1] == vals.shape[0]:
     9767                print '  rotating values !'
     9768                #madiscvarv = madiscvarv.transpose()
     9769                d2Dx = d2Dx.transpose()
     9770                d2Dy = d2Dy.transpose()
     9771
     9772            # vals have to be 2D otherwise... ?
     9773            dx = madiscvarv.shape[1]
     9774            dy = madiscvarv.shape[0]
     9775            for j in range(dy):
     9776                for i in range(dx):
     9777                    if not madiscvarv.mask[j,i]:
     9778                        discvarv.append(madiscvarv[j,i])
     9779                        discx.append(d2Dx[j,i])
     9780                        discy.append(d2Dy[j,i])
     9781            if varaxis == 'x':
     9782                discx = list(discy0)
     9783                # Getting same times
     9784                discy = gen.coincident_CFtimes(discx0, timeinfo[], varunits)
     9785            else:
     9786                discx = list(discx0)
     9787                discy = list(discy0)
     9788
     9789
     9790            Ndiscvals = len(discvarv)
     9791            print '    getting:', Ndiscvals, 'discrete values'
     9792            discv = np.zeros((Ndiscvals,3), dtype=np.float)
     9793            discv[:,0] = discx[:]
     9794            discv[:,1] = discy[:]
     9795            discv[:,2] = discvarv[:]
     9796
     9797            # Absolute xtremes for the plot
     9798            absxn = np.max([absxn,dxn])
     9799            absxx = np.min([absxx,dxx])
     9800            absyn = np.max([absyn,dyn])
     9801            absyx = np.min([absyx,dyx])
     9802
     9803        onc.close()
     9804
     9805    graphnx = [absxn, absxx, absyn, absyx]
     9806    print 'limits of the graphic:', graphnx
     9807
     9808    shading_nx = []
     9809    if shadminmax.split(',')[0][0:1] != 'S':
     9810            shading_nx.append(np.float(shadminmax.split(',')[0]))
     9811    else:
     9812        shading_nx.append(shadminmax.split(',')[0])
     9813
     9814    if shadminmax.split(',')[1][0:1] != 'S':
     9815        shading_nx.append(np.float(shadminmax.split(',')[1]))
     9816    else:
     9817        shading_nx.append(shadminmax.split(',')[1])
     9818
     9819    if discvals == 'auto':
     9820        discinfo = ['o', 5., 0.25, '#000000']
     9821    else:
     9822        discinfo = [discvals.split(',')[0], np.float(discvals.split(',')[1]),        \
     9823          np.float(discvals.split(',')[2]), discvals.split(',')[3]]
     9824
     9825    if mapvalue == 'None': mapvalue = None
     9826
     9827    if colorbarvals == 'auto': colorbarvals = 'rainbow,auto,auto'
     9828    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvals,',')
     9829    colormapv = [colbarn, fmtcolbar, colbaror]
     9830
     9831    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimxyf,',')
     9832    xaxis = [xstyl, xaxf, Nxax, xaxor]
     9833    yaxis = [ystyl, yaxf, Nyax, yaxor]
     9834
     9835    if revals == 'None':
     9836        revals = None
     9837
     9838    drw.plot_2Dshad_obssim_time(discv, contv, odxv, odyv, revals, diminfo, xaxis,    \
     9839      yaxis, colormapv, shading_nx, varinfo, discinfo, mapvalue, graphnx, figt,      \
     9840      kindfig, close)
     9841
     9842    return
     9843
     9844contfilen = '/home/lluis/estudios/ChemGBsAs/tests/199501/simout_snddiags.nc;ta;time;pres;time|-1,pres|-1'
     9845discfilen = '/home/lluis/estudios/ChemGBsAs/tests/199501/obs/snd/UWyoming_snd_87576.nc;ta;time;pres;time|-1,pres|-1'
     9846
     9847values = 'ta:time,bottom_top:Vfix,auto,3600.,auto,Vfix,auto,50.,auto:auto:Srange,Srange:o,5.:obs!&!sim!Ezeiza!airport!sounding:pdf:flip@y:None:yes'
     9848
     9849draw_2D_shad_contdisc_time(contfilen+'@'+discfilen, values, axfig=None, fig=None)
     9850
     9851
    91289852#quit()
    91299853
    … …  
    91459869
    91469870# Not checking file operation
    9147 Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_2cont',                        \
    9148   'draw_2D_shad_cont_time',                                                          \
     9871Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_contdisc',                     \
     9872  'draw_2D_shad_2cont', 'draw_2D_shad_cont_time',                                    \
    91499873  'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_2lines', 'draw_2lines_time',  \
    91509874  'draw_bar', 'draw_bar_line', 'draw_bar_line_time', 'draw_bar_time',                \
    … …  
    92049928    elif oper == 'draw_2D_shad_cont':
    92059929        draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
     9930    elif oper == 'draw_2D_shad_contdisc':
     9931        draw_2D_shad_contdisc(opts.ncfile, opts.values)
    92069932    elif oper == 'draw_2D_shad_2cont':
    92079933        draw_2D_shad_2cont(opts.ncfile, opts.values, opts.varname)

  • trunk/tools/drawing_tools.py

    r1946 r1947  
    16301630    fname = 'transform'
    16311631
    1632     print fname + ' Lluis: input values trans:', trans, 'dxv:', dxv, 'dyv:', dyv, 'dxt', dxt, 'dyt', dyt, 'dxl', dxl, \
    1633       'dyl', dyl, 'dxtit', dxtit, 'dytit', dytit, 'axxkind', axxkind, 'axykind', axykind
     1632    #print fname + ' Lluis: input values trans:', trans, 'dxv:', dxv, 'dyv:', dyv, 'dxt', dxt, 'dyt', dyt, 'dxl', dxl, \
     1633    #  'dyl', dyl, 'dxtit', dxtit, 'dytit', dytit, 'axxkind', axxkind, 'axykind', axykind
    16341634
    16351635    transforms = trans.split('|')
    … …  
    16881688        newdxtit = str(dxtit)
    16891689        wdTx = True
    1690         print fname + ': Lluis hey:', newdxtit
     1690    #    print fname + ': Lluis hey:', newdxtit
    16911691    else:
    16921692        newdxtit = None
    … …  
    16961696        newdytit = str(dytit)
    16971697        wdTy = True
    1698         print fname + ': Lluis hey:', newdytit
     1698    #    print fname + ': Lluis hey:', newdytit
    16991699    else:
    17001700        newdytit = None
    … …  
    17811781            quit(-1)
    17821782
    1783     print fname + ' Lluis: newdxv:', newdxv, 'newdyv:', newdyv, 'newdxt:', newdxt, 'newdyt:', newdyt, \
    1784       'newdxl', newdxl, 'newdyl', newdyl, 'newdxtit', newdxtit, 'newdytit', newdytit
     1783    #print fname + ' Lluis: newdxv:', newdxv, 'newdyv:', newdyv, 'newdxt:', newdxt, 'newdyt:', newdyt, \
     1784    #  'newdxl', newdxl, 'newdyl', newdyl, 'newdxtit', newdxtit, 'newdytit', newdytit
    17851785
    17861786    return newvals, newdxv, newdyv, newdxt, newdyt, newdxl, newdyl, newdxtit, newdytit
    … …  
    53695369    return lonv,latv
    53705370
    5371 def dxdy_lonlatDIMS(dxv,dyv,dnx,dny,dd):
     5371def dxdy_lonlatDIMS(dxv,dyv,dnx,dny,dd,error=True):
    53725372    """ Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values for a given
    53735373      list of values
    … …  
    53805380        [dimname]: name of the dimension
    53815381        [val]: value (-1 for all the range)
     5382      error: whether script should quit if final shapes do not coincide
    53825383    """
    53835384    fname = 'dxdy_lonlatDIMS'
    … …  
    54295430        print '  ' + fname + ': dimension size on x:', len(lonv.shape), 'and on y:', \
    54305431          len(latv.shape),'do not coincide!!'
    5431         quit(-1)
     5432        if error: quit(-1)
    54325433
    54335434    return lonv,latv
    … …  
    1146611467    dp = (psfc-pmin)/(dz-1)
    1146711468   
     11469    ws = SatMixingRatio(tsfc+273.15,psfc)
     11470    thetasfc = (tsfc+273.15)*(psfc/p0)**(R/Cp)
    1146811471    for iz in range(dz):
    1146911472        ethetapotv[iz,0] = psfc-dp*iz
    1147011473        thetapotv = (tsfc+273.15)*((psfc-dp*iz)/p0)**(R/Cp)
    11471         ws = SatMixingRatio(tsfc+273.15,psfc)
    11472         ethetapotv[iz,1] = thetapotv*np.exp(L*1000.*ws/(Cp*(tsfc+273.15))) -273.15
    11473         print iz,psfc-dp*iz,thetapotv-273.15,ethetapotv[iz,1]
    11474 
     11474        tfromthetasfc = thetasfc*(p0/(psfc-dp*iz))**(-R/Cp)
     11475        wstfromthetasfc = SatMixingRatio(tfromthetasfc,psfc-dp*iz)
     11476        #ethetapotv[iz,1] = thetapotv*np.exp(L*1000.*ws/(Cp*(tsfc+273.15))) -273.15
     11477        ethetapotsfcv = thetapotv*np.exp(L*1000.*wstfromthetasfc/(Cp*(tsfc+273.15))) -273.15
     11478        ethetapotv[iz,1] = ethetapotsfcv
     11479        print iz,psfc-dp*iz,thetapotv-273.15,tfromthetasfc-273.15,':',ethetapotv[iz,1],wstfromthetasfc,ethetapotsfcv
     11480
     11481    print thetasfc,ws
    1147511482    #quit()
    1147611483    return ethetapotv
    … …  
    1163411641        ax.semilogy(dryt[:,1], dryt[:,0]/100., color='#FFCCCC', linewidth=0.75)
    1163511642
    11636     # Including moist- adiabats
     11643    # Including moist- adiabats (not working yet)
    1163711644    #for itt in range(ddt):
    1163811645    #    moist = ethetapotct(xmin+itt*xfrqtk, pmin=10000.)
    … …  
    1171911726
    1172011727    return
     11728
     11729def plot_2Dshad_obssim(obsv, varsv, dimxv, dimyv, reva, diminf, xaxv, yaxv, cbarv,   \
     11730  vs, varinf, discinf, mapv, gmaxmin, figtitle, kfig, close):
     11731    """ Function to plot a 2D shadding plot with comparison between observations (as
     11732     filled colored circles) and simulation (shadded continuous field)
     11733      obsv= 3-column with observational values (x, y, value)
     11734      varsv= 2D matrix of simulation values
     11735      dimxv= values for x-coordinate from simulation
     11736      dimyv= values for y-coordinate from simulation
     11737      reva= in case of modification of the plotting of 2D variables ('|' can be used for combination)
     11738        * 'transpose': reverse the axes (x-->y, y-->x)
     11739        * 'flip'@[x/y]: flip the axis x or y
     11740      diminf= dictionary with the information from the dimensions in the plot
     11741        diminf['units']= [dimxu, dimyu]: units at the axes of x and y
     11742        diminf['names']= [dimxn, dimyn]: names of the dimension
     11743      xaxv= list with the x-axis paramteres [style, format, number and orientation]
     11744      yaxv= list with the y-axis paramteres [style, format, number and orientation]
     11745      cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror]
     11746        colorbar: name of the color bar to use
     11747        cbarfmt: format of the numbers in the colorbar
     11748        cbaror: orientation of the colorbar
     11749      vs= list with minmum and maximum values to plot in shadow or one strings forminimum and maximum for:
     11750        'Srange': for full range
     11751        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
     11752        'Saroundminmax@val': for min*val,max*val
     11753        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
     11754          percentile_(100-val)-median)
     11755        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
     11756        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
     11757        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
     11758           percentile_(100-val)-median)
     11759      varinf= dictionary with the information of the variable in the plot
     11760        varinf['units']= units of the variable to shadow
     11761        varinf['name']= name of the variable
     11762      discinf= [[type],[size],[lwidth],[lcol]] characteristics of the points for the discrete field
     11763        [type]: type of point. Any from matoplib must be filled !
     11764        [size]: size of point
     11765        [lwidht]: width of the line around the marker
     11766        [lcol]: color of the line around the marker
     11767      mapv= map characteristics: [proj],[res]
     11768        see full documentation: http://matplotlib.org/basemap/
     11769        [proj]: projection
     11770          * 'cyl', cilindric
     11771          * 'lcc', lambert conformal
     11772        [res]: resolution:
     11773          * 'c', crude
     11774          * 'l', low
     11775          * 'i', intermediate
     11776          * 'h', high
     11777          * 'f', full
     11778      gmaxmin= [xaxismin, xaxismax, yaxismin, yaxismax] minimum and maximum values of
     11779        the axis of the figure
     11780      figtitle= title of the figure
     11781      kifg= kind of figure file output (png, ps, pdf)
     11782      close= whether figure should be closed or not
     11783    """
     11784    fname = 'plot_2Dshad_obssim'
     11785
     11786    if len(varsv.shape) != 2:
     11787        print errormsg
     11788        print '  ' + fname + ': wrong variable shadow rank:', varsv.shape,           \
     11789          'is has to be 2D!!'
     11790        quit(-1)
     11791
     11792    # Getting the right lon values for plotting
     11793    if mapv is not None:
     11794        dimxv = np.where(dimxv > 180., dimxv-360., dimxv)
     11795
     11796    # Axis ticks
     11797    # Usually axis > x must be the lon, thus...
     11798    dimxv0 = dimxv.copy()
     11799    dimyv0 = dimyv.copy()
     11800
     11801    dxn = dimxv.min()
     11802    dxx = dimxv.max()
     11803    dyn = dimyv.min()
     11804    dyx = dimyv.max()
     11805
     11806    if xaxv[0] == 'pretty':
     11807        dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2]))
     11808    elif xaxv[0] == 'Nfix':
     11809        dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2]))
     11810    elif xaxv[0] == 'Vfix':
     11811        dimxt0 = np.arange(xaxv[2]*int(dxn/xaxv[2]),dxx+xaxv[2],xaxv[2])
     11812    if yaxv[0] == 'pretty':
     11813        dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2]))
     11814    elif yaxv[0] == 'Nfix':
     11815        dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2]))
     11816    elif yaxv[0] == 'Vfix':
     11817        dimyt0 = np.arange(yaxv[2]*int(dyn/yaxv[2]),dyx+yaxv[2],yaxv[2])
     11818
     11819    dimxl0 = []
     11820    for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], style=xaxv[1]))
     11821    dimyl0 = []
     11822    for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1]))
     11823
     11824    dimn = diminf['names']
     11825    dimxu = diminf['units'][0]
     11826    dimyu = diminf['units'][1]
     11827
     11828    dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')'
     11829    dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')'
     11830
     11831    pixkind = 'data'
     11832
     11833    print 'reva:', reva
     11834    if reva is not None:
     11835        varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT, limx, limy =      \
     11836          gen.transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0,     \
     11837          dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0)
     11838    else:
     11839        dimxv = dimxv0
     11840        dimyv = dimyv0
     11841        dimxt = dimxt0
     11842        dimyt = dimyt0
     11843        dimxl = dimxl0
     11844        dimyl = dimyl0
     11845        dimxT = dimxT0
     11846        dimyT = dimyT0
     11847        limx = [dxn, dxx]
     11848        limy = [dyn, dyx]
     11849
     11850    if not mapv is None:
     11851        map_proj=mapv.split(',')[0]
     11852        map_res=mapv.split(',')[1]
     11853
     11854        lon0 = dimxv.copy()
     11855        lat0 = dimyv.copy()
     11856
     11857        dx = lon0.shape[1]
     11858        dy = lat0.shape[0]
     11859
     11860        nlon = np.min(lon0)
     11861        xlon = np.max(lon0)
     11862        nlat = np.min(lat0)
     11863        xlat = np.max(lat0)
     11864
     11865        lon2 = lon0[dy/2,dx/2]
     11866        lat2 = lat0[dy/2,dx/2]
     11867
     11868        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
     11869          xlon, ',', xlat
     11870
     11871        if map_proj == 'cyl':
     11872            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
     11873              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
     11874        elif map_proj == 'lcc':
     11875            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
     11876              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
     11877        else:
     11878            print errormsg
     11879            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
     11880            print '    available: cyl, lcc'
     11881            quit(-1)
     11882 
     11883        x,y = m(lon0,lat0)
     11884
     11885    else:
     11886        # No following data values
     11887        x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv))
     11888        y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv))
     11889
     11890    # Changing limits of the colors
     11891    # Combining values
     11892    fvarsv = varsv.flatten()
     11893    fobsv = obsv[:,2].flatten()
     11894    allvals = np.array(list(fvarsv) + list(fobsv))
     11895
     11896    vsend = graphic_range(vs, allvals)
     11897
     11898    plt.rc('text',usetex=True)
     11899
     11900    # 2D continuous shaded field
     11901    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], vmax=vsend[1])
     11902    # in case using scatter
     11903    #cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2])
     11904
     11905    # Adding observations
     11906    Nobs = obsv.shape[0]
     11907
     11908    # Normalize into x,y space
     11909    if mapv is None:
     11910        xmin = np.min(dimxv)
     11911        xmax = np.max(dimxv)
     11912        ymin = np.min(dimyv)
     11913        ymax = np.max(dimyv)
     11914    else:
     11915        xmin = np.min(lon0)
     11916        xmax = np.max(lon0)
     11917        ymin = np.min(lat0)
     11918        ymax = np.max(lat0)
     11919
     11920    # Normalize colors
     11921    my_cmap = plt.cm.get_cmap(cbarv[0])
     11922    norm = mpl.colors.Normalize(vsend[0],vsend[1])       
     11923
     11924    for iobs in range(Nobs):
     11925        xobs = (obsv[iobs,0]-xmin)/(xmax-xmin)
     11926        yobs = (obsv[iobs,1]-ymin)/(ymax-ymin)
     11927        clv = my_cmap(norm(obsv[iobs,2]))
     11928        plt.plot(xobs, yobs, discinf[0], color=clv, markersize=discinf[1], zorder=100)
     11929        plt.plot(xobs, yobs, discinf[0], color=discinf[3], linewidth=discinf[2],     \
     11930          markersize=discinf[1], zorder=100, fillstyle='none')
     11931
     11932    # with scatter... same result!
     11933    #xobs = (obsv[:,0]-xmin)/(xmax-xmin)
     11934    #yobs = (obsv[:,1]-ymin)/(ymax-ymin)
     11935    #clv = my_cmap(norm(obsv[:,2]))
     11936    #plt.scatter(xobs, yobs, s=discinf[1]*10., c=clv, zorder=100)
     11937
     11938    if cbarv[2] == 'horizontal':
     11939        cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2])
     11940        # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib
     11941        ticklabels= cbar.ax.get_xticklabels()
     11942        Nticks = len(ticklabels)
     11943        ticklabs = []
     11944        for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text())
     11945        cbar.ax.set_xticklabels(ticklabs,rotation=90)
     11946    else:
     11947        cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2])
     11948
     11949    if not mapv is None:
     11950        if cbarv[0] == 'gist_gray':
     11951            m.drawcoastlines(color="red")
     11952        else:
     11953            m.drawcoastlines()
     11954
     11955        meridians = gen.pretty_int(nlon,xlon,xaxv[2])
     11956        m.drawmeridians(meridians,labels=[True,False,False,True])
     11957                                                                                                             
     11958        parallels = gen.pretty_int(nlat,xlat,yaxv[2])
     11959        m.drawparallels(parallels,labels=[False,True,True,False])
     11960
     11961        plt.xlabel('W-E')
     11962        plt.ylabel('S-N')
     11963    else:
     11964        plt.xlabel(dimxT)
     11965        plt.ylabel(dimyT)
     11966
     11967    if mapv is None:
     11968        ldimxt = list(dimxt)
     11969        ldimyt = list(dimyt)
     11970        for val in ldimxt:
     11971            if val < np.min(limx): ldimxt.remove(val)
     11972            if val > np.max(limx): ldimxt.remove(val)
     11973        for val in ldimyt:
     11974            if val < np.min(limy): ldimyt.remove(val)
     11975            if val > np.max(limy): ldimyt.remove(val)
     11976        dimxt = np.array(ldimxt)
     11977        dimyt = np.array(ldimyt)
     11978        dxt = gen.interpolate_locs(dimxt,dimxv[0,:],'lin')
     11979        dyt = gen.interpolate_locs(dimyt,dimyv[:,0],'lin')
     11980        dimxt = dxt/varsv.shape[1]
     11981        dimyt = dyt/varsv.shape[0]
     11982
     11983        # Using output from transform
     11984        gxn = (gmaxmin[0]-xmin)/(xmax-xmin)
     11985        gxx = (gmaxmin[1]-xmin)/(xmax-xmin)
     11986        gyn = (gmaxmin[2]-ymin)/(ymax-ymin)
     11987        gyx = (gmaxmin[3]-ymin)/(ymax-ymin)
     11988
     11989        #gxn = (limx[0]-xmin)/(xmax-xmin)
     11990        #gxx = (limx[1]-xmin)/(xmax-xmin)
     11991        #gyn = (limy[0]-ymin)/(ymax-ymin)
     11992        #gyx = (limy[1]-ymin)/(ymax-ymin)
     11993
     11994        # Are axis flipped?
     11995        signxl = np.abs(limx[0]-limx[1])/(limx[0]-limx[1])
     11996        signxt = np.abs(ldimxt[0]-ldimxt[1])/(ldimxt[0]-ldimxt[1])
     11997        signyl = np.abs(limy[0]-limy[1])/(limy[0]-limy[1])
     11998        signyt = np.abs(ldimyt[0]-ldimyt[1])/(ldimyt[0]-ldimyt[1])
     11999
     12000        if signxl != signxt:
     12001            plt.xticks(1.-dimxt, dimxl, rotation=xaxv[3])
     12002        else:
     12003            plt.xticks(dimxt, dimxl, rotation=xaxv[3])
     12004        if signyl != signyt:
     12005            plt.yticks(1.-dimyt, dimyl, rotation=yaxv[3])
     12006        else:
     12007            plt.yticks(dimyt, dimyl, rotation=yaxv[3])
     12008
     12009        # Only if gmaxmin is used to compute limits of the figure !!
     12010        if signxl != signxt:
     12011            newgxn = gxx
     12012            newgxx = gxn
     12013        else:
     12014            newgxn = gxn
     12015            newgxx = gxx
     12016        if signyl != signyt:
     12017            newgyn = gyx
     12018            newgyx = gyn
     12019        else:
     12020            newgyn = gyn
     12021            newgyx = gyx
     12022
     12023        plt.axis([newgxn, newgxx, newgyn, newgyx])
     12024        # otherwise
     12025        #plt.axis([gxn, gxx, gyn, gyx])
     12026
     12027# units labels
     12028    cbar.set_label(gen.latex_text(varinf['name']) + ' (' +                           \
     12029      units_lunits(varinf['units']) + ')')
     12030   
     12031    plt.title(gen.latex_text(figtitle))
     12032
     12033    figname = '2Dshad_obs-sim_comparison.png'
     12034    output_kind(kfig, figname, close)
     12035
     12036    return
     12037
     12038def plot_2Dshad_obssim_time(obsv, varsv, dimvv, dimtv, dataxis, reva, diminf, vaxv,  \
     12039  timeinf, cbarv, vs, varinf, discinf, mapv, gmaxmin, figtitle, kfig, close):
     12040    """ Function to plot a 2D shadding plot with comparison between observations (as
     12041     filled colored circles) and simulation (shadded continuous field)
     12042      obsv= 3-column with observational values (x, y, value)
     12043      varsv= 2D matrix of simulation values
     12044      dimvv= values for v-coordinate from simulation
     12045      dimtv= values for t-coordinate from simulation
     12046      dataxis= which axis 'x' or 'y' should be use for the non-temporal dimension
     12047      reva= in case of modification of the plotting of 2D variables ('|' can be used for combination)
     12048        * 'transpose': reverse the axes (x-->y, y-->x)
     12049        * 'flip'@[x/y]: flip the axis x or y
     12050      diminf= dictionary with the information from the dimensions in the plot
     12051        diminf['units']= [dimxu, dimyu]: units at the axes of x and y
     12052        diminf['names']= [dimxn, dimyn]: names of the dimension
     12053      vaxv= list with the v-axis paramteres [style, format, number and orientation]
     12054      timeinf= list with information for the time-axis
     12055        [timepos]: positions of the markers
     12056        [timelabs]: labels of the markers
     12057        [ttitle]: title in graphic
     12058      cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror]
     12059        colorbar: name of the color bar to use
     12060        cbarfmt: format of the numbers in the colorbar
     12061        cbaror: orientation of the colorbar
     12062      vs= list with minmum and maximum values to plot in shadow or one strings forminimum and maximum for:
     12063        'Srange': for full range
     12064        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
     12065        'Saroundminmax@val': for min*val,max*val
     12066        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
     12067          percentile_(100-val)-median)
     12068        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
     12069        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
     12070        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
     12071           percentile_(100-val)-median)
     12072      varinf= dictionary with the information of the variable in the plot
     12073        varinf['units']= units of the variable to shadow
     12074        varinf['name']= name of the variable
     12075      discinf= [[type],[size],[lwidth],[lcol]] characteristics of the points for the discrete field
     12076        [type]: type of point. Any from matoplib must be filled !
     12077        [size]: size of point
     12078        [lwidht]: width of the line around the marker
     12079        [lcol]: color of the line around the marker
     12080      mapv= map characteristics: [proj],[res]
     12081        see full documentation: http://matplotlib.org/basemap/
     12082        [proj]: projection
     12083          * 'cyl', cilindric
     12084          * 'lcc', lambert conformal
     12085        [res]: resolution:
     12086          * 'c', crude
     12087          * 'l', low
     12088          * 'i', intermediate
     12089          * 'h', high
     12090          * 'f', full
     12091      gmaxmin= [xaxismin, xaxismax, yaxismin, yaxismax] minimum and maximum values of
     12092        the axis of the figure
     12093      figtitle= title of the figure
     12094      kifg= kind of figure file output (png, ps, pdf)
     12095      close= whether figure should be closed or not
     12096    """
     12097    fname = 'plot_2Dshad_obssim'
     12098
     12099    if len(varsv.shape) != 2:
     12100        print errormsg
     12101        print '  ' + fname + ': wrong variable shadow rank:', varsv.shape,           \
     12102          'is has to be 2D!!'
     12103        quit(-1)
     12104
     12105    # Getting the right lon values for plotting
     12106    if mapv is not None:
     12107        dimxv = np.where(dimxv > 180., dimxv-360., dimxv)
     12108
     12109    # Axis ticks
     12110    # Usually axis > x must be the time, thus...
     12111    dimxv0 = dimtv.copy()
     12112    dimyv0 = dimyv.copy()
     12113
     12114    dxn = dimtv.min()
     12115    dxx = dimtv.max()
     12116    dyn = dimyv.min()
     12117    dyx = dimyv.max()
     12118
     12119    dimxt0 = timeinf[0]
     12120    dimxl0 = timeinf[1]
     12121
     12122    if vaxv[0] == 'pretty':
     12123        dimyt0 = np.array(gen.pretty_int(dyn,dyx,vaxv[2]))
     12124    elif vaxv[0] == 'Nfix':
     12125        dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*vaxv[2]))
     12126    elif vaxv[0] == 'Vfix':
     12127        dimyt0 = np.arange(vaxv[2]*int(dyn/vaxv[2]),dyx+vaxv[2],vaxv[2])
     12128
     12129    dimyl0 = []
     12130    for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1]))
     12131
     12132    dimyu = diminf['units'][1]
     12133
     12134    dimxT0 = timeinf[2]
     12135    dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')'
     12136
     12137    pixkind = 'data'
     12138
     12139    if dataxis == 'x':
     12140        if reva is not None:
     12141            reva = 'transpose|' + reva
     12142        else:
     12143            reva = 'transpose'
     12144
     12145    if reva is not None:
     12146        varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT, limx, limy =      \
     12147          gen.transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0,     \
     12148          dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0)
     12149    else:
     12150        dimxv = dimxv0
     12151        dimyv = dimyv0
     12152        dimxt = dimxt0
     12153        dimyt = dimyt0
     12154        dimxl = dimxl0
     12155        dimyl = dimyl0
     12156        dimxT = dimxT0
     12157        dimyT = dimyT0
     12158        limx = [dxn, dxx]
     12159        limy = [dyn, dyx]
     12160
     12161    if not mapv is None:
     12162        map_proj=mapv.split(',')[0]
     12163        map_res=mapv.split(',')[1]
     12164
     12165        lon0 = dimxv.copy()
     12166        lat0 = dimyv.copy()
     12167
     12168        dx = lon0.shape[1]
     12169        dy = lat0.shape[0]
     12170
     12171        nlon = np.min(lon0)
     12172        xlon = np.max(lon0)
     12173        nlat = np.min(lat0)
     12174        xlat = np.max(lat0)
     12175
     12176        lon2 = lon0[dy/2,dx/2]
     12177        lat2 = lat0[dy/2,dx/2]
     12178
     12179        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
     12180          xlon, ',', xlat
     12181
     12182        if map_proj == 'cyl':
     12183            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
     12184              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
     12185        elif map_proj == 'lcc':
     12186            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
     12187              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
     12188        else:
     12189            print errormsg
     12190            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
     12191            print '    available: cyl, lcc'
     12192            quit(-1)
     12193 
     12194        x,y = m(lon0,lat0)
     12195
     12196    else:
     12197        # No following data values
     12198        x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv))
     12199        y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv))
     12200
     12201    # Changing limits of the colors
     12202    # Combining values
     12203    fvarsv = varsv.flatten()
     12204    fobsv = obsv[:,2].flatten()
     12205    allvals = np.array(list(fvarsv) + list(fobsv))
     12206
     12207    vsend = graphic_range(vs, allvals)
     12208
     12209    plt.rc('text',usetex=True)
     12210
     12211    # 2D continuous shaded field
     12212    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], vmax=vsend[1])
     12213    # in case using scatter
     12214    #cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2])
     12215
     12216    # Adding observations
     12217    Nobs = obsv.shape[0]
     12218
     12219    # Normalize into x,y space
     12220    if mapv is None:
     12221        xmin = np.min(dimxv)
     12222        xmax = np.max(dimxv)
     12223        ymin = np.min(dimyv)
     12224        ymax = np.max(dimyv)
     12225    else:
     12226        xmin = np.min(lon0)
     12227        xmax = np.max(lon0)
     12228        ymin = np.min(lat0)
     12229        ymax = np.max(lat0)
     12230
     12231    # Normalize colors
     12232    my_cmap = plt.cm.get_cmap(cbarv[0])
     12233    norm = mpl.colors.Normalize(vsend[0],vsend[1])       
     12234
     12235    for iobs in range(Nobs):
     12236        xobs = (obsv[iobs,0]-xmin)/(xmax-xmin)
     12237        yobs = (obsv[iobs,1]-ymin)/(ymax-ymin)
     12238        clv = my_cmap(norm(obsv[iobs,2]))
     12239        plt.plot(xobs, yobs, discinf[0], color=clv, markersize=discinf[1], zorder=100)
     12240        plt.plot(xobs, yobs, discinf[0], color=discinf[3], linewidth=discinf[2],     \
     12241          markersize=discinf[1], zorder=100, fillstyle='none')
     12242
     12243    # with scatter... same result!
     12244    #xobs = (obsv[:,0]-xmin)/(xmax-xmin)
     12245    #yobs = (obsv[:,1]-ymin)/(ymax-ymin)
     12246    #clv = my_cmap(norm(obsv[:,2]))
     12247    #plt.scatter(xobs, yobs, s=discinf[1]*10., c=clv, zorder=100)
     12248
     12249    if cbarv[2] == 'horizontal':
     12250        cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2])
     12251        # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib
     12252        ticklabels= cbar.ax.get_xticklabels()
     12253        Nticks = len(ticklabels)
     12254        ticklabs = []
     12255        for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text())
     12256        cbar.ax.set_xticklabels(ticklabs,rotation=90)
     12257    else:
     12258        cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2])
     12259
     12260    if not mapv is None:
     12261        if cbarv[0] == 'gist_gray':
     12262            m.drawcoastlines(color="red")
     12263        else:
     12264            m.drawcoastlines()
     12265
     12266        meridians = gen.pretty_int(nlon,xlon,xaxv[2])
     12267        m.drawmeridians(meridians,labels=[True,False,False,True])
     12268                                                                                                             
     12269        parallels = gen.pretty_int(nlat,xlat,yaxv[2])
     12270        m.drawparallels(parallels,labels=[False,True,True,False])
     12271
     12272        plt.xlabel('W-E')
     12273        plt.ylabel('S-N')
     12274    else:
     12275        plt.xlabel(dimxT)
     12276        plt.ylabel(dimyT)
     12277
     12278    if mapv is None:
     12279        ldimxt = list(dimxt)
     12280        ldimyt = list(dimyt)
     12281        for val in ldimxt:
     12282            if val < np.min(limx): ldimxt.remove(val)
     12283            if val > np.max(limx): ldimxt.remove(val)
     12284        for val in ldimyt:
     12285            if val < np.min(limy): ldimyt.remove(val)
     12286            if val > np.max(limy): ldimyt.remove(val)
     12287        dimxt = np.array(ldimxt)
     12288        dimyt = np.array(ldimyt)
     12289        dxt = gen.interpolate_locs(dimxt,dimxv[0,:],'lin')
     12290        dyt = gen.interpolate_locs(dimyt,dimyv[:,0],'lin')
     12291        dimxt = dxt/varsv.shape[1]
     12292        dimyt = dyt/varsv.shape[0]
     12293
     12294        # Using output from transform
     12295        gxn = (gmaxmin[0]-xmin)/(xmax-xmin)
     12296        gxx = (gmaxmin[1]-xmin)/(xmax-xmin)
     12297        gyn = (gmaxmin[2]-ymin)/(ymax-ymin)
     12298        gyx = (gmaxmin[3]-ymin)/(ymax-ymin)
     12299
     12300        #gxn = (limx[0]-xmin)/(xmax-xmin)
     12301        #gxx = (limx[1]-xmin)/(xmax-xmin)
     12302        #gyn = (limy[0]-ymin)/(ymax-ymin)
     12303        #gyx = (limy[1]-ymin)/(ymax-ymin)
     12304
     12305        # Are axis flipped?
     12306        signxl = np.abs(limx[0]-limx[1])/(limx[0]-limx[1])
     12307        signxt = np.abs(ldimxt[0]-ldimxt[1])/(ldimxt[0]-ldimxt[1])
     12308        signyl = np.abs(limy[0]-limy[1])/(limy[0]-limy[1])
     12309        signyt = np.abs(ldimyt[0]-ldimyt[1])/(ldimyt[0]-ldimyt[1])
     12310
     12311        if signxl != signxt:
     12312            plt.xticks(1.-dimxt, dimxl, rotation=xaxv[3])
     12313        else:
     12314            plt.xticks(dimxt, dimxl, rotation=xaxv[3])
     12315        if signyl != signyt:
     12316            plt.yticks(1.-dimyt, dimyl, rotation=yaxv[3])
     12317        else:
     12318            plt.yticks(dimyt, dimyl, rotation=yaxv[3])
     12319
     12320        # Only if gmaxmin is used to compute limits of the figure !!
     12321        if signxl != signxt:
     12322            newgxn = gxx
     12323            newgxx = gxn
     12324        else:
     12325            newgxn = gxn
     12326            newgxx = gxx
     12327        if signyl != signyt:
     12328            newgyn = gyx
     12329            newgyx = gyn
     12330        else:
     12331            newgyn = gyn
     12332            newgyx = gyx
     12333
     12334        plt.axis([newgxn, newgxx, newgyn, newgyx])
     12335        # otherwise
     12336        #plt.axis([gxn, gxx, gyn, gyx])
     12337
     12338# units labels
     12339    cbar.set_label(gen.latex_text(varinf['name']) + ' (' +                           \
     12340      units_lunits(varinf['units']) + ')')
     12341   
     12342    plt.title(gen.latex_text(figtitle))
     12343
     12344    figname = '2Dshad_obs-sim_comparison_time.png'
     12345    output_kind(kfig, figname, close)
     12346
     12347    return

  • trunk/tools/generic_tools.py

    r1934 r1947  
    57635763    return angle
    57645764
    5765 def lonlat2D(lon,lat):
     5765def lonlat2D(lon, lat, error=True):
    57665766    """ Function to return lon, lat 2D matrices from any lon,lat matrix
    57675767      lon= matrix with longitude values
    57685768      lat= matrix with latitude values
     5769      error= whether should quit if there is an error
    57695770    """
    57705771    fname = 'lonlat2D'
    57715772
    5772     if len(lon.shape) != len(lat.shape):
     5773    if len(lon.shape) == len(lat.shape):
     5774        if len(lon.shape) == 3:
     5775            lonvv = lon[0,:,:]
     5776            latvv = lat[0,:,:]
     5777        elif len(lon.shape) == 2:
     5778            lonvv = lon[:]
     5779            latvv = lat[:]
     5780        elif len(lon.shape) == 1:
     5781            lonlatv = np.meshgrid(lon[:],lat[:])
     5782            lonvv = lonlatv[0]
     5783            latvv = lonlatv[1]
     5784    else:
    57735785        print errormsg
    57745786        print '  ' + fname + ': longitude values with shape:', lon.shape,            \
    57755787          'is different than latitude values with shape:', lat.shape, '(dif. size) !!'
    5776         quit(-1)
    5777 
    5778     if len(lon.shape) == 3:
    5779         lonvv = lon[0,:,:]
    5780         latvv = lat[0,:,:]
    5781     elif len(lon.shape) == 2:
    5782         lonvv = lon[:]
    5783         latvv = lat[:]
    5784     elif len(lon.shape) == 1:
    5785         lonlatv = np.meshgrid(lon[:],lat[:])
    5786         lonvv = lonlatv[0]
    5787         latvv = lonlatv[1]
     5788        if error: quit(-1)
     5789
     5790        # Not so weird case !!
     5791        print '    Trying to fix it!'
     5792        if len(lon.shape) == 3:
     5793            lonv0 = lon[0,:,:]
     5794        else:
     5795            lonv0 = lon.copy()
     5796        if len(lat.shape) == 3:
     5797            latv0 = lat[0,:,:]
     5798        else:
     5799            latv0 = lat.copy()
     5800        if len(lonv0.shape) == 1 and len(latv0.shape) == 2:
     5801            lonvv = np.zeros(latv0.shape, dtype=np.float)
     5802            if lonv0.shape[0] == latv0.shape[0]:
     5803                for k in range(latv0.shape[1]):
     5804                    lonvv[:,k] = lonv0[:]
     5805            else:
     5806                for k in range(latv0.shape[0]):
     5807                    lonvv[k,:] = lonv0[:]
     5808            latvv = latv0.copy()
     5809
     5810        if len(lonv0.shape) == 2 and len(latv0.shape) == 1:
     5811            lonvv = lonv0.copy()
     5812            latvv = np.zeros(lonv0.shape, dtype=np.float)
     5813            if lonv0.shape[0] == latv0.shape[0]:
     5814                for k in range(lonv0.shape[1]):
     5815                    latvv[:,k] = latv0[:]
     5816            else:
     5817                for k in range(lonv0.shape[0]):
     5818                    latvv[k,:] = latv0[:]
    57885819
    57895820    return lonvv, latvv
    57905821
    5791 def interpolate_locs(locs,coords,kinterp):
     5822def interpolate_locs_old(locs,coords,kinterp):
    57925823    """ Function to provide interpolate locations on a given axis
    57935824    interpolate_locs(locs,axis,kinterp)
    … …  
    58235854
    58245855    for iloc in range(Nlocs):
     5856        found = False
     5857        a = None
     5858        b = None
     5859        c = None
    58255860        for icor in range(Ncoords-1):
    58265861            if locs[iloc] < minc and dcoords > 0.:
    58275862                a = 0.
    5828                 b = 1. / (coords[1] - coords[0])
     5863                b = coords[1] - coords[0]
    58295864                c = coords[0]
    58305865            elif locs[iloc] > maxc and dcoords > 0.:
    58315866                a = (Ncoords-1)*1.
    5832                 b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
     5867                b = coords[Ncoords-1] - coords[Ncoords-2]
    58335868                c = coords[Ncoords-2]
    58345869            elif locs[iloc] < minc and dcoords < 0.:
    58355870                a = (Ncoords-1)*1.
    5836                 b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
     5871                b = coords[Ncoords-1] - coords[Ncoords-2]
    58375872                c = coords[Ncoords-2]
    58385873            elif locs[iloc] > maxc and dcoords < 0.:
    58395874                a = 0.
    5840                 b = 1. / (coords[1] - coords[0])
     5875                b = (coords[1] - coords[0])
    58415876                c = coords[0]
    58425877            elif locs[iloc] >= coords[icor] and locs[iloc] < coords[icor+1] and dcoords > 0.:
    58435878                a = icor*1.
    5844                 b = 1. / (coords[icor+1] - coords[icor])
     5879                b = coords[icor+1] - coords[icor]
    58455880                c = coords[icor]
    5846                 print coords[icor], locs[iloc], coords[icor+1], ':', icor, '->', a, b
    58475881            elif locs[iloc] <= coords[icor] and locs[iloc] > coords[icor+1] and dcoords < 0.:
    58485882                a = icor*1.
    5849                 b = 1. / (coords[icor+1] - coords[icor])
     5883                b = coords[icor+1] - coords[icor]
    58505884                c = coords[icor]
     5885        if a is not None: found = True
     5886        if not found:
     5887            print errormsg
     5888            print '  ' + fname + ": case not prepared !!"
     5889            print '    locs[iloc]:', locs[iloc], 'minc:', minc, 'maxc:', maxc, 'dcoords:', dcoords
     5890            quit(-1)
    58515891
    58525892        if kinterp == 'lin':
    5853             intlocs[iloc] = a + (locs[iloc] - c)*b
     5893            intlocs[iloc] = a + (locs[iloc] - c)/b
    58545894        else:
    58555895            print errormsg
    58565896            print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
    58575897            quit(-1)
     5898
     5899    return intlocs
     5900
     5901def interpolate_locs(locs,coords,kinterp):
     5902    """ Function to provide interpolate locations on a given axis
     5903    interpolate_locs(locs,axis,kinterp)
     5904      locs= locations to interpolate
     5905      coords= axis values with the reference of coordinates
     5906      kinterp: kind of interpolation
     5907        'lin': linear
     5908    >>> coordinates = np.arange((10), dtype=np.float)
     5909    >>> values = np.array([-1.2, 2.4, 5.6, 7.8, 12.0])
     5910    >>> interpolate_locs(values,coordinates,'lin')
     5911    [ -1.2   2.4   5.6   7.8  12. ]
     5912    >>> coordinates[0] = 0.5
     5913    >>> coordinates[2] = 2.5
     5914    >>> interpolate_locs(values,coordinates,'lin')
     5915    [ -3.4          1.93333333   5.6          7.8         12.        ]
     5916    >>> coordinates = -10+np.arange((10), dtype=np.float)
     5917    >>> values = -np.array([-1.2, 2.4, 5.6, 7.8, 12.0])
     5918    >>> interpolate_locs(values,coordinates,'lin')
     5919    [ 11.2   7.6   4.4   2.2  -2. ]
     5920    """
     5921
     5922    fname = 'interpolate_locs'
     5923
     5924    if type(locs) == type('S') and locs == 'h':
     5925        print fname + '_____________________________________________________________'
     5926        print interpolate_locs.__doc__
     5927        quit()
     5928
     5929    Nlocs = locs.shape[0]
     5930    Ncoords = coords.shape[0]
     5931
     5932    origsing = np.abs(coords[Ncoords-1]-coords[0])/(coords[Ncoords-1]-coords[0])
     5933
     5934    intlocs = np.zeros((Nlocs), dtype=np.float)
     5935    minc = np.min(coords)
     5936    maxc = np.max(coords)
     5937
     5938    sortc = list(coords)
     5939    sortc.sort()
     5940
     5941    for iloc in range(Nlocs):
     5942        a = None
     5943        refi = None
     5944        if locs[iloc] < minc:
     5945            a = 0
     5946            refi = 0
     5947            sing = -1.
     5948        elif locs[iloc] > maxc:
     5949            a = Ncoords-1
     5950            refi = Ncoords-2
     5951            sing = 1.
     5952        else:
     5953            for icor in range(Ncoords-1):
     5954                if locs[iloc] >= sortc[icor] and locs[iloc] < sortc[icor+1]:
     5955                    a = icor
     5956                    refi = icor
     5957                    sing = 1.
     5958        if a is None:
     5959            print errormsg
     5960            print '  ' + fname + ': locs[iloc]:', locs[iloc], 'minc:', minc, 'maxc:',\
     5961              maxc, 'not ready!!'
     5962            quit(-1)
     5963        b = sortc[refi+1]-sortc[refi]
     5964        #print 'locs[iloc]:', locs[iloc], 'a:', a, 'refi:', refi, 'b:', b, 'np.abs(sortc[refi]-locs[iloc]):', np.abs(sortc[refi]-locs[iloc])
     5965
     5966        if kinterp == 'lin':
     5967            intlocs[iloc] = a + sing*np.abs(sortc[a]-locs[iloc])/b
     5968        else:
     5969            print errormsg
     5970            print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
     5971            quit(-1)
     5972
     5973    if origsing < 0.: intlocs = Ncoords - intlocs
    58585974
    58595975    return intlocs
    … …  
    1230112417
    1230212418    # axis limits. They will be used to flip the axis in plot if necessary
    12303     if len(dxv.shape) == 1:
    12304         newxaxislim = np.array([dxv[0], dxv[len(dxv)-1]])
    12305     elif len(dxv.shape) == 2:
    12306         newxaxislim = np.array([dxv[0,0], dxv[0,dxv.shape[1]-1]])
    12307     else:
    12308         print errormsg
    12309         print '  ' + fname + ": wrong rank of data for axis 'x':", dxv.shape, '!!'
    12310         print '    must be of rank 2!!'
    12311         quit(-1)
    12312 
    12313     if len(dyv.shape) == 1:
    12314         newyaxislim = np.array([dyv[0], dyv[len(dyv)-1]])
    12315     elif len(dyv.shape) == 2:
    12316         newyaxislim = np.array([dyv[0,0], dyv[dyv.shape[0]-1],0])
    12317     else:
    12318         print errormsg
    12319         print '  ' + fname + ": wrong rank of data for axis 'y':", dyv.shape, '!!'
    12320         print '    must be of rank 2!!'
    12321         quit(-1)
     12419    newxaxislim = [np.min(dxv), np.max(dxv)]
     12420    newyaxislim = [np.min(dyv), np.max(dyv)]
    1232212421
    1232312422    for transform in transforms:
    … …  
    1237612475            flip = transform.split('@')[1]
    1237712476            if flip == 'x':
    12378                 if len(newxaxislim.shape) == 1:
    12379                     newxaxislim = newxaxislim[::-1]
    12380                 else:
    12381                     for iy in len(newxaxislim.shape[0]):
    12382                         newxaxislim[iy,:] = newxaxislim[iy,::-1]
     12477                oldv = list(newxaxislim)
     12478                newxaxislim[0] = oldv[1]
     12479                newxaxislim[1] = oldv[0]
    1238312480            elif flip == 'y':
    12384                 if len(newyaxislim.shape) == 1:
    12385                     newyaxislim = newyaxislim[::-1]
    12386                 else:
    12387                     for ix in len(newyaxislim.shape[1]):
    12388                         newyaxislim[:,ix] = newyaxislim[::-1,ix]
     12481                oldv = list(newyaxislim)
     12482                newyaxislim[0] = oldv[1]
     12483                newyaxislim[1] = oldv[0]
    1238912484            elif flip == 'z':
    1239012485                newvals = newvals[...,::-1,:,:]

  • trunk/tools/variables_values.dat

    r1940 r1947  
    335335presnivs, pres, air_pressure, 0., 103000., air|pressure, Pa, Blues, $pres$, pres
    336336pres, pres, air_pressure, 0., 103000., air|pressure, Pa, Blues, $pres$, pres
     337air_pressure, pres, air_pressure, 0., 103000., air|pressure, Pa, Blues, $pres$, pres
    337338lpres, pres, air_pressure, 0., 103000., air|pressure, Pa, Blues, $pres$, pres
    338339LPRES, pres, air_pressure, 0., 103000., air|pressure, Pa, Blues, $pres$, pres
Note: See TracChangeset for help on using the changeset viewer.