source: lmdz_wrf/trunk/tools/drawing_tools.py @ 425

Last change on this file since 425 was 422, checked in by lfita, 10 years ago

Changing title and dx,dy dimensions on 'plot_Trajectories'

File size: 181.4 KB
Line 
1# -*- coding: iso-8859-15 -*-
2#import pylab as plt
3# From http://stackoverflow.com/questions/13336823/matplotlib-python-error
4import numpy as np
5import matplotlib as mpl
6mpl.use('Agg')
7import matplotlib.pyplot as plt
8from mpl_toolkits.basemap import Basemap
9import os
10from netCDF4 import Dataset as NetCDFFile
11import nc_var_tools as ncvar
12
13errormsg = 'ERROR -- error -- ERROR -- error'
14warnmsg = 'WARNING -- waring -- WARNING -- warning'
15
16fillValue = 1.e20
17
18####### Funtions
19# searchInlist:
20# datetimeStr_datetime:
21# dateStr_date:
22# numVector_String:
23# timeref_datetime:
24# slice_variable:
25# interpolate_locs:
26# datetimeStr_conversion:
27# percendone:
28# netCDFdatetime_realdatetime:
29# file_nlines:
30# variables_values:
31# check_colorBar:
32# units_lunits:
33# ASCII_LaTeX:
34# pretty_int:
35# DegGradSec_deg:
36# intT2dt:
37# lonlat_values:
38# date_CFtime:
39# pot_values:
40# CFtimes_plot:
41# color_lines:
42# output_kind:
43# check_arguments:
44# Str_Bool:
45# plot_points:
46# plot_2Dfield:
47# plot_2Dfield_easy:
48# plot_topo_geogrid:
49# plot_topo_geogrid_boxes:
50# plot_2D_shadow:
51# plot_2D_shadow_time: Plotting a 2D field with one of the axes being time
52# plot_Neighbourghood_evol:Plotting neighbourghood evolution# plot_Trajectories
53# plot_2D_shadow_contour:
54# plot_2D_shadow_contour_time:
55# dxdy_lonlat: Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values
56# plot_2D_shadow_line:
57# plot_lines: Function to plot a collection of lines
58
59# From nc_var_tools.py
60def searchInlist(listname, nameFind):
61    """ Function to search a value within a list
62    listname = list
63    nameFind = value to find
64    >>> searInlist(['1', '2', '3', '5'], '5')
65    True
66    """
67    for x in listname:
68      if x == nameFind:
69        return True
70    return False
71
72def datetimeStr_datetime(StringDT):
73    """ Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
74    >>> datetimeStr_datetime('1976-02-17_00:00:00')
75    1976-02-17 00:00:00
76    """
77    import datetime as dt
78
79    fname = 'datetimeStr_datetime'
80
81    dateD = np.zeros((3), dtype=int)
82    timeT = np.zeros((3), dtype=int)
83
84    dateD[0] = int(StringDT[0:4])
85    dateD[1] = int(StringDT[5:7])
86    dateD[2] = int(StringDT[8:10])
87
88    trefT = StringDT.find(':')
89    if not trefT == -1:
90#        print '  ' + fname + ': refdate with time!'
91        timeT[0] = int(StringDT[11:13])
92        timeT[1] = int(StringDT[14:16])
93        timeT[2] = int(StringDT[17:19])
94
95    if int(dateD[0]) == 0:
96        print warnmsg
97        print '    ' + fname + ': 0 reference year!! changing to 1'
98        dateD[0] = 1 
99 
100    newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])
101
102    return newdatetime
103
104def dateStr_date(StringDate):
105  """ Function to transform a string date ([YYYY]-[MM]-[DD] format) to a date object
106  >>> dateStr_date('1976-02-17')
107  1976-02-17
108  """
109  import datetime as dt
110
111  dateD = StringDate.split('-')
112  if int(dateD[0]) == 0:
113    print warnmsg
114    print '    dateStr_date: 0 reference year!! changing to 1'
115    dateD[0] = 1
116  newdate = dt.date(int(dateD[0]), int(dateD[1]), int(dateD[2]))
117  return newdate
118
119def numVector_String(vec,char):
120    """ Function to transform a vector of numbers to a single string [char] separated
121    numVector_String(vec,char)
122      vec= vector with the numerical values
123      char= single character to split the values
124    >>> print numVector_String(np.arange(10),' ')
125    0 1 2 3 4 5 6 7 8 9
126    """
127    fname = 'numVector_String'
128
129    if vec == 'h':
130        print fname + '_____________________________________________________________'
131        print numVector_String.__doc__
132        quit()
133
134    Nvals = len(vec)
135
136    string=''
137    for i in range(Nvals):
138        if i == 0:
139            string = str(vec[i])
140        else:
141            string = string + char + str(vec[i])
142
143    return string
144
145def timeref_datetime(refd, timeval, tu):
146    """ Function to transform from a [timeval] in [tu] units from the time referece [tref] to datetime object
147    refd: time of reference (as datetime object)
148    timeval: time value (as [tu] from [tref])
149    tu: time units
150    >>> timeref = date(1949,12,1,0,0,0)
151    >>> timeref_datetime(timeref, 229784.36, hours)
152    1976-02-17 08:21:36
153    """
154    import datetime as dt
155    import numpy as np
156
157## Not in timedelta
158#    if tu == 'years':
159#        realdate = refdate + dt.timedelta(years=float(timeval))
160#    elif tu == 'months':
161#        realdate = refdate + dt.timedelta(months=float(timeval))
162    if tu == 'weeks':
163        realdate = refd + dt.timedelta(weeks=float(timeval))
164    elif tu == 'days':
165        realdate = refd + dt.timedelta(days=float(timeval))
166    elif tu == 'hours':
167        realdate = refd + dt.timedelta(hours=float(timeval))
168    elif tu == 'minutes':
169        realdate = refd + dt.timedelta(minutes=float(timeval))
170    elif tu == 'seconds':
171        realdate = refd + dt.timedelta(seconds=float(timeval))
172    elif tu == 'milliseconds':
173        realdate = refd + dt.timedelta(milliseconds=float(timeval))
174    else:
175          print errormsg
176          print '    timeref_datetime: time units "' + tu + '" not ready!!!!'
177          quit(-1)
178
179    return realdate
180
181def slice_variable(varobj, dimslice):
182    """ Function to return a slice of a given variable according to values to its
183      dimensions
184    slice_variable(varobj, dims)
185      varobj= object wit the variable
186      dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension
187        [value]:
188          * [integer]: which value of the dimension
189          * -1: all along the dimension
190          * [beg]:[end] slice from [beg] to [end]
191    """
192    fname = 'slice_variable'
193
194    if varobj == 'h':
195        print fname + '_____________________________________________________________'
196        print slice_variable.__doc__
197        quit()
198
199    vardims = varobj.dimensions
200    Ndimvar = len(vardims)
201
202    Ndimcut = len(dimslice.split('|'))
203    if Ndimcut == 0:
204        Ndimcut = 1
205        dimcut = list(dimslice)
206
207    dimsl = dimslice.split('|')
208
209    varvalsdim = []
210    dimnslice = []
211
212    for idd in range(Ndimvar):
213        found = False
214        for idc in range(Ndimcut):
215            dimcutn = dimsl[idc].split(':')[0]
216            dimcutv = dimsl[idc].split(':')[1]
217            if vardims[idd] == dimcutn:
218                posfrac = dimcutv.find('@')
219                if posfrac != -1:
220                    inifrac = int(dimcutv.split('@')[0])
221                    endfrac = int(dimcutv.split('@')[1])
222                    varvalsdim.append(slice(inifrac,endfrac))
223                    dimnslice.append(vardims[idd])
224                else:
225                    if int(dimcutv) == -1:
226                        varvalsdim.append(slice(0,varobj.shape[idd]))
227                        dimnslice.append(vardims[idd])
228                    elif int(dimcutv) == -9:
229                        varvalsdim.append(int(varobj.shape[idd])-1)
230                    else:
231                        varvalsdim.append(int(dimcutv))
232                found = True
233                break
234        if not found and not searchInlist(dimnslice,vardims[idd]):
235            varvalsdim.append(slice(0,varobj.shape[idd]))
236            dimnslice.append(vardims[idd])
237
238    varvalues = varobj[tuple(varvalsdim)]
239
240    return varvalues, dimnslice
241
242def interpolate_locs(locs,coords,kinterp):
243    """ Function to provide interpolate locations on a given axis
244    interpolate_locs(locs,axis,kinterp)
245      locs= locations to interpolate
246      coords= axis values with the reference of coordinates
247      kinterp: kind of interpolation
248        'lin': linear
249    >>> coordinates = np.arange((10), dtype=np.float)
250    >>> values = np.array([-1.2, 2.4, 5.6, 7.8, 12.0])
251    >>> interpolate_locs(values,coordinates,'lin')
252    [ -1.2   2.4   5.6   7.8  13. ]
253    >>> coordinates[0] = 0.5
254    >>> coordinates[2] = 2.5
255    >>> interpolate_locs(values,coordinates,'lin')
256    [ -3.4          1.93333333   5.6          7.8         13.        ]
257    """
258
259    fname = 'interpolate_locs'
260
261    if locs == 'h':
262        print fname + '_____________________________________________________________'
263        print interpolate_locs.__doc__
264        quit()
265
266    Nlocs = locs.shape[0]
267    Ncoords = coords.shape[0]
268
269    dcoords = coords[Ncoords-1] - coords[0]
270
271    intlocs = np.zeros((Nlocs), dtype=np.float)
272    minc = np.min(coords)
273    maxc = np.max(coords)
274
275    for iloc in range(Nlocs):
276        for icor in range(Ncoords-1):
277            if locs[iloc] < minc and dcoords > 0.:
278                a = 0.
279                b = 1. / (coords[1] - coords[0])
280                c = coords[0]
281            elif locs[iloc] > maxc and dcoords > 0.:
282                a = (Ncoords-1)*1.
283                b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
284                c = coords[Ncoords-2]
285            elif locs[iloc] < minc and dcoords < 0.:
286                a = (Ncoords-1)*1.
287                b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
288                c = coords[Ncoords-2]
289            elif locs[iloc] > maxc and dcoords < 0.:
290                a = 0.
291                b = 1. / (coords[1] - coords[0])
292                c = coords[0]
293            elif locs[iloc] >= coords[icor] and locs[iloc] < coords[icor+1] and dcoords > 0.:
294                a = icor*1.
295                b = 1. / (coords[icor+1] - coords[icor])
296                c = coords[icor]
297                print coords[icor], locs[iloc], coords[icor+1], ':', icor, '->', a, b
298            elif locs[iloc] <= coords[icor] and locs[iloc] > coords[icor+1] and dcoords < 0.:
299                a = icor*1.
300                b = 1. / (coords[icor+1] - coords[icor])
301                c = coords[icor]
302
303        if kinterp == 'lin':
304            intlocs[iloc] = a + (locs[iloc] - c)*b
305        else:
306            print errormsg
307            print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
308            quit(-1)
309
310    return intlocs
311
312def datetimeStr_conversion(StringDT,typeSi,typeSo):
313    """ Function to transform a string date to an another date object
314    StringDT= string with the date and time
315    typeSi= type of datetime string input
316    typeSo= type of datetime string output
317      [typeSi/o]
318        'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
319        'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
320        'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
321        'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
322        'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
323        'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
324        'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H],
325          [H], ':', [M], [M], ':', [S], [S]
326    >>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
327    [1976    2   17    8   32    5]
328    >>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
329    1980/03/05 18-00-00
330    """
331    import datetime as dt
332
333    fname = 'datetimeStr_conversion'
334
335    if StringDT[0:1] == 'h':
336        print fname + '_____________________________________________________________'
337        print datetimeStr_conversion.__doc__
338        quit()
339
340    if typeSi == 'cfTime':
341        timeval = np.float(StringDT.split(',')[0])
342        tunits = StringDT.split(',')[1].split(' ')[0]
343        Srefdate = StringDT.split(',')[1].split(' ')[2]
344
345# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
346##
347        yrref=Srefdate[0:4]
348        monref=Srefdate[5:7]
349        dayref=Srefdate[8:10]
350
351        trefT = Srefdate.find(':')
352        if not trefT == -1:
353#            print '  ' + fname + ': refdate with time!'
354            horref=Srefdate[11:13]
355            minref=Srefdate[14:16]
356            secref=Srefdate[17:19]
357            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
358              '_' + horref + ':' + minref + ':' + secref)
359        else:
360            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
361              + '_00:00:00')
362
363        if tunits == 'weeks':
364            newdate = refdate + dt.timedelta(weeks=float(timeval))
365        elif tunits == 'days':
366            newdate = refdate + dt.timedelta(days=float(timeval))
367        elif tunits == 'hours':
368            newdate = refdate + dt.timedelta(hours=float(timeval))
369        elif tunits == 'minutes':
370            newdate = refdate + dt.timedelta(minutes=float(timeval))
371        elif tunits == 'seconds':
372            newdate = refdate + dt.timedelta(seconds=float(timeval))
373        elif tunits == 'milliseconds':
374            newdate = refdate + dt.timedelta(milliseconds=float(timeval))
375        else:
376              print errormsg
377              print '    timeref_datetime: time units "' + tunits + '" not ready!!!!'
378              quit(-1)
379
380        yr = newdate.year
381        mo = newdate.month
382        da = newdate.day
383        ho = newdate.hour
384        mi = newdate.minute
385        se = newdate.second
386    elif typeSi == 'matYmdHMS':
387        yr = StringDT[0]
388        mo = StringDT[1]
389        da = StringDT[2]
390        ho = StringDT[3]
391        mi = StringDT[4]
392        se = StringDT[5]
393    elif typeSi == 'YmdHMS':
394        yr = int(StringDT[0:4])
395        mo = int(StringDT[4:6])
396        da = int(StringDT[6:8])
397        ho = int(StringDT[8:10])
398        mi = int(StringDT[10:12])
399        se = int(StringDT[12:14])
400    elif typeSi == 'Y-m-d_H:M:S':
401        dateDT = StringDT.split('_')
402        dateD = dateDT[0].split('-')
403        timeT = dateDT[1].split(':')
404        yr = int(dateD[0])
405        mo = int(dateD[1])
406        da = int(dateD[2])
407        ho = int(timeT[0])
408        mi = int(timeT[1])
409        se = int(timeT[2])
410    elif typeSi == 'Y-m-d H:M:S':
411        dateDT = StringDT.split(' ')
412        dateD = dateDT[0].split('-')
413        timeT = dateDT[1].split(':')
414        yr = int(dateD[0])
415        mo = int(dateD[1])
416        da = int(dateD[2])
417        ho = int(timeT[0])
418        mi = int(timeT[1])
419        se = int(timeT[2])
420    elif typeSi == 'Y/m/d H-M-S':
421        dateDT = StringDT.split(' ')
422        dateD = dateDT[0].split('/')
423        timeT = dateDT[1].split('-')
424        yr = int(dateD[0])
425        mo = int(dateD[1])
426        da = int(dateD[2])
427        ho = int(timeT[0])
428        mi = int(timeT[1])
429        se = int(timeT[2])
430    elif typeSi == 'WRFdatetime':
431        yr = int(StringDT[0])*1000 + int(StringDT[1])*100 + int(StringDT[2])*10 +    \
432          int(StringDT[3])
433        mo = int(StringDT[5])*10 + int(StringDT[6])
434        da = int(StringDT[8])*10 + int(StringDT[9])
435        ho = int(StringDT[11])*10 + int(StringDT[12])
436        mi = int(StringDT[14])*10 + int(StringDT[15])
437        se = int(StringDT[17])*10 + int(StringDT[18])
438    else:
439        print errormsg
440        print '  ' + fname + ': type of String input date "' + typeSi +              \
441          '" not ready !!!!'
442        quit(-1)
443
444    if typeSo == 'matYmdHMS':
445        dateYmdHMS = np.zeros((6), dtype=int)
446        dateYmdHMS[0] =  yr
447        dateYmdHMS[1] =  mo
448        dateYmdHMS[2] =  da
449        dateYmdHMS[3] =  ho
450        dateYmdHMS[4] =  mi
451        dateYmdHMS[5] =  se
452    elif typeSo == 'YmdHMS':
453        dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) +        \
454          str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
455    elif typeSo == 'Y-m-d_H:M:S':
456        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
457          str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
458          str(se).zfill(2)
459    elif typeSo == 'Y-m-d H:M:S':
460        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
461          str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
462          str(se).zfill(2)
463    elif typeSo == 'Y/m/d H-M-S':
464        dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' +               \
465          str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
466          str(se).zfill(2)
467    elif typeSo == 'WRFdatetime':
468        dateYmdHMS = []
469        yM = yr/1000
470        yC = (yr-yM*1000)/100
471        yD = (yr-yM*1000-yC*100)/10
472        yU = yr-yM*1000-yC*100-yD*10
473
474        mD = mo/10
475        mU = mo-mD*10
476       
477        dD = da/10
478        dU = da-dD*10
479
480        hD = ho/10
481        hU = ho-hD*10
482
483        miD = mi/10
484        miU = mi-miD*10
485
486        sD = se/10
487        sU = se-sD*10
488
489        dateYmdHMS.append(str(yM))
490        dateYmdHMS.append(str(yC))
491        dateYmdHMS.append(str(yD))
492        dateYmdHMS.append(str(yU))
493        dateYmdHMS.append('-')
494        dateYmdHMS.append(str(mD))
495        dateYmdHMS.append(str(mU))
496        dateYmdHMS.append('-')
497        dateYmdHMS.append(str(dD))
498        dateYmdHMS.append(str(dU))
499        dateYmdHMS.append('_')
500        dateYmdHMS.append(str(hD))
501        dateYmdHMS.append(str(hU))
502        dateYmdHMS.append(':')
503        dateYmdHMS.append(str(miD))
504        dateYmdHMS.append(str(miU))
505        dateYmdHMS.append(':')
506        dateYmdHMS.append(str(sD))
507        dateYmdHMS.append(str(sU))
508    else:
509        print errormsg
510        print '  ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
511        quit(-1)
512
513    return dateYmdHMS
514
515def percendone(nvals,tot,percen,msg):
516    """ Function to provide the percentage of an action across the matrix
517    nvals=number of values
518    tot=total number of values
519    percen=percentage frequency for which the message is wanted
520    msg= message
521    """
522    from sys import stdout
523
524    num = int(tot * percen/100)
525    if (nvals%num == 0): 
526        print '\r        ' + msg + '{0:8.3g}'.format(nvals*100./tot) + ' %',
527        stdout.flush()
528
529    return ''
530
531def netCDFdatetime_realdatetime(units, tcalendar, times):
532    """ Function to transfrom from netCDF CF-compilant times to real time
533    """
534    import datetime as dt
535
536    txtunits = units.split(' ')
537    tunits = txtunits[0]
538    Srefdate = txtunits[len(txtunits) - 1]
539
540# Calendar type
541##
542    is360 = False
543    if tcalendar is not None:
544      print '  netCDFdatetime_realdatetime: There is a calendar attribute'
545      if tcalendar == '365_day' or tcalendar == 'noleap':
546          print '    netCDFdatetime_realdatetime: No leap years!'
547          isleapcal = False
548      elif tcalendar == 'proleptic_gregorian' or tcalendar == 'standard' or tcalendar == 'gregorian':
549          isleapcal = True
550      elif tcalendar == '360_day':
551          is360 = True
552          isleapcal = False
553      else:
554          print errormsg
555          print '    netCDFdatetime_realdatetime: Calendar "' + tcalendar + '" not prepared!'
556          quit(-1)
557
558# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
559##
560    timeval = Srefdate.find(':')
561
562    if not timeval == -1:
563        print '  netCDFdatetime_realdatetime: refdate with time!'
564        refdate = datetimeStr_datetime(Srefdate)
565    else:
566        refdate = dateStr_date(Srefdate + '_00:00:00')
567
568    dimt = len(times)
569#    datetype = type(dt.datetime(1972,02,01))
570#    realdates = np.array(dimt, datetype)
571#    print realdates
572
573## Not in timedelta
574#  if tunits == 'years':
575#    for it in range(dimt):
576#      realdate = refdate + dt.timedelta(years=float(times[it]))
577#      realdates[it] = int(realdate.year)
578#  elif tunits == 'months':
579#    for it in range(dimt):
580#      realdate = refdate + dt.timedelta(months=float(times[it]))
581#      realdates[it] = int(realdate.year)
582#    realdates = []
583    realdates = np.zeros((dimt, 6), dtype=int)
584    if tunits == 'weeks':
585        for it in range(dimt):
586            realdate = refdate + dt.timedelta(weeks=float(times[it]))
587            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
588    elif tunits == 'days':
589        for it in range(dimt):
590            realdate = refdate + dt.timedelta(days=float(times[it]))
591            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
592    elif tunits == 'hours':
593        for it in range(dimt):
594            realdate = refdate + dt.timedelta(hours=float(times[it]))
595#            if not isleapcal:
596#                Nleapdays = cal.leapdays(int(refdate.year), int(realdate.year))
597#                realdate = realdate - dt.timedelta(days=Nleapdays)
598#            if is360:
599#                Nyears360 = int(realdate.year) - int(refdate.year) + 1
600#                realdate = realdate -dt.timedelta(days=Nyears360*5)
601#            realdates[it] = realdate
602#        realdates = refdate + dt.timedelta(hours=float(times))
603            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
604    elif tunits == 'minutes':
605        for it in range(dimt):
606            realdate = refdate + dt.timedelta(minutes=float(times[it]))
607            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
608    elif tunits == 'seconds':
609        for it in range(dimt):
610            realdate = refdate + dt.timedelta(seconds=float(times[it]))
611            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
612    elif tunits == 'milliseconds':
613        for it in range(dimt):
614            realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
615            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
616    elif tunits == 'microseconds':
617        for it in range(dimt):
618            realdate = refdate + dt.timedelta(microseconds=float(times[it]))
619            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
620    else:
621        print errormsg
622        print '  netCDFdatetime_realdatetime: time units "' + tunits + '" is not ready!!!'
623        quit(-1)
624
625    return realdates
626
627def file_nlines(filen):
628    """ Function to provide the number of lines of a file
629    filen= name of the file
630    >>> file_nlines('trajectory.dat')
631    49
632    """
633    fname = 'file_nlines'
634
635    if not os.path.isfile(filen):
636        print errormsg
637        print '  ' + fname + ' file: "' + filen + '" does not exist !!'
638        quit(-1)
639
640    fo = open(filen,'r')
641
642    nlines=0
643    for line in fo: nlines = nlines + 1
644
645    fo.close()
646
647    return nlines
648
649def realdatetime1_CFcompilant(time, Srefdate, tunits):
650    """ Function to transform a matrix with a real time value ([year, month, day,
651      hour, minute, second]) to a netCDF one
652        time= matrix with time
653        Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
654        tunits= units of time respect to Srefdate
655    >>> realdatetime1_CFcompilant([1976, 2, 17, 8, 20, 0], '19491201000000', 'hours')
656    229784.33333333
657    """ 
658
659    import datetime as dt
660    yrref=int(Srefdate[0:4])
661    monref=int(Srefdate[4:6])
662    dayref=int(Srefdate[6:8])
663    horref=int(Srefdate[8:10])
664    minref=int(Srefdate[10:12])
665    secref=int(Srefdate[12:14])
666 
667    refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref)
668
669    if tunits == 'weeks':
670        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5])-refdate
671        cfdates = (cfdate.days + cfdate.seconds/(3600.*24.))/7.
672    elif tunits == 'days':
673        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
674        cfdates = cfdate.days + cfdate.seconds/(3600.*24.)
675    elif tunits == 'hours':
676        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
677        cfdates = cfdate.days*24. + cfdate.seconds/3600.
678    elif tunits == 'minutes':
679        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
680        cfdates = cfdate.days*24.*60. + cfdate.seconds/60.
681    elif tunits == 'seconds':
682        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
683        cfdates = cfdate.days*24.*3600. + cfdate.seconds
684    elif tunits == 'milliseconds':
685        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
686        cfdates = cfdate.days*1000.*24.*3600. + cfdate.seconds*1000.
687    elif tunits == 'microseconds':
688        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],times[5]) - refdate
689        cfdates = cfdate.days*1000000.*24.*3600. + cfdate.seconds*1000000.
690    else:
691        print errormsg
692        print '  ' + fname + ': time units "' + tunits + '" is not ready!!!'
693        quit(-1)
694
695    return cfdates
696
697def basicvardef(varobj, vstname, vlname, vunits):
698    """ Function to give the basic attributes to a variable
699    varobj= netCDF variable object
700    vstname= standard name of the variable
701    vlname= long name of the variable
702    vunits= units of the variable
703    """
704    attr = varobj.setncattr('standard_name', vstname)
705    attr = varobj.setncattr('long_name', vlname)
706    attr = varobj.setncattr('units', vunits)
707
708    return 
709
710def variables_values(varName):
711    """ Function to provide values to plot the different variables values from ASCII file
712      'variables_values.dat'
713    variables_values(varName)
714      [varName]= name of the variable
715        return: [var name], [std name], [minimum], [maximum],
716          [long name]('|' for spaces), [units], [color palette] (following:
717          http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
718     [varn]: original name of the variable
719       NOTE: It might be better doing it with an external ASII file. But then we
720         got an extra dependency...
721    >>> variables_values('WRFght')
722    ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
723    """
724    import subprocess as sub
725
726    fname='variables_values'
727
728    if varName == 'h':
729        print fname + '_____________________________________________________________'
730        print variables_values.__doc__
731        quit()
732
733# This does not work....
734#    folderins = sub.Popen(["pwd"], stdout=sub.PIPE)
735#    folder = list(folderins.communicate())[0].replace('\n','')
736# From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
737    folder = os.path.dirname(os.path.realpath(__file__))
738
739    infile = folder + '/variables_values.dat'
740
741    if not os.path.isfile(infile):
742        print errormsg
743        print '  ' + fname + ": File '" + infile + "' does not exist !!"
744        quit(-1)
745
746# Variable name might come with a statistical surname...
747    stats=['min','max','mean','stdv', 'sum']
748
749# Variables with a statistical section on their name...
750    NOstatsvars = ['zmaxth', 'zmax_th', 'lmax_th', 'lmaxth']
751
752    ifst = False
753    if not searchInlist(NOstatsvars, varName.lower()):
754        for st in stats:
755            if varName.find(st) > -1:
756                print '    '+ fname + ": varibale '" + varName + "' with a " +       \
757                  "statistical surname: '",st,"' !!"
758                Lst = len(st)
759                LvarName = len(varName)
760                varn = varName[0:LvarName - Lst]
761                ifst = True
762                break
763    if not ifst:
764        varn = varName
765
766    ncf = open(infile, 'r')
767
768    for line in ncf:
769        if line[0:1] != '#':
770            values = line.replace('\n','').split(',')
771            if len(values) != 8:
772                print errormsg
773                print "problem in varibale:'", values[0],                            \
774                  'it should have 8 values and it has',len(values)
775                quit(-1)
776
777            if varn[0:6] == 'varDIM': 
778# Variable from a dimension (all with 'varDIM' prefix)
779                Lvarn = len(varn)
780                varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                 \
781                  "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1',  \
782                   'rainbow']
783            else:
784                varvals = [values[1].replace(' ',''), values[2].replace(' ',''),     \
785                  np.float(values[3]), np.float(values[4]),values[5].replace(' ',''),\
786                  values[6].replace(' ',''), values[7].replace(' ','')]
787            if values[0] == varn:
788                ncf.close()
789                return varvals
790                break
791
792    print errormsg
793    print '  ' + fname + ": variable '" + varn + "' not defined !!!"
794    ncf.close()
795    quit(-1)
796
797    return 
798
799def variables_values_old(varName):
800    """ Function to provide values to plot the different variables
801    variables_values(varName)
802      [varName]= name of the variable
803        return: [var name], [std name], [minimum], [maximum],
804          [long name]('|' for spaces), [units], [color palette] (following:
805          http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
806     [varn]: original name of the variable
807       NOTE: It might be better doing it with an external ASII file. But then we
808         got an extra dependency...
809    >>> variables_values('WRFght')
810    ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
811    """
812    fname='variables_values'
813
814    if varName == 'h':
815        print fname + '_____________________________________________________________'
816        print variables_values.__doc__
817        quit()
818
819# Variable name might come with a statistical surname...
820    stats=['min','max','mean','stdv', 'sum']
821
822    ifst = False
823    for st in stats:
824        if varName.find(st) > -1:
825            print '    '+ fname + ": varibale '" + varName + "' with a statistical "+\
826              " surname: '",st,"' !!"
827            Lst = len(st)
828            LvarName = len(varName)
829            varn = varName[0:LvarName - Lst]
830            ifst = True
831            break
832    if not ifst:
833        varn = varName
834
835    if varn[0:6] == 'varDIM': 
836# Variable from a dimension (all with 'varDIM' prefix)
837        Lvarn = len(varn)
838        varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                         \
839          "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1', 'rainbox']
840    elif varn == 'a_tht' or varn == 'LA_THT':
841        varvals = ['ath', 'total_thermal_plume_cover', 0., 1.,                       \
842        'total|column|thermal|plume|cover', '1', 'YlGnBu']
843    elif varn == 'acprc' or varn == 'RAINC':
844        varvals = ['acprc', 'accumulated_cmulus_precipitation', 0., 3.e4,            \
845          'accumulated|cmulus|precipitation', 'mm', 'Blues']
846    elif varn == 'acprnc' or varn == 'RAINNC':
847        varvals = ['acprnc', 'accumulated_non-cmulus_precipitation', 0., 3.e4,       \
848          'accumulated|non-cmulus|precipitation', 'mm', 'Blues']
849    elif varn == 'bils' or varn == 'LBILS':
850        varvals = ['bils', 'surface_total_heat_flux', -100., 100.,                   \
851          'surface|total|heat|flux', 'Wm-2', 'seismic']
852    elif varn == 'landcat' or varn == 'category':
853        varvals = ['landcat', 'land_categories', 0., 22., 'land|categories', '1',    \
854          'rainbow']
855    elif varn == 'c' or varn == 'QCLOUD' or varn == 'oliq' or varn == 'OLIQ':
856        varvals = ['c', 'condensed_water_mixing_ratio', 0., 3.e-4,                   \
857          'condensed|water|mixing|ratio', 'kgkg-1', 'BuPu']
858    elif varn == 'ci' or varn == 'iwcon' or varn == 'LIWCON':
859        varvals = ['ci', 'cloud_iced_water_mixing_ratio', 0., 0.0003,                \
860         'cloud|iced|water|mixing|ratio', 'kgkg-1', 'Purples']
861    elif varn == 'cl' or varn == 'lwcon' or varn == 'LLWCON':
862        varvals = ['cl', 'cloud_liquidwater_mixing_ratio', 0., 0.0003,               \
863         'cloud|liquid|water|mixing|ratio', 'kgkg-1', 'Blues']
864    elif varn == 'cld' or varn == 'CLDFRA' or varn == 'rneb' or varn == 'lrneb' or   \
865      varn == 'LRNEB':
866        varvals = ['cld', 'cloud_area_fraction', 0., 1., 'cloud|fraction', '1',      \
867          'gist_gray']
868    elif varn == 'cldc' or varn == 'rnebcon' or varn == 'lrnebcon' or                \
869      varn == 'LRNEBCON':
870        varvals = ['cldc', 'convective_cloud_area_fraction', 0., 1.,                 \
871          'convective|cloud|fraction', '1', 'gist_gray']
872    elif varn == 'cldl' or varn == 'rnebls' or varn == 'lrnebls' or varn == 'LRNEBLS':
873        varvals = ['cldl', 'large_scale_cloud_area_fraction', 0., 1.,                \
874          'large|scale|cloud|fraction', '1', 'gist_gray']
875    elif varn == 'clt' or varn == 'CLT' or varn == 'cldt' or                         \
876      varn == 'Total cloudiness':
877        varvals = ['clt', 'cloud_area_fraction', 0., 1., 'total|cloud|cover', '1',   \
878          'gist_gray']
879    elif varn == 'cll' or varn == 'cldl' or varn == 'LCLDL' or                       \
880      varn == 'Low-level cloudiness':
881        varvals = ['cll', 'low_level_cloud_area_fraction', 0., 1.,                   \
882          'low|level|(p|>|680|hPa)|cloud|fraction', '1', 'gist_gray']
883    elif varn == 'clm' or varn == 'cldm' or varn == 'LCLDM' or                       \
884      varn == 'Mid-level cloudiness':
885        varvals = ['clm', 'mid_level_cloud_area_fraction', 0., 1.,                   \
886          'medium|level|(440|<|p|<|680|hPa)|cloud|fraction', '1', 'gist_gray']
887    elif varn == 'clh' or varn == 'cldh' or varn == 'LCLDH' or                       \
888      varn == 'High-level cloudiness':
889        varvals = ['clh', 'high_level_cloud_area_fraction', 0., 1.,                  \
890          'high|level|(p|<|440|hPa)|cloud|fraction', '1', 'gist_gray']
891    elif varn == 'clmf' or varn == 'fbase' or varn == 'LFBASE':
892        varvals = ['clmf', 'cloud_base_max_flux', -0.3, 0.3, 'cloud|base|max|flux',  \
893          'kgm-2s-1', 'seismic']
894    elif varn == 'clp' or varn == 'pbase' or varn == 'LPBASE':
895        varvals = ['clp', 'cloud_base_pressure', -0.3, 0.3, 'cloud|base|pressure',   \
896          'Pa', 'Reds']
897    elif varn == 'cpt' or varn == 'ptconv' or varn == 'LPTCONV':
898        varvals = ['cpt', 'convective_point', 0., 1., 'convective|point', '1',       \
899          'seismic']
900    elif varn == 'dqajs' or varn == 'LDQAJS':
901        varvals = ['dqajs', 'dry_adjustment_water_vapor_tendency', -0.0003, 0.0003,  \
902        'dry|adjustment|water|vapor|tendency', 'kg/kg/s', 'seismic']
903    elif varn == 'dqcon' or varn == 'LDQCON':
904        varvals = ['dqcon', 'convective_water_vapor_tendency', -3e-8, 3.e-8,         \
905        'convective|water|vapor|tendency', 'kg/kg/s', 'seismic']
906    elif varn == 'dqdyn' or varn == 'LDQDYN':
907        varvals = ['dqdyn', 'dynamics_water_vapor_tendency', -3.e-7, 3.e-7,          \
908        'dynamics|water|vapor|tendency', 'kg/kg/s', 'seismic']
909    elif varn == 'dqeva' or varn == 'LDQEVA':
910        varvals = ['dqeva', 'evaporation_water_vapor_tendency', -3.e-6, 3.e-6,       \
911        'evaporation|water|vapor|tendency', 'kg/kg/s', 'seismic']
912    elif varn == 'dqlscst' or varn == 'LDQLSCST':
913        varvals = ['dqlscst', 'stratocumulus_water_vapor_tendency', -3.e-7, 3.e-7,   \
914        'stratocumulus|water|vapor|tendency', 'kg/kg/s', 'seismic']
915    elif varn == 'dqlscth' or varn == 'LDQLSCTH': 
916        varvals = ['dqlscth', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7,        \
917        'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic']
918    elif varn == 'dqlsc' or varn == 'LDQLSC':
919        varvals = ['dqlsc', 'condensation_water_vapor_tendency', -3.e-6, 3.e-6,      \
920        'condensation|water|vapor|tendency', 'kg/kg/s', 'seismic']
921    elif varn == 'dqphy' or varn == 'LDQPHY':
922        varvals = ['dqphy', 'physics_water_vapor_tendency', -3.e-7, 3.e-7,           \
923        'physics|water|vapor|tendency', 'kg/kg/s', 'seismic']
924    elif varn == 'dqthe' or varn == 'LDQTHE':
925        varvals = ['dqthe', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7,          \
926        'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic']
927    elif varn == 'dqvdf' or varn == 'LDQVDF':
928        varvals = ['dqvdf', 'vertical_difussion_water_vapor_tendency', -3.e-8, 3.e-8,\
929        'vertical|difussion|water|vapor|tendency', 'kg/kg/s', 'seismic']
930    elif varn == 'dqwak' or varn == 'LDQWAK':
931        varvals = ['dqwak', 'wake_water_vapor_tendency', -3.e-7, 3.e-7,              \
932        'wake|water|vapor|tendency', 'kg/kg/s', 'seismic']
933    elif varn == 'dta' or varn == 'tnt' or varn == 'LTNT':
934        varvals = ['dta', 'tendency_air_temperature', -3.e-3, 3.e-3,                 \
935        'tendency|of|air|temperature', 'K/s', 'seismic']
936    elif varn == 'dtac' or varn == 'tntc' or varn == 'LTNTC':
937        varvals = ['dtac', 'moist_convection_tendency_air_temperature', -3.e-3,      \
938        3.e-3, 'moist|convection|tendency|of|air|temperature', 'K/s', 'seismic']
939    elif varn == 'dtar' or varn == 'tntr' or varn == 'LTNTR':
940        varvals = ['dtar', 'radiative_heating_tendency_air_temperature', -3.e-3,     \
941          3.e-3, 'radiative|heating|tendency|of|air|temperature', 'K/s', 'seismic']
942    elif varn == 'dtascpbl' or varn == 'tntscpbl' or varn == 'LTNTSCPBL':
943        varvals = ['dtascpbl',                                                       \
944          'stratiform_cloud_precipitation_BL_mixing_tendency_air_temperature',       \
945          -3.e-6, 3.e-6,                                                             \
946          'stratiform|cloud|precipitation|Boundary|Layer|mixing|tendency|air|'       +
947          'temperature', 'K/s', 'seismic']
948    elif varn == 'dtajs' or varn == 'LDTAJS':
949        varvals = ['dtajs', 'dry_adjustment_thermal_tendency', -3.e-5, 3.e-5,        \
950        'dry|adjustment|thermal|tendency', 'K/s', 'seismic']
951    elif varn == 'dtcon' or varn == 'LDTCON':
952        varvals = ['dtcon', 'convective_thermal_tendency', -3.e-5, 3.e-5,            \
953        'convective|thermal|tendency', 'K/s', 'seismic']
954    elif varn == 'dtdyn' or varn == 'LDTDYN':
955        varvals = ['dtdyn', 'dynamics_thermal_tendency', -3.e-4, 3.e-4,              \
956        'dynamics|thermal|tendency', 'K/s', 'seismic']
957    elif varn == 'dteva' or varn == 'LDTEVA':
958        varvals = ['dteva', 'evaporation_thermal_tendency', -3.e-3, 3.e-3,           \
959        'evaporation|thermal|tendency', 'K/s', 'seismic']
960    elif varn == 'dtlscst' or varn == 'LDTLSCST':
961        varvals = ['dtlscst', 'stratocumulus_thermal_tendency', -3.e-4, 3.e-4,       \
962        'stratocumulus|thermal|tendency', 'K/s', 'seismic']
963    elif varn == 'dtlscth' or varn == 'LDTLSCTH':
964        varvals = ['dtlscth', 'thermals_thermal_tendency', -3.e-4, 3.e-4,            \
965        'thermal|plumes|thermal|tendency', 'K/s', 'seismic']
966    elif varn == 'dtlsc' or varn == 'LDTLSC':
967        varvals = ['dtlsc', 'condensation_thermal_tendency', -3.e-3, 3.e-3,          \
968        'condensation|thermal|tendency', 'K/s', 'seismic']
969    elif varn == 'dtlwr' or varn == 'LDTLWR':
970        varvals = ['dtlwr', 'long_wave_thermal_tendency', -3.e-3, 3.e-3, \
971        'long|wave|radiation|thermal|tendency', 'K/s', 'seismic']
972    elif varn == 'dtphy' or varn == 'LDTPHY':
973        varvals = ['dtphy', 'physics_thermal_tendency', -3.e-4, 3.e-4,               \
974        'physics|thermal|tendency', 'K/s', 'seismic']
975    elif varn == 'dtsw0' or varn == 'LDTSW0':
976        varvals = ['dtsw0', 'cloudy_sky_short_wave_thermal_tendency', -3.e-4, 3.e-4, \
977        'cloudy|sky|short|wave|radiation|thermal|tendency', 'K/s', 'seismic']
978    elif varn == 'dtthe' or varn == 'LDTTHE':
979        varvals = ['dtthe', 'thermals_thermal_tendency', -3.e-4, 3.e-4,              \
980        'thermal|plumes|thermal|tendency', 'K/s', 'seismic']
981    elif varn == 'dtvdf' or varn == 'LDTVDF':
982        varvals = ['dtvdf', 'vertical_difussion_thermal_tendency', -3.e-5, 3.e-5,    \
983        'vertical|difussion|thermal|tendency', 'K/s', 'seismic']
984    elif varn == 'dtwak' or varn == 'LDTWAK':
985        varvals = ['dtwak', 'wake_thermal_tendency', -3.e-4, 3.e-4,                  \
986        'wake|thermal|tendency', 'K/s', 'seismic']
987    elif varn == 'ducon' or varn == 'LDUCON':
988        varvals = ['ducon', 'convective_eastward_wind_tendency', -3.e-3, 3.e-3,      \
989        'convective|eastward|wind|tendency', 'ms-2', 'seismic']
990    elif varn == 'dudyn' or varn == 'LDUDYN':
991        varvals = ['dudyn', 'dynamics_eastward_wind_tendency', -3.e-3, 3.e-3,        \
992        'dynamics|eastward|wind|tendency', 'ms-2', 'seismic']
993    elif varn == 'duvdf' or varn == 'LDUVDF':
994        varvals = ['duvdf', 'vertical_difussion_eastward_wind_tendency', -3.e-3,     \
995         3.e-3, 'vertical|difussion|eastward|wind|tendency', 'ms-2', 'seismic']
996    elif varn == 'dvcon' or varn == 'LDVCON':
997        varvals = ['dvcon', 'convective_difussion_northward_wind_tendency', -3.e-3,  \
998         3.e-3, 'convective|northward|wind|tendency', 'ms-2', 'seismic']
999    elif varn == 'dvdyn' or varn == 'LDVDYN':
1000        varvals = ['dvdyn', 'dynamics_northward_wind_tendency', -3.e-3,              \
1001         3.e-3, 'dynamics|difussion|northward|wind|tendency', 'ms-2', 'seismic']
1002    elif varn == 'dvvdf' or varn == 'LDVVDF':
1003        varvals = ['dvvdf', 'vertical_difussion_northward_wind_tendency', -3.e-3,    \
1004         3.e-3, 'vertical|difussion|northward|wind|tendency', 'ms-2', 'seismic']
1005    elif varn == 'etau' or varn == 'ZNU':
1006        varvals = ['etau', 'etau', 0., 1, 'eta values on half (mass) levels', '-',   \
1007        'reds']
1008    elif varn == 'evspsbl' or varn == 'LEVAP' or varn == 'evap' or varn == 'SFCEVPde':
1009        varvals = ['evspsbl', 'water_evaporation_flux', 0., 1.5e-4,                  \
1010          'water|evaporation|flux', 'kgm-2s-1', 'Blues']
1011    elif varn == 'evspsbl' or varn == 'SFCEVPde':
1012        varvals = ['evspsblac', 'water_evaporation_flux_ac', 0., 1.5e-4,             \
1013          'accumulated|water|evaporation|flux', 'kgm-2', 'Blues']
1014    elif varn == 'g' or varn == 'QGRAUPEL':
1015        varvals = ['g', 'grauepl_mixing_ratio', 0., 0.0003, 'graupel|mixing|ratio',  \
1016          'kgkg-1', 'Purples']
1017    elif varn == 'h2o' or varn == 'LH2O':
1018        varvals = ['h2o', 'water_mass_fraction', 0., 3.e-2,                          \
1019          'mass|fraction|of|water', '1', 'Blues']
1020    elif varn == 'h' or varn == 'QHAIL':
1021        varvals = ['h', 'hail_mixing_ratio', 0., 0.0003, 'hail|mixing|ratio',        \
1022          'kgkg-1', 'Purples']
1023    elif varn == 'hfls' or varn == 'LH' or varn == 'LFLAT' or varn == 'flat':
1024        varvals = ['hfls', 'surface_upward_latent_heat_flux', -400., 400.,           \
1025          'upward|latnt|heat|flux|at|the|surface', 'Wm-2', 'seismic']
1026    elif varn == 'hfss' or varn == 'LSENS' or varn == 'sens' or varn == 'HFX':
1027        varvals = ['hfss', 'surface_upward_sensible_heat_flux', -150., 150.,         \
1028          'upward|sensible|heat|flux|at|the|surface', 'Wm-2', 'seismic']
1029    elif varn == 'hfso' or varn == 'GRDFLX':
1030        varvals = ['hfso', 'downward_heat_flux_in_soil', -150., 150.,                \
1031          'Downward|soil|heat|flux', 'Wm-2', 'seismic']
1032    elif varn == 'hus' or varn == 'WRFrh' or varn == 'LMDZrh' or varn == 'rhum' or   \
1033      varn == 'LRHUM':
1034        varvals = ['hus', 'specific_humidity', 0., 1., 'specific|humidty', '1',      \
1035          'BuPu']
1036    elif varn == 'huss' or varn == 'WRFrhs' or varn == 'LMDZrhs' or varn == 'rh2m' or\
1037      varn == 'LRH2M':
1038        varvals = ['huss', 'specific_humidity', 0., 1., 'specific|humidty|at|2m',    \
1039          '1', 'BuPu']
1040    elif varn == 'i' or varn == 'QICE':
1041        varvals = ['i', 'iced_water_mixing_ratio', 0., 0.0003,                       \
1042         'iced|water|mixing|ratio', 'kgkg-1', 'Purples']
1043    elif varn == 'lat' or varn == 'XLAT' or varn == 'XLAT_M' or varn == 'latitude':
1044        varvals = ['lat', 'latitude', -90., 90., 'latitude', 'degrees North',        \
1045          'seismic']
1046    elif varn == 'lcl' or varn == 's_lcl' or varn == 'ls_lcl' or varn == 'LS_LCL':
1047        varvals = ['lcl', 'condensation_level', 0., 2500., 'level|of|condensation',  \
1048          'm', 'Greens']
1049    elif varn == 'lambdath' or varn == 'lambda_th' or varn == 'LLAMBDA_TH':
1050        varvals = ['lambdath', 'thermal_plume_vertical_velocity', -30., 30.,         \
1051          'thermal|plume|vertical|velocity', 'm/s', 'seismic']
1052    elif varn == 'lmaxth' or varn == 'LLMAXTH':
1053        varvals = ['lmaxth', 'upper_level_thermals', 0., 100., 'upper|level|thermals'\
1054          , '1', 'Greens']
1055    elif varn == 'lon' or varn == 'XLONG' or varn == 'XLONG_M':
1056        varvals = ['lon', 'longitude', -180., 180., 'longitude', 'degrees East',     \
1057          'seismic']
1058    elif varn == 'longitude':
1059        varvals = ['lon', 'longitude', 0., 360., 'longitude', 'degrees East',        \
1060          'seismic']
1061    elif varn == 'orog' or varn == 'HGT' or varn == 'HGT_M':
1062        varvals = ['orog', 'orography',  0., 3000., 'surface|altitude', 'm','terrain']
1063    elif varn == 'pfc' or varn == 'plfc' or varn == 'LPLFC':
1064        varvals = ['pfc', 'pressure_free_convection', 100., 1100.,                   \
1065          'pressure|free|convection', 'hPa', 'BuPu']
1066    elif varn == 'plcl' or varn == 'LPLCL':
1067        varvals = ['plcl', 'pressure_lifting_condensation_level', 700., 1100.,       \
1068          'pressure|lifting|condensation|level', 'hPa', 'BuPu']
1069    elif varn == 'pr' or varn == 'RAINTOT' or varn == 'precip' or                    \
1070      varn == 'LPRECIP' or varn == 'Precip Totale liq+sol':
1071        varvals = ['pr', 'precipitation_flux', 0., 1.e-4, 'precipitation|flux',      \
1072          'kgm-2s-1', 'BuPu']
1073    elif varn == 'prprof' or varn == 'vprecip' or varn == 'LVPRECIP':
1074        varvals = ['prprof', 'precipitation_profile', 0., 1.e-3,                     \
1075          'precipitation|profile', 'kg/m2/s', 'BuPu']
1076    elif varn == 'prprofci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
1077        varvals = ['prprofci', 'precipitation_profile_convective_i', 0., 1.e-3,      \
1078          'precipitation|profile|convective|i', 'kg/m2/s', 'BuPu']
1079    elif varn == 'prprofcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
1080        varvals = ['prprofcl', 'precipitation_profile_convective_l', 0., 1.e-3,      \
1081          'precipitation|profile|convective|l', 'kg/m2/s', 'BuPu']
1082    elif varn == 'prprofli' or varn == 'pr_lsc_i' or varn == 'LPR_LSC_I':
1083        varvals = ['prprofli', 'precipitation_profile_large_scale_i', 0., 1.e-3,     \
1084          'precipitation|profile|large|scale|i', 'kg/m2/s', 'BuPu']
1085    elif varn == 'prprofll' or varn == 'pr_lsc_l' or varn == 'LPR_LSC_L':
1086        varvals = ['prprofll', 'precipitation_profile_large_scale_l', 0., 1.e-3,     \
1087          'precipitation|profile|large|scale|l', 'kg/m2/s', 'BuPu']
1088    elif varn == 'pracc' or varn == 'ACRAINTOT':
1089        varvals = ['pracc', 'precipitation_amount', 0., 100.,                        \
1090          'accumulated|precipitation', 'kgm-2', 'BuPu']
1091    elif varn == 'prc' or varn == 'LPLUC' or varn == 'pluc' or varn == 'WRFprc' or   \
1092      varn == 'RAINCde':
1093        varvals = ['prc', 'convective_precipitation_flux', 0., 2.e-4,                \
1094          'convective|precipitation|flux', 'kgm-2s-1', 'Blues']
1095    elif varn == 'prci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
1096        varvals = ['prci', 'convective_ice_precipitation_flux', 0., 0.003,           \
1097          'convective|ice|precipitation|flux', 'kgm-2s-1', 'Purples']
1098    elif varn == 'prcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
1099        varvals = ['prcl', 'convective_liquid_precipitation_flux', 0., 0.003,        \
1100          'convective|liquid|precipitation|flux', 'kgm-2s-1', 'Blues']
1101    elif varn == 'pres' or varn == 'presnivs' or varn == 'pressure' or               \
1102      varn == 'lpres' or varn == 'LPRES':
1103        varvals = ['pres', 'air_pressure', 0., 103000., 'air|pressure', 'Pa',        \
1104          'Blues']
1105    elif varn == 'prls' or varn == 'WRFprls' or varn == 'LPLUL' or varn == 'plul' or \
1106       varn == 'RAINNCde':
1107        varvals = ['prls', 'large_scale_precipitation_flux', 0., 2.e-4,              \
1108          'large|scale|precipitation|flux', 'kgm-2s-1', 'Blues']
1109    elif varn == 'prsn' or varn == 'SNOW' or varn == 'snow' or varn == 'LSNOW':
1110        varvals = ['prsn', 'snowfall', 0., 1.e-4, 'snowfall|flux', 'kgm-2s-1', 'BuPu']
1111    elif varn == 'prw' or varn == 'WRFprh':
1112        varvals = ['prw', 'atmosphere_water_vapor_content', 0., 10.,                 \
1113          'water|vapor"path', 'kgm-2', 'Blues']
1114    elif varn == 'ps' or varn == 'psfc' or varn =='PSFC' or varn == 'psol' or        \
1115      varn == 'Surface Pressure':
1116        varvals=['ps', 'surface_air_pressure', 85000., 105400., 'surface|pressure',  \
1117          'hPa', 'cool']
1118    elif varn == 'psl' or varn == 'mslp' or varn =='WRFmslp':
1119        varvals=['psl', 'air_pressure_at_sea_level', 85000., 104000.,                \
1120          'mean|sea|level|pressure', 'Pa', 'Greens']
1121    elif varn == 'qth' or varn == 'q_th' or varn == 'LQ_TH':
1122        varvals = ['qth', 'thermal_plume_total_water_content', 0., 25.,              \
1123          'total|water|cotent|in|thermal|plume', 'mm', 'YlOrRd']
1124    elif varn == 'r' or varn == 'QVAPOR' or varn == 'ovap' or varn == 'LOVAP':
1125        varvals = ['r', 'water_mixing_ratio', 0., 0.03, 'water|mixing|ratio',        \
1126          'kgkg-1', 'BuPu']
1127    elif varn == 'r2' or varn == 'Q2':
1128        varvals = ['r2', 'water_mixing_ratio_at_2m', 0., 0.03, 'water|mixing|' +     \
1129          'ratio|at|2|m','kgkg-1', 'BuPu']
1130    elif varn == 'rsds' or varn == 'SWdnSFC' or varn == 'SWdn at surface' or         \
1131      varn == 'SWDOWN':
1132        varvals=['rsds', 'surface_downwelling_shortwave_flux_in_air',  0., 1200.,    \
1133          'downward|SW|surface|radiation', 'Wm-2' ,'Reds']
1134    elif varn == 'rsdsacc':
1135        varvals=['rsdsacc', 'accumulated_surface_downwelling_shortwave_flux_in_air', \
1136          0., 1200., 'accumulated|downward|SW|surface|radiation', 'Wm-2' ,'Reds']
1137    elif varn == 'rvor' or varn == 'WRFrvor':
1138        varvals = ['rvor', 'air_relative_vorticity', -2.5E-3, 2.5E-3,                \
1139          'air|relative|vorticity', 's-1', 'seismic']
1140    elif varn == 'rvors' or varn == 'WRFrvors':
1141        varvals = ['rvors', 'surface_air_relative_vorticity', -2.5E-3, 2.5E-3,       \
1142          'surface|air|relative|vorticity', 's-1', 'seismic']
1143    elif varn == 's' or varn == 'QSNOW':
1144        varvals = ['s', 'snow_mixing_ratio', 0., 0.0003, 'snow|mixing|ratio',        \
1145          'kgkg-1', 'Purples']
1146    elif varn == 'stherm' or varn == 'LS_THERM':
1147        varvals = ['stherm', 'thermals_excess', 0., 0.8, 'thermals|excess', 'K',     \
1148          'Reds']
1149    elif varn == 'ta' or varn == 'WRFt' or varn == 'temp' or varn == 'LTEMP' or      \
1150      varn == 'Air temperature':
1151        varvals = ['ta', 'air_temperature', 195., 320., 'air|temperature', 'K',      \
1152          'YlOrRd']
1153    elif varn == 'tah' or varn == 'theta' or varn == 'LTHETA':
1154        varvals = ['tah', 'potential_air_temperature', 195., 320.,                   \
1155          'potential|air|temperature', 'K', 'YlOrRd']
1156    elif varn == 'tas' or varn == 'T2' or varn == 't2m' or varn == 'T2M' or          \
1157      varn == 'Temperature 2m':
1158        varvals = ['tas', 'air_temperature', 240., 310., 'air|temperature|at|2m', \
1159          K', 'YlOrRd']
1160    elif varn == 'tds' or varn == 'TH2':
1161        varvals = ['tds', 'air_dew_point_temperature', 240., 310.,                   \
1162          'air|dew|point|temperature|at|2m', 'K', 'YlGnBu']
1163    elif varn == 'tke' or varn == 'TKE' or varn == 'tke' or varn == 'LTKE':
1164        varvals = ['tke', 'turbulent_kinetic_energy', 0., 0.003,                     \
1165          'turbulent|kinetic|energy', 'm2/s2', 'Reds']
1166    elif varn == 'time'or varn == 'time_counter':
1167        varvals = ['time', 'time', 0., 1000., 'time',                                \
1168          'hours|since|1949/12/01|00:00:00', 'Reds']
1169    elif varn == 'tmla' or varn == 's_pblt' or varn == 'LS_PBLT':
1170        varvals = ['tmla', 'atmosphere_top_boundary_layer_temperature', 250., 330.,  \
1171          'atmosphere|top|boundary|layer|temperature', 'K', 'Reds']
1172    elif varn == 'ua' or varn == 'vitu' or varn == 'U' or varn == 'Zonal wind' or    \
1173      varn == 'LVITU':
1174        varvals = ['ua', 'eastward_wind', -30., 30., 'eastward|wind', 'ms-1',        \
1175          'seismic']
1176    elif varn == 'uas' or varn == 'u10m' or varn == 'U10' or varn =='Vent zonal 10m':
1177        varvals = ['uas', 'eastward_wind', -30., 30., 'eastward|2m|wind',    \
1178          'ms-1', 'seismic']
1179    elif varn == 'va' or varn == 'vitv' or varn == 'V' or varn == 'Meridional wind'  \
1180      or varn == 'LVITV':
1181        varvals = ['va', 'northward_wind', -30., 30., 'northward|wind', 'ms-1',      \
1182          'seismic']
1183    elif varn == 'vas' or varn == 'v10m' or varn == 'V10' or                         \
1184      varn =='Vent meridien 10m':
1185        varvals = ['vas', 'northward_wind', -30., 30., 'northward|2m|wind', 'ms-1',  \
1186          'seismic']
1187    elif varn == 'wakedeltaq' or varn == 'wake_deltaq' or varn == 'lwake_deltaq' or  \
1188      varn == 'LWAKE_DELTAQ':
1189        varvals = ['wakedeltaq', 'wake_delta_vapor', -0.003, 0.003,                  \
1190          'wake|delta|mixing|ratio', '-', 'seismic']
1191    elif varn == 'wakedeltat' or varn == 'wake_deltat' or varn == 'lwake_deltat' or  \
1192      varn == 'LWAKE_DELTAT':
1193        varvals = ['wakedeltat', 'wake_delta_temp', -0.003, 0.003,                   \
1194          'wake|delta|temperature', '-', 'seismic']
1195    elif varn == 'wakeh' or varn == 'wake_h' or varn == 'LWAKE_H':
1196        varvals = ['wakeh', 'wake_height', 0., 1000., 'height|of|the|wakes', 'm',    \
1197          'YlOrRd']
1198    elif varn == 'wakeomg' or varn == 'wake_omg' or varn == 'lwake_omg' or           \
1199      varn == 'LWAKE_OMG':
1200        varvals = ['wakeomg', 'wake_omega', 0., 3., 'wake|omega', \
1201          '-', 'BuGn']
1202    elif varn == 'wakes' or varn == 'wake_s' or varn == 'LWAKE_S':
1203        varvals = ['wakes', 'wake_area_fraction', 0., 0.5, 'wake|spatial|fraction',  \
1204          '1', 'BuGn']
1205    elif varn == 'wa' or varn == 'W' or varn == 'Vertical wind':
1206        varvals = ['wa', 'upward_wind', -10., 10., 'upward|wind', 'ms-1',            \
1207          'seismic']
1208    elif varn == 'wap' or varn == 'vitw' or varn == 'LVITW':
1209        varvals = ['wap', 'upward_wind', -3.e-10, 3.e-10, 'upward|wind', 'mPa-1',    \
1210          'seismic']
1211    elif varn == 'wss' or varn == 'SPDUV':
1212        varvals = ['wss', 'air_velocity',  0., 30., 'surface|horizontal|wind|speed', \
1213          'ms-1', 'Reds']
1214# Water budget
1215# Water budget de-accumulated
1216    elif varn == 'ccond' or varn == 'CCOND' or varn == 'ACCCONDde':
1217        varvals = ['ccond', 'cw_cond',  0., 30.,                                     \
1218          'cloud|water|condensation', 'mm', 'Reds']
1219    elif varn == 'wbr' or varn == 'ACQVAPORde':
1220        varvals = ['wbr', 'wbr',  0., 30., 'Water|Budget|water|wapor', 'mm', 'Blues']
1221    elif varn == 'diabh' or varn == 'DIABH' or varn == 'ACDIABHde':
1222        varvals = ['diabh', 'diabh',  0., 30., 'diabatic|heating', 'K', 'Reds']
1223    elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde':
1224        varvals = ['wbpw', 'water_budget_pw',  0., 30., 'Water|Budget|water|content',\
1225           'mms-1', 'Reds']
1226    elif varn == 'wbf' or varn == 'WBACF' or varn == 'WBACFde':
1227        varvals = ['wbf', 'water_budget_hfcqv',  0., 30.,                       \
1228          'Water|Budget|horizontal|convergence|of|water|vapour|(+,|' +   \
1229          'conv.;|-,|div.)', 'mms-1', 'Reds']
1230    elif varn == 'wbfc' or varn == 'WBFC' or varn == 'WBACFCde':
1231        varvals = ['wbfc', 'water_budget_fc',  0., 30.,                         \
1232          'Water|Budget|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\
1233          'div.)', 'mms-1', 'Reds']
1234    elif varn == 'wbfp' or varn == 'WBFP' or varn == 'WBACFPde':
1235        varvals = ['wbfp', 'water_budget_cfp',  0., 30.,                        \
1236          'Water|Budget|horizontal|convergence|of|precipitation|(+,|' +  \
1237          'conv.;|-,|div.)', 'mms-1', 'Reds']
1238    elif varn == 'wbz' or varn == 'WBZ' or varn == 'WBACZde':
1239        varvals = ['wbz', 'water_budget_z',  0., 30.,                           \
1240          'Water|Budget|vertical|convergence|of|water|vapour|(+,|conv.' +\
1241          ';|-,|div.)', 'mms-1', 'Reds']
1242    elif varn == 'wbc' or varn == 'WBC' or varn == 'WBACCde':
1243        varvals = ['wbc', 'water_budget_c',  0., 30.,                           \
1244          'Water|Budget|Cloud|water|species','mms-1', 'Reds']
1245    elif varn == 'wbqvd' or varn == 'WBQVD' or varn == 'WBACQVDde':
1246        varvals = ['wbqvd', 'water_budget_qvd',  0., 30.,                       \
1247          'Water|Budget|water|vapour|divergence', 'mms-1', 'Reds']
1248    elif varn == 'wbqvblten' or varn == 'WBQVBLTEN' or varn == 'WBACQVBLTENde':
1249        varvals = ['wbqvblten', 'water_budget_qv_blten',  0., 30.,              \
1250          'Water|Budget|QV|tendency|due|to|pbl|parameterization',        \
1251          'kg kg-1 s-1', 'Reds']
1252    elif varn == 'wbqvcuten' or varn == 'WBQVCUTEN' or varn == 'WBACQVCUTENde':
1253        varvals = ['wbqvcuten', 'water_budget_qv_cuten',  0., 30.,              \
1254          'Water|Budget|QV|tendency|due|to|cu|parameterization',         \
1255          'kg kg-1 s-1', 'Reds']
1256    elif varn == 'wbqvshten' or varn == 'WBQVSHTEN' or varn == 'WBACQVSHTENde':
1257        varvals = ['wbqvshten', 'water_budget_qv_shten',  0., 30.,              \
1258          'Water|Budget|QV|tendency|due|to|shallow|cu|parameterization', \
1259          'kg kg-1 s-1', 'Reds']
1260    elif varn == 'wbpr' or varn == 'WBP' or varn == 'WBACPde':
1261        varvals = ['wbpr', 'water_budget_pr',  0., 30.,                         \
1262          'Water|Budget|recipitation', 'mms-1', 'Reds']
1263    elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde':
1264        varvals = ['wbpw', 'water_budget_pw',  0., 30.,                         \
1265          'Water|Budget|water|content', 'mms-1', 'Reds']
1266    elif varn == 'wbcondt' or varn == 'WBCONDT' or varn == 'WBACCONDTde':
1267        varvals = ['wbcondt', 'water_budget_condt',  0., 30.,                   \
1268          'Water|Budget|condensation|and|deposition', 'mms-1', 'Reds']
1269    elif varn == 'wbqcm' or varn == 'WBQCM' or varn == 'WBACQCMde':
1270        varvals = ['wbqcm', 'water_budget_qcm',  0., 30.,                       \
1271          'Water|Budget|hydrometeor|change|and|convergence', 'mms-1', 'Reds']
1272    elif varn == 'wbsi' or varn == 'WBSI' or varn == 'WBACSIde':
1273        varvals = ['wbsi', 'water_budget_si',  0., 30.,                         \
1274          'Water|Budget|hydrometeor|sink', 'mms-1', 'Reds']
1275    elif varn == 'wbso' or varn == 'WBSO' or varn == 'WBACSOde':
1276        varvals = ['wbso', 'water_budget_so',  0., 30.,                         \
1277          'Water|Budget|hydrometeor|source', 'mms-1', 'Reds']
1278# Water Budget accumulated
1279    elif varn == 'ccondac' or varn == 'ACCCOND':
1280        varvals = ['ccondac', 'cw_cond_ac',  0., 30.,                                \
1281          'accumulated|cloud|water|condensation', 'mm', 'Reds']
1282    elif varn == 'rac' or varn == 'ACQVAPOR':
1283        varvals = ['rac', 'ac_r',  0., 30., 'accumualted|water|wapor', 'mm', 'Blues']
1284    elif varn == 'diabhac' or varn == 'ACDIABH':
1285        varvals = ['diabhac', 'diabh_ac',  0., 30., 'accumualted|diabatic|heating',  \
1286          'K', 'Reds']
1287    elif varn == 'wbpwac' or varn == 'WBACPW':
1288        varvals = ['wbpwac', 'water_budget_pw_ac',  0., 30.,                         \
1289          'Water|Budget|accumulated|water|content', 'mm', 'Reds']
1290    elif varn == 'wbfac' or varn == 'WBACF':
1291        varvals = ['wbfac', 'water_budget_hfcqv_ac',  0., 30.,                       \
1292          'Water|Budget|accumulated|horizontal|convergence|of|water|vapour|(+,|' +   \
1293          'conv.;|-,|div.)', 'mm', 'Reds']
1294    elif varn == 'wbfcac' or varn == 'WBACFC':
1295        varvals = ['wbfcac', 'water_budget_fc_ac',  0., 30.,                         \
1296          'Water|Budget|accumulated|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\
1297          'div.)', 'mm', 'Reds']
1298    elif varn == 'wbfpac' or varn == 'WBACFP':
1299        varvals = ['wbfpac', 'water_budget_cfp_ac',  0., 30.,                        \
1300          'Water|Budget|accumulated|horizontal|convergence|of|precipitation|(+,|' +  \
1301          'conv.;|-,|div.)', 'mm', 'Reds']
1302    elif varn == 'wbzac' or varn == 'WBACZ':
1303        varvals = ['wbzac', 'water_budget_z_ac',  0., 30.,                           \
1304          'Water|Budget|accumulated|vertical|convergence|of|water|vapour|(+,|conv.' +\
1305          ';|-,|div.)', 'mm', 'Reds']
1306    elif varn == 'wbcac' or varn == 'WBACC':
1307        varvals = ['wbcac', 'water_budget_c_ac',  0., 30.,                           \
1308          'Water|Budget|accumulated|Cloud|water|species','mm', 'Reds']
1309    elif varn == 'wbqvdac' or varn == 'WBACQVD':
1310        varvals = ['wbqvdac', 'water_budget_qvd_ac',  0., 30.,                       \
1311          'Water|Budget|accumulated|water|vapour|divergence', 'mm', 'Reds']
1312    elif varn == 'wbqvbltenac' or varn == 'WBACQVBLTEN':
1313        varvals = ['wbqvbltenac', 'water_budget_qv_blten_ac',  0., 30.,              \
1314          'Water|Budget|accumulated|QV|tendency|due|to|pbl|parameterization',        \
1315          'kg kg-1 s-1', 'Reds']
1316    elif varn == 'wbqvcutenac' or varn == 'WBACQVCUTEN':
1317        varvals = ['wbqvcutenac', 'water_budget_qv_cuten_ac',  0., 30.,              \
1318          'Water|Budget|accumulated|QV|tendency|due|to|cu|parameterization',         \
1319          'kg kg-1 s-1', 'Reds']
1320    elif varn == 'wbqvshtenac' or varn == 'WBACQVSHTEN':
1321        varvals = ['wbqvshtenac', 'water_budget_qv_shten_ac',  0., 30.,              \
1322          'Water|Budget|accumulated|QV|tendency|due|to|shallow|cu|parameterization', \
1323          'kg kg-1 s-1', 'Reds']
1324    elif varn == 'wbprac' or varn == 'WBACP':
1325        varvals = ['wbprac', 'water_budget_pr_ac',  0., 30.,                         \
1326          'Water|Budget|accumulated|precipitation', 'mm', 'Reds']
1327    elif varn == 'wbpwac' or varn == 'WBACPW':
1328        varvals = ['wbpwac', 'water_budget_pw_ac',  0., 30.,                         \
1329          'Water|Budget|accumulated|water|content', 'mm', 'Reds']
1330    elif varn == 'wbcondtac' or varn == 'WBACCONDT':
1331        varvals = ['wbcondtac', 'water_budget_condt_ac',  0., 30.,                   \
1332          'Water|Budget|accumulated|condensation|and|deposition', 'mm', 'Reds']
1333    elif varn == 'wbqcmac' or varn == 'WBACQCM':
1334        varvals = ['wbqcmac', 'water_budget_qcm_ac',  0., 30.,                       \
1335          'Water|Budget|accumulated|hydrometeor|change|and|convergence', 'mm', 'Reds']
1336    elif varn == 'wbsiac' or varn == 'WBACSI':
1337        varvals = ['wbsiac', 'water_budget_si_ac',  0., 30.,                         \
1338          'Water|Budget|accumulated|hydrometeor|sink', 'mm', 'Reds']
1339    elif varn == 'wbsoac' or varn == 'WBACSO':
1340        varvals = ['wbsoac', 'water_budget_so_ac',  0., 30.,                         \
1341          'Water|Budget|accumulated|hydrometeor|source', 'mm', 'Reds']
1342
1343    elif varn == 'xtime' or varn == 'XTIME':
1344        varvals = ['xtime', 'time',  0., 1.e5, 'time',                               \
1345          'minutes|since|simulation|start', 'Reds']
1346    elif varn == 'x' or varn == 'X':
1347        varvals = ['x', 'x',  0., 100., 'x', '-', 'Reds']
1348    elif varn == 'y' or varn == 'Y':
1349        varvals = ['y', 'y',  0., 100., 'y', '-', 'Blues']
1350    elif varn == 'z' or varn == 'Z':
1351        varvals = ['z', 'z',  0., 100., 'z', '-', 'Greens']
1352    elif varn == 'zg' or varn == 'WRFght' or varn == 'Geopotential height' or        \
1353      varn == 'geop' or varn == 'LGEOP':
1354        varvals = ['zg', 'geopotential_height', 0., 80000., 'geopotential|height',   \
1355          'm2s-2', 'rainbow']
1356    elif varn == 'zmaxth' or varn == 'zmax_th'  or varn == 'LZMAX_TH':
1357        varvals = ['zmaxth', 'thermal_plume_height', 0., 4000.,                     \
1358          'maximum|thermals|plume|height', 'm', 'YlOrRd']
1359    elif varn == 'zmla' or varn == 's_pblh' or varn == 'LS_PBLH':
1360        varvals = ['zmla', 'atmosphere_boundary_layer_thickness', 0., 2500.,         \
1361          'atmosphere|boundary|layer|thickness', 'm', 'Blues']
1362    else:
1363        print errormsg
1364        print '  ' + fname + ": variable '" + varn + "' not defined !!!"
1365        quit(-1)
1366
1367    return varvals
1368
1369#######    #######    #######    #######    #######    #######    #######    #######    #######    #######
1370
1371def check_colorBar(cbarn):
1372    """ Check if the given colorbar exists in matplotlib
1373    """
1374    fname = 'check_colorBar'
1375
1376# Possible color bars
1377    colorbars = ['binary', 'Blues', 'BuGn', 'BuPu', 'gist_yarg', 'GnBu', 'Greens',   \
1378      'Greys', 'Oranges', 'OrRd', 'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',       \
1379      'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'afmhot', 'autumn', 'bone',      \
1380      'cool', 'copper', 'gist_gray', 'gist_heat', 'gray', 'hot', 'pink', 'spring',   \
1381      'summer', 'winter', 'BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr', 'RdBu', \
1382      'RdGy', 'RdYlBu', 'RdYlGn', 'seismic', 'Accent', 'Dark2', 'hsv', 'Paired',     \
1383      'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3', 'spectral', 'gist_earth',        \
1384      'gist_ncar', 'gist_rainbow', 'gist_stern', 'jet', 'brg', 'CMRmap', 'cubehelix',\
1385      'gnuplot', 'gnuplot2', 'ocean', 'rainbow', 'terrain', 'flag', 'prism']
1386
1387    if not searchInlist(colorbars,cbarn):
1388        print warnmsg
1389        print '  ' + fname + ' color bar: "' + cbarn + '" does not exist !!'
1390        print '  a standard one will be use instead !!'
1391
1392    return
1393
1394def units_lunits(u):
1395    """ Fucntion to provide LaTeX equivalences from a given units
1396      u= units to transform
1397    >>> units_lunits('kgkg-1')
1398    '$kgkg^{-1}$'   
1399    """
1400    fname = 'units_lunits'
1401
1402    if u == 'h':
1403        print fname + '_____________________________________________________________'
1404        print units_lunits.__doc__
1405        quit()
1406
1407# Units which does not change
1408    same = ['1', 'category', 'day', 'degrees East', 'degrees Nord', 'degrees North', \
1409      'g', 'hour', 'hPa', 'K', 'Km', 'kg', 'km', 'm', 'minute', 'mm', 'month', 'Pa', \
1410      's', 'second', 'um', 'year', '-']
1411
1412    if searchInlist(same,u):
1413        lu = '$' + u + '$'
1414    elif len(u.split(' ')) > 1 and u.split(' ')[1] == 'since':
1415        uparts = u.split(' ')
1416        ip=0
1417        for up in uparts:
1418            if ip == 0:
1419               lu = '$' + up
1420            else:
1421               lu = lu + '\ ' + up
1422            ip=ip+1
1423        lu = lu + '$'
1424    else:
1425        if u == '': lu='-'
1426        elif u == 'C': lu='$^{\circ}C$'
1427        elif u == 'days': lu='$day$'
1428        elif u == 'degrees_east': lu='$degrees\ East$'
1429        elif u == 'degree_east': lu='$degrees\ East$'
1430        elif u == 'degrees longitude': lu='$degrees\ East$'
1431        elif u == 'degrees latitude': lu='$degrees\ North$'
1432        elif u == 'degrees_north': lu='$degrees\ North$'
1433        elif u == 'degree_north': lu='$degrees\ North$'
1434        elif u == 'deg C': lu='$^{\circ}C$'
1435        elif u == 'degC': lu='$^{\circ}C$'
1436        elif u == 'deg K': lu='$K$'
1437        elif u == 'degK': lu='$K$'
1438        elif u == 'hours': lu='$hour$'
1439        elif u == 'J/kg': lu='$Jkg^{-1}$'
1440        elif u == 'Jkg-1': lu='$Jkg^{-1}$'
1441        elif u == 'K/m': lu='$Km^{-1}$'
1442        elif u == 'Km-1': lu='$Km^{-1}$'
1443        elif u == 'K/s': lu='$Ks^{-1}$'
1444        elif u == 'Ks-1': lu='$Ks^{-1}$'
1445        elif u == 'K s-1': lu='$Ks^{-1}$'
1446        elif u == 'kg/kg': lu='$kgkg^{-1}$'
1447        elif u == 'kgkg-1': lu='$kgkg^{-1}$'
1448        elif u == 'kg kg-1': lu='$kgkg^{-1}$'
1449        elif u == '(kg/kg)/s': lu='$kgkg^{-1}s^{-1}$'
1450        elif u == 'kgkg-1s-1': lu='$kgkg^{-1}s^{-1}$'
1451        elif u == 'kg kg-1 s-1': lu='$kgkg^{-1}s^{-1}$'
1452        elif u == 'kg/m2': lu='$kgm^{-2}$'
1453        elif u == 'kgm-2': lu='$kgm^{-2}$'
1454        elif u == 'kg m-2': lu='$kgm^{-2}$'
1455        elif u == 'Kg m-2': lu='$kgm^{-2}$'
1456        elif u == 'kg/m2/s': lu='$kgm^{-2}s^{-1}$'
1457        elif u == 'kg/(m2*s)': lu='$kgm^{-2}s^{-1}$'
1458        elif u == 'kg/(s*m2)': lu='$kgm^{-2}s^{-1}$'
1459        elif u == 'kgm-2s-1': lu='$kgm^{-2}s^{-1}$'
1460        elif u == 'kg m-2 s-1': lu='$kgm^{-2}s^{-1}$'
1461        elif u == '1/m': lu='$m^{-1}$'
1462        elif u == 'm-1': lu='$m^{-1}$'
1463        elif u == 'm2/s': lu='$m2s^{-1}$'
1464        elif u == 'm2s-1': lu='$m2s^{-1}$'
1465        elif u == 'm2/s2': lu='$m2s^{-2}$'
1466        elif u == 'm/s': lu='$ms^{-1}$'
1467        elif u == 'mmh-3': lu='$mmh^{-3}$'
1468        elif u == 'ms-1': lu='$ms^{-1}$'
1469        elif u == 'm s-1': lu='$ms^{-1}$'
1470        elif u == 'm/s2': lu='$ms^{-2}$'
1471        elif u == 'ms-2': lu='$ms^{-2}$'
1472        elif u == 'minutes': lu='$minute$'
1473        elif u == 'Pa/s': lu='$Pas^{-1}$'
1474        elif u == 'Pas-1': lu='$Pas^{-1}$'
1475        elif u == 'W m-2': lu='$Wm^{-2}$'
1476        elif u == 'Wm-2': lu='$Wm^{-2}$'
1477        elif u == 'W/m2': lu='$Wm^{-2}$'
1478        elif u == '1/s': lu='$s^{-1}$'
1479        elif u == 's-1': lu='$s^{-1}$'
1480        elif u == 'seconds': lu='$second$'
1481        elif u == '%': lu='\%'
1482        else:
1483            print errormsg
1484            print '  ' + fname + ': units "' + u + '" not ready!!!!'
1485            quit(-1)
1486
1487    return lu
1488
1489def ASCII_LaTeX(ln):
1490    """ Function to transform from an ASCII line to LaTeX codification
1491      >>> ASCII_LaTeX('Laboratoire de Météorologie Dynamique però Hovmöller')
1492      Laboratoire de M\'et\'eorologie Dynamique per\`o Hovm\"oller
1493    """
1494    fname='ASCII_LaTeX'
1495
1496    if ln == 'h':
1497        print fname + '_____________________________________________________________'
1498        print ASCII_LaTeX.__doc__
1499        quit()
1500
1501    newln = ln.replace('\\', '\\textbackslash')
1502
1503    newln = newln.replace('á', "\\'a")
1504    newln = newln.replace('é', "\\'e")
1505    newln = newln.replace('í', "\\'i")
1506    newln = newln.replace('ó', "\\'o")
1507    newln = newln.replace('ú', "\\'u")
1508
1509    newln = newln.replace('à', "\\`a")
1510    newln = newln.replace('Ú', "\\`e")
1511    newln = newln.replace('ì', "\\`i")
1512    newln = newln.replace('ò', "\\`o")
1513    newln = newln.replace('ù', "\\`u")
1514
1515    newln = newln.replace('â', "\\^a")
1516    newln = newln.replace('ê', "\\^e")
1517    newln = newln.replace('î', "\\^i")
1518    newln = newln.replace('ÃŽ', "\\^o")
1519    newln = newln.replace('û', "\\^u")
1520
1521    newln = newln.replace('À', '\\"a')
1522    newln = newln.replace('ë', '\\"e')
1523    newln = newln.replace('ï', '\\"i')
1524    newln = newln.replace('ö', '\\"o')
1525    newln = newln.replace('ÃŒ', '\\"u')
1526
1527    newln = newln.replace('ç', '\c{c}')
1528    newln = newln.replace('ñ', '\~{n}')
1529
1530    newln = newln.replace('Á', "\\'A")
1531    newln = newln.replace('É', "\\'E")
1532    newln = newln.replace('Í', "\\'I")
1533    newln = newln.replace('Ó', "\\'O")
1534    newln = newln.replace('Ú', "\\'U")
1535
1536    newln = newln.replace('À', "\\`A")
1537    newln = newln.replace('È', "\\`E")
1538    newln = newln.replace('Ì', "\\`I")
1539    newln = newln.replace('Ò', "\\`O")
1540    newln = newln.replace('Ù', "\\`U")
1541
1542    newln = newln.replace('Â', "\\^A")
1543    newln = newln.replace('Ê', "\\^E")
1544    newln = newln.replace('Î', "\\^I")
1545    newln = newln.replace('Ô', "\\^O")
1546    newln = newln.replace('Û', "\\^U")
1547
1548    newln = newln.replace('Ä', '\\"A')
1549    newln = newln.replace('Ë', '\\"E')
1550    newln = newln.replace('Ï', '\\"I')
1551    newln = newln.replace('Ö', '\\"O')
1552    newln = newln.replace('Ü', '\\"U')
1553
1554    newln = newln.replace('Ç', '\\c{C}')
1555    newln = newln.replace('Ñ', '\\~{N}')
1556
1557    newln = newln.replace('¡', '!`')
1558    newln = newln.replace('¿', '¿`')
1559    newln = newln.replace('%', '\\%')
1560    newln = newln.replace('#', '\\#')
1561    newln = newln.replace('&', '\\&')
1562    newln = newln.replace('$', '\\$')
1563    newln = newln.replace('_', '\\_')
1564    newln = newln.replace('·', '\\textperiodcentered')
1565    newln = newln.replace('<', '$<$')
1566    newln = newln.replace('>', '$>$')
1567    newln = newln.replace('', '*')
1568#    newln = newln.replace('º', '$^{\\circ}$')
1569    newln = newln.replace('ª', '$^{a}$')
1570    newln = newln.replace('º', '$^{o}$')
1571    newln = newln.replace('°', '$^{\\circ}$')
1572    newln = newln.replace('\n', '\\\\\n')
1573    newln = newln.replace('\t', '\\medskip')
1574
1575    return newln
1576
1577def pretty_int(minv,maxv,Nint):
1578    """ Function to plot nice intervals
1579      minv= minimum value
1580      maxv= maximum value
1581      Nint= number of intervals
1582    >>> pretty_int(23.50,67.21,5)
1583    [ 25.  30.  35.  40.  45.  50.  55.  60.  65.]
1584    >>> pretty_int(-23.50,67.21,15)
1585    [  0.  20.  40.  60.]
1586    pretty_int(14.75,25.25,5)
1587    [ 16.  18.  20.  22.  24.]
1588    """ 
1589    fname = 'pretty_int'
1590    nice_int = [1,2,5]
1591
1592#    print 'minv: ',minv,'maxv:',maxv,'Nint:',Nint
1593
1594    interval = np.abs(maxv - minv)
1595
1596    potinterval = np.log10(interval)
1597    Ipotint = int(potinterval)
1598    intvalue = np.float(interval / np.float(Nint))
1599
1600# new
1601    potinterval = np.log10(intvalue)
1602    Ipotint = int(potinterval)
1603
1604#    print 'interval:', interval, 'intavlue:', intvalue, 'potinterval:', potinterval, \
1605#     'Ipotint:', Ipotint, 'intvalue:', intvalue
1606
1607    mindist = 10.e15
1608    for inice in nice_int:
1609#        print inice,':',inice*10.**Ipotint,np.abs(inice*10.**Ipotint - intvalue),mindist
1610        if np.abs(inice*10.**Ipotint - intvalue) < mindist:
1611            mindist = np.abs(inice*10.**Ipotint - intvalue)
1612            closestint = inice
1613
1614    Ibeg = int(minv / (closestint*10.**Ipotint))
1615
1616    values = []
1617    val = closestint*(Ibeg)*10.**(Ipotint)
1618
1619#    print 'closestint:',closestint,'Ibeg:',Ibeg,'val:',val
1620
1621    while val < maxv:
1622        values.append(val)
1623        val = val + closestint*10.**Ipotint
1624
1625    return np.array(values, dtype=np.float)
1626
1627def DegGradSec_deg(grad,deg,sec):
1628    """ Function to transform from a coordinate in grad deg sec to degrees (decimal)
1629    >>> DegGradSec_deg(39.,49.,26.)
1630    39.8238888889
1631    """
1632    fname = 'DegGradSec_deg'
1633
1634    if grad == 'h':
1635        print fname + '_____________________________________________________________'
1636        print DegGradSec_deg.__doc__
1637        quit()
1638
1639    deg = grad + deg/60. + sec/3600.
1640
1641    return deg
1642
1643def intT2dt(intT,tu):
1644    """ Function to provide an 'timedelta' object from a given interval value
1645      intT= interval value
1646      tu= interval units, [tu]= 'd': day, 'w': week, 'h': hour, 'i': minute, 's': second,
1647        'l': milisecond
1648
1649      >>> intT2dt(3.5,'s')
1650      0:00:03.500000
1651
1652      >>> intT2dt(3.5,'w')
1653      24 days, 12:00:00
1654    """
1655    import datetime as dt 
1656
1657    fname = 'intT2dt'
1658
1659    if tu == 'w':
1660        dtv = dt.timedelta(weeks=np.float(intT))
1661    elif tu == 'd':
1662        dtv = dt.timedelta(days=np.float(intT))
1663    elif tu == 'h':
1664        dtv = dt.timedelta(hours=np.float(intT))
1665    elif tu == 'i':
1666        dtv = dt.timedelta(minutes=np.float(intT))
1667    elif tu == 's':
1668        dtv = dt.timedelta(seconds=np.float(intT))
1669    elif tu == 'l':
1670        dtv = dt.timedelta(milliseconds=np.float(intT))
1671    else:
1672        print errormsg
1673        print '  ' + fname + ': time units "' + tu + '" not ready!!!!'
1674        quit(-1)
1675
1676    return dtv
1677
1678def lonlat_values(mapfile,lonvn,latvn):
1679    """ Function to obtain the lon/lat matrices from a given netCDF file
1680    lonlat_values(mapfile,lonvn,latvn)
1681      [mapfile]= netCDF file name
1682      [lonvn]= variable name with the longitudes
1683      [latvn]= variable name with the latitudes
1684    """
1685
1686    fname = 'lonlat_values'
1687
1688    if mapfile == 'h':
1689        print fname + '_____________________________________________________________'
1690        print lonlat_values.__doc__
1691        quit()
1692
1693    if not os.path.isfile(mapfile):
1694        print errormsg
1695        print '  ' + fname + ": map file '" + mapfile + "' does not exist !!"
1696        quit(-1)
1697
1698    ncobj = NetCDFFile(mapfile, 'r')
1699    lonobj = ncobj.variables[lonvn]
1700    latobj = ncobj.variables[latvn]
1701
1702    if len(lonobj.shape) == 3:
1703        lonv = lonobj[0,:,:]
1704        latv = latobj[0,:,:]
1705    elif len(lonobj.shape) == 2:
1706        lonv = lonobj[:,:]
1707        latv = latobj[:,:]
1708    elif len(lonobj.shape) == 1:
1709        lon0 = lonobj[:]
1710        lat0 = latobj[:]
1711        lonv = np.zeros( (len(lat0),len(lon0)), dtype=np.float )
1712        latv = np.zeros( (len(lat0),len(lon0)), dtype=np.float )
1713        for iy in range(len(lat0)):
1714            lonv[iy,:] = lon0
1715        for ix in range(len(lon0)):
1716            latv[:,ix] = lat0
1717    else:
1718        print errormsg
1719        print '  ' + fname + ': lon/lat variables shape:',lonobj.shape,'not ready!!'
1720        quit(-1)
1721
1722    return lonv, latv
1723
1724def date_CFtime(ind,refd,tunits):
1725    """ Function to transform from a given date object a CF-convention time
1726      ind= date object to transform
1727      refd= reference date
1728      tunits= units for time
1729        >>> date_CFtime(dt.datetime(1976,02,17,08,30,00), dt.datetime(1949,12,01,00,00,00), 'seconds')
1730        827224200.0
1731    """
1732    import datetime as dt
1733
1734    fname = 'date_CFtime'
1735
1736    dt = ind - refd
1737
1738    if tunits == 'weeks':
1739        value = dt.days/7. + dt.seconds/(3600.*24.*7.)
1740    elif tunits == 'days':
1741        value = dt.days + dt.seconds/(3600.*24.)
1742    elif tunits == 'hours':
1743        value = dt.days*24. + dt.seconds/(3600.)
1744    elif tunits == 'minutes':
1745        value = dt.days*24.*60. + dt.seconds/(60.)
1746    elif tunits == 'seconds':
1747        value = dt.days*24.*3600. + dt.seconds
1748    elif tunits == 'milliseconds':
1749        value = dt.days*24.*3600.*1000. + dt.seconds*1000.
1750    else:
1751        print errormsg
1752        print '  ' + fname + ': reference time units "' + trefu + '" not ready!!!!'
1753        quit(-1)
1754
1755    return value
1756
1757def pot_values(values, uvals):
1758    """ Function to modify a seies of values by their potency of 10
1759      pot_values(values, uvals)
1760      values= values to modify
1761      uvals= units of the values
1762      >>> vals = np.sin(np.arange(20)*np.pi/5.+0.01)*10.e-5
1763      >>> pot_values(vals,'ms-1')
1764      (array([  0.00000000e+00,   5.87785252e-01,   9.51056516e-01,
1765         9.51056516e-01,   5.87785252e-01,   1.22464680e-16,
1766        -5.87785252e-01,  -9.51056516e-01,  -9.51056516e-01,
1767        -5.87785252e-01,  -2.44929360e-16,   5.87785252e-01,
1768         9.51056516e-01,   9.51056516e-01,   5.87785252e-01,
1769         3.67394040e-16,  -5.87785252e-01,  -9.51056516e-01,
1770        -9.51056516e-01,  -5.87785252e-01]), -4, 'x10e-4 ms-1', 'x10e-4')
1771    """
1772
1773    fname = 'pot_values'
1774
1775    if np.min(values) != 0.:
1776        potmin = int( np.log10( np.abs(np.min(values)) ) )
1777    else:
1778        potmin = 0
1779
1780    if np.max(values) != 0.:
1781        potmax = int( np.log10( np.abs(np.max(values)) ) )
1782    else:
1783        potmax = 0
1784
1785    if potmin * potmax > 9:
1786        potval = -np.min([np.abs(potmin), np.abs(potmax)]) * np.abs(potmin) / potmin
1787
1788        newvalues = values*10.**potval
1789        potvalue = - potval
1790        potS = 'x10e' + str(potvalue)
1791        newunits = potS + ' ' + uvals
1792    else:
1793        newvalues = values
1794        potvalue = None
1795        potS = ''
1796        newunits = uvals
1797
1798    return newvalues, potvalue, newunits, potS
1799
1800def CFtimes_plot(timev,units,kind,tfmt):
1801    """ Function to provide a list of string values from a CF time values in order
1802      to use them in a plot, according to the series of characteristics.
1803      String outputs will be suited to the 'human-like' output
1804        timev= time values (real values)
1805        units= units string according to CF conventions ([tunits] since
1806          [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]])
1807        kind= kind of output
1808          'Nval': according to a given number of values as 'Nval',[Nval]
1809          'exct': according to an exact time unit as 'exct',[tunit];
1810            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1811              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1812              'l': milisecond
1813        tfmt= desired format
1814          >>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'Nval,5',"%Y/%m/%d %H:%M:%S")
1815            0.0 1979/12/01 00:00:00
1816            24.75 1979/12/02 00:45:00
1817            49.5 1979/12/03 01:30:00
1818            74.25 1979/12/04 02:15:00
1819            99.0 1979/12/05 03:00:00
1820          >>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'exct,2,d',"%Y/%m/%d %H:%M:%S")
1821            0.0 1979/12/01 00:00:00
1822            48.0 1979/12/03 00:00:00
1823            96.0 1979/12/05 00:00:00
1824            144.0 1979/12/07 00:00:00
1825    """ 
1826    import datetime as dt 
1827
1828# Seconds between 0001 and 1901 Jan - 01
1829    secs0001_1901=59958144000.
1830
1831    fname = 'CFtimes_plot'
1832
1833    if timev == 'h':
1834        print fname + '_____________________________________________________________'
1835        print CFtimes_plot.__doc__
1836        quit()
1837
1838# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
1839##
1840    trefT = units.find(':')
1841    txtunits = units.split(' ')
1842    Ntxtunits = len(txtunits)
1843
1844    if Ntxtunits == 3:
1845        Srefdate = txtunits[Ntxtunits - 1]
1846    else:
1847        Srefdate = txtunits[Ntxtunits - 2]
1848
1849    if not trefT == -1:
1850#        print '  ' + fname + ': refdate with time!'
1851        if Ntxtunits == 3:
1852            refdate = datetimeStr_datetime(Srefdate)
1853        else:
1854            refdate = datetimeStr_datetime(Srefdate + '_' + txtunits[Ntxtunits - 1])
1855    else:
1856        refdate = datetimeStr_datetime(Srefdate + '_00:00:00')
1857
1858    trefunits=units.split(' ')[0]
1859    if trefunits == 'weeks':
1860        trefu = 'w'
1861    elif trefunits == 'days':
1862        trefu = 'd'
1863    elif trefunits == 'hours':
1864        trefu = 'h'
1865    elif trefunits == 'minutes':
1866        trefu = 'm'
1867    elif trefunits == 'seconds':
1868        trefu = 's'
1869    elif trefunits == 'milliseconds':
1870        trefu = 'l'
1871    else:
1872        print errormsg
1873        print '  ' + fname + ': reference time units "' + trefu + '" not ready!!!!'
1874        quit(-1)
1875
1876    okind=kind.split(',')[0]
1877    dtv = len(timev)
1878 
1879    if refdate.year == 1:
1880        print warnmsg
1881        print '  ' + fname + ': changing reference date: ',refdate,                  \
1882          'to 1901-01-01_00:00:00 !!!'
1883        refdate = datetimeStr_datetime('1901-01-01_00:00:00')
1884        if trefu == 'w': timev = timev - secs0001_1901/(7.*24.*3600.)
1885        if trefu == 'd': timev = timev - secs0001_1901/(24.*3600.)
1886        if trefu == 'h': timev = timev - secs0001_1901/(3600.)
1887        if trefu == 'm': timev = timev - secs0001_1901/(60.)
1888        if trefu == 's': timev = timev - secs0001_1901
1889        if trefu == 'l': timev = timev - secs0001_1901*1000.
1890
1891    firstT = timev[0]
1892    lastT = timev[dtv-1]
1893
1894# First and last times as datetime objects
1895    firstTdt = timeref_datetime(refdate, firstT, trefunits)
1896    lastTdt = timeref_datetime(refdate, lastT, trefunits)
1897
1898# First and last times as [year, mon, day, hour, minut, second] vectors
1899    firstTvec = np.zeros((6), dtype= np.float)
1900    lastTvec = np.zeros((6), dtype= np.float)
1901    chTvec = np.zeros((6), dtype= bool)
1902
1903    firstTvec = np.array([firstTdt.year, firstTdt.month, firstTdt.day, firstTdt.hour,\
1904      firstTdt.minute, firstTdt.second])
1905    lastTvec = np.array([lastTdt.year, lastTdt.month, lastTdt.day, lastTdt.hour,     \
1906      lastTdt.minute, lastTdt.second])
1907
1908    chdate= lastTvec - firstTvec
1909    chTvec = np.where (chdate != 0., True, False)
1910
1911    timeout = []
1912    if okind == 'Nval':
1913        nvalues = int(kind.split(',')[1])
1914        intervT = (lastT - firstT)/(nvalues-1)
1915        dtintervT = intT2dt(intervT, trefu)
1916
1917        for it in range(nvalues):
1918            timeout.append(firstTdt + dtintervT*it)
1919    elif okind == 'exct':
1920        Nunits = int(kind.split(',')[1])
1921        tu = kind.split(',')[2]
1922
1923# Generic incremental dt [seconds] according to all the possibilities ['c', 'y', 'm',
1924#   'w', 'd', 'h', 'i', 's', 'l'], some of them approximated (because they are not
1925#   already necessary!)
1926        basedt = np.zeros((9), dtype=np.float)
1927        basedt[0] = (365.*100. + 25.)*24.*3600.
1928        basedt[1] = 365.*24.*3600.
1929        basedt[2] = 31.*24.*3600.
1930        basedt[3] = 7.*24.*3600.
1931        basedt[4] = 24.*3600.
1932        basedt[5] = 3600.
1933        basedt[6] = 60.
1934        basedt[7] = 1.
1935        basedt[8] = 1000.
1936
1937# Increment according to the units of the CF dates
1938        if trefunits == 'weeks':
1939            basedt = basedt/(7.*24.*3600.)
1940        elif trefunits == 'days':
1941            basedt = basedt/(24.*3600.)
1942        elif trefunits == 'hours':
1943            basedt = basedt/(3600.)
1944        elif trefunits == 'minutes':
1945            basedt = basedt/(60.)
1946        elif trefunits == 'seconds':
1947            basedt = basedt
1948        elif trefunits == 'milliseconds':
1949            basedt = basedt*1000.
1950
1951        if tu == 'c':
1952            ti = firstTvec[0]
1953            tf = lastTvec[0] 
1954            centi = firstTvec[0] / 100
1955
1956            for it in range((tf - ti)/(Nunits*100) + 1):
1957                timeout.append(dt.datetime(centi+it*Nunits*100, 1, 1, 0, 0, 0))
1958        elif tu == 'y':
1959            ti = firstTvec[0]
1960            tf = lastTvec[0]
1961            yeari = firstTvec[0]
1962
1963            for it in range((tf - ti)/(Nunits) + 1):
1964                timeout.append(dt.datetime(yeari+it*Nunits, 1, 1, 0, 0, 0))
1965        elif tu == 'm':
1966            ti = firstTvec[1]
1967            tf = lastTvec[1]
1968            yr = firstTvec[0]
1969            mon = firstTvec[1]
1970
1971            for it in range((tf - ti)/(Nunits) + 1):
1972                mon = mon+it*Nunits
1973                if mon > 12:
1974                    yr = yr + 1
1975                    mon = 1
1976
1977                timeout.append(dt.datetime(yr, mon, 1, 0, 0, 0))
1978        elif tu == 'w':
1979            datev=firstTdt
1980            it=0
1981            while datev <= lastTdt:
1982                datev = firstTdt + dt.timedelta(days=7*Nunits*it)
1983                timeout.append(datev)
1984                it = it + 1
1985        elif tu == 'd':
1986#            datev=firstTdt
1987            yr = firstTvec[0]
1988            mon = firstTvec[1]
1989            day = firstTvec[2]
1990
1991            if np.sum(firstTvec[2:5]) > 0:
1992                firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0)
1993                datev = dt.datetime(yr, mon, day+1, 0, 0, 0)
1994            else:
1995                firstTdt = dt.datetime(yr, mon, day, 0, 0, 0)
1996                datev = dt.datetime(yr, mon, day, 0, 0, 0)
1997
1998            it=0
1999            while datev <= lastTdt:
2000                datev = firstTdt + dt.timedelta(days=Nunits*it)
2001                timeout.append(datev)
2002                it = it + 1
2003        elif tu == 'h':
2004            datev=firstTdt
2005            yr = firstTvec[0]
2006            mon = firstTvec[1]
2007            day = firstTvec[2]
2008            hour = firstTvec[3]
2009
2010            if np.sum(firstTvec[4:5]) > 0 or np.mod(hour,Nunits) != 0:
2011                tadvance = 2*Nunits
2012                if tadvance >= 24:
2013                    firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0)
2014                    datev = dt.datetime(yr, mon, day+1, 0, 0, 0)
2015                else:
2016                    firstTdt = dt.datetime(yr, mon, day, Nunits, 0, 0)
2017                    datev = dt.datetime(yr, mon, day, Nunits, 0, 0)
2018            else:
2019                firstTdt = dt.datetime(yr, mon, day, hour, 0, 0)
2020                datev = dt.datetime(yr, mon, day, hour, 0, 0)
2021
2022            it=0
2023            while datev <= lastTdt:
2024                datev = firstTdt + dt.timedelta(seconds=Nunits*3600*it)
2025                timeout.append(datev)
2026                it = it + 1                 
2027        elif tu == 'i':
2028            datev=firstTdt
2029            yr = firstTvec[0]
2030            mon = firstTvec[1]
2031            day = firstTvec[2]
2032            hour = firstTvec[3]
2033            minu = firstTvec[4]
2034
2035            if firstTvec[5] > 0 or np.mod(minu,Nunits) != 0:
2036                tadvance = 2*Nunits
2037                if tadvance >= 60:
2038                    firstTdt = dt.datetime(yr, mon, day, hour, 0, 0)
2039                    datev = dt.datetime(yr, mon, day, hour, 0, 0)
2040                else:
2041                    firstTdt = dt.datetime(yr, mon, day, hour, Nunits, 0)
2042                    datev = dt.datetime(yr, mon, day, hour, Nunits, 0)
2043            else:
2044                firstTdt = dt.datetime(yr, mon, day, hour, minu, 0)
2045                datev = dt.datetime(yr, mon, day, hour, minu, 0)
2046            it=0
2047            while datev <= lastTdt:
2048                datev = firstTdt + dt.timedelta(seconds=Nunits*60*it)
2049                timeout.append(datev)
2050                it = it + 1                 
2051        elif tu == 's':
2052            datev=firstTdt
2053            it=0
2054            while datev <= lastTdt:
2055                datev = firstTdt + dt.timedelta(seconds=Nunits)
2056                timeout.append(datev)
2057                it = it + 1                 
2058        elif tu == 'l':
2059            datev=firstTdt
2060            it=0
2061            while datev <= lastTdt:
2062                datev = firstTdt + dt.timedelta(seconds=Nunits*it/1000.)
2063                timeout.append(datev)
2064                it = it + 1
2065        else:
2066            print errormsg
2067            print '  '  + fname + ': exact units "' + tu + '" not ready!!!!!'
2068            quit(-1)
2069
2070    else:
2071        print errormsg
2072        print '  ' + fname + ': output kind "' + okind + '" not ready!!!!'
2073        quit(-1)
2074
2075    dtout = len(timeout)
2076
2077    timeoutS = []
2078    timeoutv = np.zeros((dtout), dtype=np.float)
2079
2080    for it in range(dtout):
2081        timeoutS.append(timeout[it].strftime(tfmt))
2082        timeoutv[it] = date_CFtime(timeout[it], refdate, trefunits)
2083       
2084#        print it,':',timeoutv[it], timeoutS[it]
2085
2086    if len(timeoutv) < 1 or len(timeoutS) < 1:
2087        print errormsg
2088        print '  ' + fname + ': no time values are generated!'
2089        print '    values passed:',timev
2090        print '    units:',units
2091        print '    desired kind:',kind
2092        print '    format:',tfmt
2093        print '    function values ___ __ _'
2094        print '      reference date:',refdate
2095        print '      first date:',firstTdt
2096        print '      last date:',lastTdt
2097        print '      icrement:',basedt,trefunits
2098
2099        quit(-1)
2100
2101    return timeoutv, timeoutS
2102
2103def color_lines(Nlines):
2104    """ Function to provide a color list to plot lines
2105    color_lines(Nlines)
2106      Nlines= number of lines
2107    """
2108
2109    fname = 'color_lines'
2110
2111    colors = ['r', 'b', 'g', 'p', 'g']
2112
2113    colorv = []
2114
2115    colorv.append('k')
2116    for icol in range(Nlines):
2117        colorv.append(colors[icol])
2118
2119
2120    return colorv
2121
2122def output_kind(kindf, namef, close):
2123    """ Function to generate the output of the figure
2124      kindf= kind of the output
2125        null: show in screen
2126        [jpg/pdf/png/ps]: standard output types
2127      namef= name of the figure (without extension)
2128      close= if the graph has to be close or not [True/False]
2129    """
2130    fname = 'output_kind'
2131
2132    if kindf == 'h':
2133        print fname + '_____________________________________________________________'
2134        print output_kind.__doc__
2135        quit()
2136
2137    if kindf == 'null':
2138        print 'showing figure...'
2139        plt.show()
2140    elif kindf == 'gif':
2141        plt.savefig(namef + ".gif")
2142        if close: print "Successfully generation of figure '" + namef + ".jpg' !!!"
2143    elif kindf == 'jpg':
2144        plt.savefig(namef + ".jpg")
2145        if close: print "Successfully generation of figure '" + namef + ".jpg' !!!"
2146    elif kindf == 'pdf':
2147        plt.savefig(namef + ".pdf")
2148        if close: print "Successfully generation of figure '" + namef + ".pdf' !!!"
2149    elif kindf == 'png':
2150        plt.savefig(namef + ".png")
2151        if close: print "Successfully generation of figure '" + namef + ".png' !!!"
2152    elif kindf == 'ps':
2153        plt.savefig(namef + ".ps")
2154        if close: print "Successfully generation of figure '" + namef + ".ps' !!!"
2155    else:
2156        print errormsg
2157        print '  ' + fname + ' output format: "' + kindf + '" not ready !!'
2158        print errormsg
2159        quit(-1)
2160
2161    if close:
2162        plt.close()
2163
2164    return
2165
2166def check_arguments(funcname,Nargs,args,char,expectargs):
2167    """ Function to check the number of arguments if they are coincident
2168    check_arguments(funcname,Nargs,args,char)
2169      funcname= name of the function/program to check
2170      Nargs= theoretical number of arguments
2171      args= passed arguments
2172      char= character used to split the arguments
2173    """
2174
2175    fname = 'check_arguments'
2176
2177    Nvals = len(args.split(char))
2178    if Nvals != Nargs:
2179        print errormsg
2180        print '  ' + fname + ': wrong number of arguments:',Nvals," passed to  '",   \
2181          funcname, "' which requires:",Nargs,'!!'
2182        print '    given arguments:',args.split(char)
2183        print '    expected arguments:',expectargs
2184        quit(-1)
2185
2186    return
2187
2188def Str_Bool(val):
2189    """ Function to transform from a String value to a boolean one
2190    >>> Str_Bool('True')
2191    True
2192    >>> Str_Bool('0')
2193    False
2194    >>> Str_Bool('no')
2195    False
2196    """
2197
2198    fname = 'Str_Bool'
2199
2200    if val == 'True' or val == '1' or val == 'yes': 
2201        boolv = True
2202    elif val == 'False' or val == '0' or val== 'no':
2203        boolv = False
2204    else:
2205        print errormsg
2206        print '  ' + fname + ": value '" + val + "' not ready!!"
2207        quit(-1)
2208
2209    return boolv
2210
2211def plot_TimeSeries(valtimes, vunits, tunits, hfileout, vtit, ttit, tkind, tformat,  \
2212  tit, linesn, lloc, kfig):
2213    """ Function to draw time-series
2214      valtimes= list of arrays to plot [vals1[1values, 1times], [...,valsM[Mvals,Mtimes]])
2215      vunits= units of the values
2216      tunits= units of the times
2217      hfileout= header of the output figure. Final name: [hfileout]_[vtit].[kfig]
2218      vtit= variable title to be used in the graph
2219      ttit= time title to be used in the graph
2220      tkind= kind of time values to appear in the x-axis
2221        'Nval': according to a given number of values as 'Nval',[Nval]
2222        'exct': according to an exact time unit as 'exct',[tunit];
2223          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2224            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2225            'l': milisecond
2226      tformat= desired format of times
2227      tit= title of the graph
2228      linesn= list of values fot the legend
2229      lloc= location of the legend (-1, autmoatic)
2230        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2231        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2232        9: 'upper center', 10: 'center'
2233      kfig= type of figure: jpg, png, pds, ps
2234    """
2235    fname = 'plot_TimeSeries'
2236
2237    if valtimes == 'h':
2238        print fname + '_____________________________________________________________'
2239        print plot_TimeSeries.__doc__
2240        quit()
2241
2242
2243# Canging line kinds every 7 lines (end of standard colors)
2244    linekinds=['.-','x-','o-']
2245
2246    Nlines = len(valtimes)
2247
2248    Nvalues = []
2249    Ntimes = []
2250
2251    for il in range(Nlines):
2252        array = valtimes[il]
2253
2254        if Nlines == 1:
2255            print warnmsg
2256            print '  ' + fname + ': drawing only one line!'
2257
2258            Nvalues.append(array.shape[1])
2259            Ntimes.append(array.shape[0])
2260            tmin = np.min(array[1])
2261            tmax = np.max(array[1])
2262            vmin = np.min(array[0])
2263            vmax = np.max(array[0])
2264        else:
2265            Nvalues.append(array.shape[1])
2266            Ntimes.append(array.shape[0])
2267            tmin = np.min(array[1,:])
2268            tmax = np.max(array[1,:])
2269            vmin = np.min(array[0,:])
2270            vmax = np.max(array[0,:])
2271
2272        if il == 0:
2273            xmin = tmin
2274            xmax = tmax
2275            ymin = vmin
2276            ymax = vmax
2277        else:
2278            if tmin < xmin: xmin = tmin
2279            if tmax > xmax: xmax = tmax
2280            if vmin < ymin: ymin = vmin
2281            if vmax > ymax: ymax = vmax
2282
2283    dx = np.max(Ntimes)
2284    dy = np.min(Nvalues)
2285
2286    plt.rc('text', usetex=True)
2287
2288    print vtit
2289    if vtit == 'ps':
2290        plt.ylim(98000.,ymax)
2291    else:
2292        plt.ylim(ymin,ymax)
2293
2294    plt.xlim(xmin,xmax)
2295#    print 'x lim:',xmin,xmax
2296#    print 'y lim:',ymin,ymax
2297
2298    N7lines=0
2299    for il in range(Nlines):
2300        array = valtimes[il]
2301        if vtit == 'ps':
2302            array[0,:] = np.where(array[0,:] < 98000., None, array[0,:])
2303        plt.plot(array[1,:],array[0,:], linekinds[N7lines], label= linesn[il])
2304        if il == 6: N7lines = N7lines + 1
2305
2306    timevals = np.arange(xmin,xmax)*1.
2307
2308    tpos, tlabels = CFtimes_plot(timevals, tunits, tkind, tformat)
2309
2310    if len(tpos) > 10:
2311        print warnmsg
2312        print '  ' + fname + ': with "' + tkind + '" there are', len(tpos), 'xticks !'
2313
2314    plt.xticks(tpos, tlabels)
2315#    plt.Axes.set_xticklabels(tlabels)
2316
2317    plt.legend(loc=lloc)
2318    plt.xlabel(ttit)
2319    plt.ylabel(vtit + " (" + vunits + ")")
2320    plt.title(tit).replace('_','\_').replace('&','\&')
2321
2322    figname = hfileout + '_' + vtit
2323   
2324    output_kind(kfig, figname, True)
2325
2326    return
2327
2328#Nt = 50
2329#Nlines = 3
2330
2331#vtvalsv = []
2332
2333#valsv = np.zeros((2,Nt), dtype=np.float)
2334## First
2335#valsv[0,:] = np.arange(Nt)
2336#valsv[1,:] = np.arange(Nt)*180.
2337#vtvalsv.append(valsv)
2338#del(valsv)
2339
2340#valsv = np.zeros((2,Nt/2), dtype=np.float)
2341## Second
2342#valsv[0,:] = np.arange(Nt/2)
2343#valsv[1,:] = np.arange(Nt/2)*180.*2.
2344#vtvalsv.append(valsv)
2345#del(valsv)
2346
2347#valsv = np.zeros((2,Nt/4), dtype=np.float)
2348## Third
2349#valsv[0,:] = np.arange(Nt/4)
2350#valsv[1,:] = np.arange(Nt/4)*180.*4.
2351#vtvalsv.append(valsv)
2352#del(valsv)
2353
2354#varu='mm'
2355#timeu='seconds'
2356
2357#title='test'
2358#linesname = ['line 1', 'line 2', 'line 3']
2359
2360#plot_TimeSeries(vtvalsv, units_lunits(varu), timeu, 'test', 'vartest', 'time', title, linesname, 'png')
2361#quit()
2362
2363def plot_points(xval, yval, ifile, vtit, kfig, mapv):
2364    """ plotting points
2365      [x/yval]: x,y values to plot
2366      vn,vm= minmum and maximum values to plot
2367      unit= units of the variable
2368      ifile= name of the input file
2369      vtit= title of the variable
2370      kfig= kind of figure (jpg, pdf, png)
2371      mapv= map characteristics: [proj],[res]
2372        see full documentation: http://matplotlib.org/basemap/
2373        [proj]: projection
2374          * 'cyl', cilindric
2375        [res]: resolution:
2376          * 'c', crude
2377          * 'l', low
2378          * 'i', intermediate
2379          * 'h', high
2380          * 'f', full
2381    """
2382    fname = 'plot_points'
2383
2384    dx=xval.shape[0]
2385    dy=yval.shape[0]
2386
2387    plt.rc('text', usetex=True)
2388
2389    if not mapv is None:
2390        lon00 = np.where(xval[:] < 0., 360. + olon[:], olon[:])
2391        lat00 = yval[:]
2392        lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2393        lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2394
2395        for iy in range(len(lat00)):
2396            lon0[iy,:] = lon00
2397        for ix in range(len(lon00)):
2398            lat0[:,ix] = lat00
2399
2400        map_proj=mapv.split(',')[0]
2401        map_res=mapv.split(',')[1]
2402
2403        nlon = lon0[0,0]
2404        xlon = lon0[dy-1,dx-1]
2405        nlat = lat0[0,0]
2406        xlat = lat0[dy-1,dx-1]
2407
2408        lon2 = lon0[dy/2,dx/2]
2409        lat2 = lat0[dy/2,dx/2]
2410
2411        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
2412          xlon, ',', xlat
2413
2414        if map_proj == 'cyl':
2415            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
2416              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2417        elif map_proj == 'lcc':
2418            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
2419              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2420
2421        lons, lats = np.meshgrid(olon[:], olat[:])
2422        lons = np.where(lons < 0., lons + 360., lons)
2423
2424        x,y = m(lons,lats)
2425        plt.plot(x,y)
2426
2427        m.drawcoastlines()
2428
2429        meridians = pretty_int(nlon,xlon,5)
2430        m.drawmeridians(meridians,labels=[True,False,False,True])
2431
2432        parallels = pretty_int(nlat,xlat,5)
2433        m.drawparallels(parallels,labels=[False,True,True,False])
2434    else:
2435#        plt.xlim(0,dx-1)
2436#        plt.ylim(0,dy-1)
2437
2438        plt.plot(xval, yval, '.')
2439
2440    figname = ifile.replace('.','_') + '_' + vtit
2441    graphtit = vtit.replace('_','\_').replace('&','\&')
2442
2443    plt.title(graphtit)
2444   
2445    output_kind(kfig, figname, True)
2446
2447    return
2448
2449def plot_2Dfield(varv,dimn,colorbar,vn,vx,unit,olon,olat,ifile,vtit,zvalue,time,tk,  \
2450   tkt,tobj,tvals,tind,kfig,mapv,reva):
2451    """ Adding labels and other staff to the graph
2452      varv= 2D values to plot
2453      dimn= dimension names to plot
2454      colorbar= name of the color bar to use
2455      vn,vm= minmum and maximum values to plot
2456      unit= units of the variable
2457      olon,olat= longitude, latitude objects
2458      ifile= name of the input file
2459      vtit= title of the variable
2460      zvalue= value on the z axis
2461      time= value on the time axis
2462      tk= kind of time (WRF)
2463      tkt= kind of time taken
2464      tobj= tim object
2465      tvals= values of the time variable
2466      tind= time index
2467      kfig= kind of figure (jpg, pdf, png)
2468      mapv= map characteristics: [proj],[res]
2469        see full documentation: http://matplotlib.org/basemap/
2470        [proj]: projection
2471          * 'cyl', cilindric
2472        [res]: resolution:
2473          * 'c', crude
2474          * 'l', low
2475          * 'i', intermediate
2476          * 'h', high
2477          * 'f', full
2478      reva= reverse the axes (x-->y, y-->x)
2479    """
2480##    import matplotlib as mpl
2481##    mpl.use('Agg')
2482##    import matplotlib.pyplot as plt
2483
2484    if reva:
2485        print '  reversing the axes of the figure (x-->y, y-->x)!!'
2486        varv = np.transpose(varv)
2487        dimn0 = []
2488        dimn0.append(dimn[1] + '')
2489        dimn0.append(dimn[0] + '')
2490        dimn = dimn0
2491
2492    fname = 'plot_2Dfield'
2493    dx=varv.shape[1]
2494    dy=varv.shape[0]
2495
2496    plt.rc('text', usetex=True)
2497#    plt.rc('font', family='serif')
2498
2499    if not mapv is None:
2500        if len(olon[:].shape) == 3:
2501            lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,])
2502            lat0 = olat[0,]
2503        elif len(olon[:].shape) == 2:
2504            lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
2505            lat0 = olat[:]
2506        elif len(olon[:].shape) == 1:
2507            lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
2508            lat00 = olat[:]
2509            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2510            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2511
2512            for iy in range(len(lat00)):
2513                lon0[iy,:] = lon00
2514            for ix in range(len(lon00)):
2515                lat0[:,ix] = lat00
2516
2517        map_proj=mapv.split(',')[0]
2518        map_res=mapv.split(',')[1]
2519
2520        nlon = lon0[0,0]
2521        xlon = lon0[dy-1,dx-1]
2522        nlat = lat0[0,0]
2523        xlat = lat0[dy-1,dx-1]
2524
2525        lon2 = lon0[dy/2,dx/2]
2526        lat2 = lat0[dy/2,dx/2]
2527
2528        print '    lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
2529          xlon, ',', xlat
2530
2531        if map_proj == 'cyl':
2532            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
2533              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2534        elif map_proj == 'lcc':
2535            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
2536              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2537
2538        if len(olon[:].shape) == 1:
2539            lons, lats = np.meshgrid(olon[:], olat[:])
2540        else:
2541            lons = olon[0,:]
2542            lats = olat[:,0]
2543 
2544        lons = np.where(lons < 0., lons + 360., lons)
2545
2546        x,y = m(lons,lats)
2547        plt.pcolormesh(x,y,varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
2548        cbar = plt.colorbar()
2549
2550        m.drawcoastlines()
2551#        if (nlon > 180. or xlon > 180.):
2552#            nlon0 = nlon
2553#            xlon0 = xlon
2554#            if (nlon > 180.): nlon0 = nlon - 360.
2555#            if (xlon > 180.): xlon0 = xlon - 360.
2556#            meridians = pretty_int(nlon0,xlon0,5)           
2557#            meridians = np.where(meridians < 0., meridians + 360., meridians)
2558#        else:
2559#            meridians = pretty_int(nlon,xlon,5)
2560
2561        meridians = pretty_int(nlon,xlon,5)
2562
2563        m.drawmeridians(meridians,labels=[True,False,False,True])
2564        parallels = pretty_int(nlat,xlat,5)
2565        m.drawparallels(parallels,labels=[False,True,True,False])
2566
2567    else:
2568        plt.xlim(0,dx-1)
2569        plt.ylim(0,dy-1)
2570
2571        plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
2572        cbar = plt.colorbar()
2573
2574        plt.xlabel(dimn[1].replace('_','\_'))
2575        plt.ylabel(dimn[0].replace('_','\_'))
2576
2577# set the limits of the plot to the limits of the data
2578#    plt.axis([x.min(), x.max(), y.min(), y.max()])
2579
2580#    plt.plot(varv)
2581    cbar.set_label(unit)
2582
2583    figname = ifile.replace('.','_') + '_' + vtit
2584    graphtit = vtit.replace('_','\_').replace('&','\&')
2585
2586    if zvalue != 'null':
2587        graphtit = graphtit + ' at z= ' + zvalue
2588        figname = figname + '_z' + zvalue
2589    if tkt == 'tstep':
2590        graphtit = graphtit + ' at time-step= ' + time.split(',')[1]
2591        figname = figname + '_t' + time.split(',')[1].zfill(4)
2592    elif tkt == 'CFdate':
2593        graphtit = graphtit + ' at ' + tobj.strfmt("%Y%m%d%H%M%S")
2594        figname = figname + '_t' + tobj.strfmt("%Y%m%d%H%M%S")
2595
2596    if tk == 'WRF':
2597#        datev = str(timevals[timeind][0:9])
2598        datev = tvals[tind][0] + tvals[tind][1] + tvals[tind][2] + \
2599          timevals[timeind][3] + timevals[timeind][4] + timevals[timeind][5] +       \
2600          timevals[timeind][6] + timevals[timeind][7] + timevals[timeind][8] +       \
2601          timevals[timeind][9]
2602#        timev = str(timevals[timeind][11:18])
2603        timev = timevals[timeind][11] + timevals[timeind][12] +                      \
2604          timevals[timeind][13] + timevals[timeind][14] + timevals[timeind][15] +    \
2605          timevals[timeind][16] + timevals[timeind][17] + timevals[timeind][18]
2606        graphtit = vtit.replace('_','\_') + ' (' + datev + ' ' + timev + ')'
2607
2608    plt.title(graphtit)
2609   
2610    output_kind(kfig, figname, True)
2611
2612    return
2613
2614def plot_2Dfield_easy(varv,dimxv,dimyv,dimn,colorbar,vn,vx,unit,ifile,vtit,kfig,reva):
2615    """ Adding labels and other staff to the graph
2616      varv= 2D values to plot
2617      dim[x/y]v = values at the axes of x and y
2618      dimn= dimension names to plot
2619      colorbar= name of the color bar to use
2620      vn,vm= minmum and maximum values to plot
2621      unit= units of the variable
2622      ifile= name of the input file
2623      vtit= title of the variable
2624      kfig= kind of figure (jpg, pdf, png)
2625      reva= reverse the axes (x-->y, y-->x)
2626    """
2627##    import matplotlib as mpl
2628##    mpl.use('Agg')
2629##    import matplotlib.pyplot as plt
2630    fname = 'plot_2Dfield'
2631
2632    if reva:
2633        print '  reversing the axes of the figure (x-->y, y-->x)!!'
2634        varv = np.transpose(varv)
2635        dimn0 = []
2636        dimn0.append(dimn[1] + '')
2637        dimn0.append(dimn[0] + '')
2638        dimn = dimn0
2639        if len(dimyv.shape) == 2:
2640            x = np.transpose(dimyv)
2641        else:
2642            if len(dimxv.shape) == 2:
2643                ddx = len(dimyv)
2644                ddy = dimxv.shape[1]
2645            else:
2646                ddx = len(dimyv)
2647                ddy = len(dimxv)
2648
2649            x = np.zeros((ddy,ddx), dtype=np.float)
2650            for j in range(ddy):
2651                x[j,:] = dimyv
2652
2653        if len(dimxv.shape) == 2:
2654            y = np.transpose(dimxv)
2655        else:
2656            if len(dimyv.shape) == 2:
2657                ddx = dimyv.shape[0]
2658                ddy = len(dimxv)
2659            else:
2660                ddx = len(dimyv)
2661                ddy = len(dimxv)
2662
2663            y = np.zeros((ddy,ddx), dtype=np.float)
2664            for i in range(ddx):
2665                y[:,i] = dimxv
2666    else:
2667        if len(dimxv.shape) == 2:
2668            x = dimxv
2669        else:
2670            x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
2671            for j in range(len(dimyv)):
2672                x[j,:] = dimxv
2673
2674        if len(dimyv.shape) == 2:
2675            y = dimyv
2676        else:
2677            y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
2678            for i in range(len(dimxv)):
2679                x[:,i] = dimyv
2680
2681    dx=varv.shape[1]
2682    dy=varv.shape[0]
2683
2684    plt.rc('text', usetex=True)
2685    plt.xlim(0,dx-1)
2686    plt.ylim(0,dy-1)
2687
2688    plt.pcolormesh(x, y, varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
2689#    plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
2690    cbar = plt.colorbar()
2691
2692    plt.xlabel(dimn[1].replace('_','\_'))
2693    plt.ylabel(dimn[0].replace('_','\_'))
2694
2695# set the limits of the plot to the limits of the data
2696    plt.axis([x.min(), x.max(), y.min(), y.max()])
2697#    if varv.shape[1] / varv.shape[0] > 10:
2698#        plt.axes().set_aspect(0.001)
2699#    else:
2700#        plt.axes().set_aspect(np.float(varv.shape[0])/np.float(varv.shape[1]))
2701
2702    cbar.set_label(unit)
2703
2704    figname = ifile.replace('.','_') + '_' + vtit
2705    graphtit = vtit.replace('_','\_').replace('&','\&')
2706
2707    plt.title(graphtit)
2708
2709    output_kind(kfig, figname, True)
2710   
2711    return
2712
2713def plot_Trajectories(lonval, latval, linesn, olon, olat, lonlatLims, gtit, kfig,    \
2714  mapv, obsname):
2715    """ plotting points
2716      [lon/latval]= lon,lat values to plot (as list of vectors)
2717      linesn: name of the lines
2718      o[lon/lat]= object with the longitudes and the latitudes of the map to plot
2719      lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
2720      gtit= title of the graph
2721      kfig= kind of figure (jpg, pdf, png)
2722      mapv= map characteristics: [proj],[res]
2723        see full documentation: http://matplotlib.org/basemap/
2724        [proj]: projection
2725          * 'cyl', cilindric
2726          * 'lcc', lambert conformal
2727        [res]: resolution:
2728          * 'c', crude
2729          * 'l', low
2730          * 'i', intermediate
2731          * 'h', high
2732          * 'f', full
2733      obsname= name of the observations in graph (can be None for without).
2734        Observational trajectory would be the last one
2735    """
2736    fname = 'plot_Trajectories'
2737
2738    if lonval == 'h':
2739        print fname + '_____________________________________________________________'
2740        print plot_Trajectories.__doc__
2741        quit()
2742
2743# Canging line kinds every 7 lines (end of standard colors)
2744    linekinds=['.-','x-','o-']
2745
2746    Ntraj = len(lonval)
2747
2748    if obsname is not None:
2749        Ntraj = Ntraj - 1
2750
2751    N7lines = 0
2752
2753    plt.rc('text', usetex=True)
2754
2755    if not mapv is None:
2756        if len(olon[:].shape) == 3:
2757#            lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,])
2758            lon0 = olon[0,]
2759            lat0 = olat[0,]
2760        elif len(olon[:].shape) == 2:
2761#            lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
2762            lon0 = olon[:]
2763            lat0 = olat[:]
2764        elif len(olon[:].shape) == 1:
2765#            lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
2766            lon00 = olon[:]
2767            lat00 = olat[:]
2768            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2769            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2770
2771            for iy in range(len(lat00)):
2772                lon0[iy,:] = lon00
2773            for ix in range(len(lon00)):
2774                lat0[:,ix] = lat00
2775
2776        map_proj=mapv.split(',')[0]
2777        map_res=mapv.split(',')[1]
2778
2779        dx = lon0.shape[1]
2780        dy = lon0.shape[0]
2781
2782        nlon = lon0[0,0]
2783        xlon = lon0[dy-1,dx-1]
2784        nlat = lat0[0,0]
2785        xlat = lat0[dy-1,dx-1]
2786
2787        lon2 = lon0[dy/2,dx/2]
2788        lat2 = lat0[dy/2,dx/2]
2789
2790        if lonlatLims is not None:
2791            plt.xlim(lonlatLims[0], lonlatLims[2])
2792            plt.ylim(lonlatLims[1], lonlatLims[3])
2793            if map_proj == 'cyl':
2794                nlon = lonlatLims[0]
2795                nlat = lonlatLims[1]
2796                xlon = lonlatLims[2]
2797                xlat = lonlatLims[3]
2798
2799        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
2800          xlon, ',', xlat
2801
2802        if map_proj == 'cyl':
2803            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
2804              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2805        elif map_proj == 'lcc':
2806            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
2807              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2808
2809        if len(olon.shape) == 3:
2810#            lons, lats = np.meshgrid(olon[0,:,:], olat[0,:,:])
2811            lons = olon[0,:,:]
2812            lats = olat[0,:,:]
2813
2814        elif len(olon.shape) == 2:
2815#            lons, lats = np.meshgrid(olon[:,:], olat[:,:])
2816            lons = olon[:,:]
2817            lats = olat[:,:]
2818        else:
2819            dx = olon.shape
2820            dy = olat.shape
2821#            print errormsg
2822#            print '  ' + fname + ': shapes of lon/lat objects', olon.shape,          \
2823#              'not ready!!!'
2824
2825        for il in range(Ntraj):
2826            plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il])
2827            if il == 6: N7lines = N7lines + 1
2828
2829        m.drawcoastlines()
2830
2831        meridians = pretty_int(nlon,xlon,5)
2832        m.drawmeridians(meridians,labels=[True,False,False,True])
2833
2834        parallels = pretty_int(nlat,xlat,5)
2835        m.drawparallels(parallels,labels=[False,True,True,False])
2836
2837        plt.xlabel('W-E')
2838        plt.ylabel('S-N')
2839
2840    else:
2841        if len(olon.shape) == 3:
2842            dx = olon.shape[2]
2843            dy = olon.shape[1]
2844        elif len(olon.shape) == 2:
2845            dx = olon.shape[1]
2846            dy = olon.shape[0]
2847        else:
2848            dx = olon.shape
2849            dy = olat.shape
2850#            print errormsg
2851#            print '  ' + fname + ': shapes of lon/lat objects', olon.shape,          \
2852#              'not ready!!!'
2853
2854        if lonlatLims is not None:
2855            plt.xlim(lonlatLims[0], lonlatLims[2])
2856            plt.ylim(lonlatLims[1], lonlatLims[3])
2857        else:
2858            plt.xlim(np.min(olon[:]),np.max(olon[:]))
2859            plt.ylim(np.min(olat[:]),np.max(olat[:]))
2860
2861        for il in range(Ntraj):
2862            plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il])
2863            if il == 6: N7lines = N7lines + 1
2864
2865        plt.xlabel('x-axis')
2866        plt.ylabel('y-axis')
2867
2868    figname = 'trajectories'
2869    graphtit = gtit
2870
2871    if obsname is not None:
2872        plt.plot(lonval[Ntraj], latval[Ntraj], linestyle='-', color='k',             \
2873          linewidth=3, label= obsname)
2874
2875    plt.title(graphtit.replace('_','\_').replace('&','\&'))
2876    plt.legend()
2877   
2878    output_kind(kfig, figname, True)
2879
2880    return
2881
2882def plot_topo_geogrid(varv, olon, olat, mint, maxt, lonlatLims, gtit, kfig, mapv,    \
2883  closeif):
2884    """ plotting geo_em.d[nn].nc topography from WPS files
2885    plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv)
2886      varv= topography values
2887      o[lon/lat]= longitude and latitude objects
2888      [min/max]t: minimum and maximum values of topography to draw
2889      lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
2890      gtit= title of the graph
2891      kfig= kind of figure (jpg, pdf, png)
2892      mapv= map characteristics: [proj],[res]
2893        see full documentation: http://matplotlib.org/basemap/
2894        [proj]: projection
2895          * 'cyl', cilindric
2896          * 'lcc', lamvbert conformal
2897        [res]: resolution:
2898          * 'c', crude
2899          * 'l', low
2900          * 'i', intermediate
2901          * 'h', high
2902          * 'f', full
2903      closeif= Boolean value if the figure has to be closed
2904    """
2905    fname = 'plot_topo_geogrid'
2906
2907    if varv == 'h':
2908        print fname + '_____________________________________________________________'
2909        print plot_topo_geogrid.__doc__
2910        quit()
2911
2912    dx=varv.shape[1]
2913    dy=varv.shape[0]
2914
2915    plt.rc('text', usetex=True)
2916#    plt.rc('font', family='serif')
2917
2918    if not mapv is None:
2919        if len(olon[:].shape) == 3:
2920            lon0 = olon[0,]
2921            lat0 = olat[0,]
2922        elif len(olon[:].shape) == 2:
2923            lon0 = olon[:]
2924            lat0 = olat[:]
2925        elif len(olon[:].shape) == 1:
2926            lon00 = olon[:]
2927            lat00 = olat[:]
2928            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2929            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
2930
2931            for iy in range(len(lat00)):
2932                lon0[iy,:] = lon00
2933            for ix in range(len(lon00)):
2934                lat0[:,ix] = lat00
2935
2936        map_proj=mapv.split(',')[0]
2937        map_res=mapv.split(',')[1]
2938        dx = lon0.shape[1]
2939        dy = lon0.shape[0]
2940
2941        if lonlatLims is not None:
2942            print '  ' + fname + ': cutting the domain to plot !!!!'
2943            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
2944            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
2945            nlon =  lonlatLims[0]
2946            xlon =  lonlatLims[2]
2947            nlat =  lonlatLims[1]
2948            xlat =  lonlatLims[3]
2949
2950            if map_proj == 'lcc':
2951                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
2952                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
2953        else:
2954            nlon = lon0[0,0]
2955            xlon = lon0[dy-1,dx-1]
2956            nlat = lat0[0,0]
2957            xlat = lat0[dy-1,dx-1]
2958            lon2 = lon0[dy/2,dx/2]
2959            lat2 = lat0[dy/2,dx/2]
2960
2961        plt.xlim(nlon, xlon)
2962        plt.ylim(nlat, xlat)
2963        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
2964          xlon, ',', xlat
2965
2966        if map_proj == 'cyl':
2967            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
2968              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2969        elif map_proj == 'lcc':
2970            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
2971              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
2972        else:
2973            print errormsg
2974            print '  ' + fname + ": map projection '" + map_proj + "' not ready !!"
2975            quit(-1)
2976
2977        if len(olon[:].shape) == 1:
2978            lons, lats = np.meshgrid(olon[:], olat[:])
2979        else:
2980            if len(olon[:].shape) == 3:
2981                lons = olon[0,:,:]
2982                lats = olat[0,:,:]
2983            else:
2984                lons = olon[:]
2985                lats = olat[:]
2986 
2987        x,y = m(lons,lats)
2988
2989        plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt)
2990        cbar = plt.colorbar()
2991
2992        m.drawcoastlines()
2993
2994        meridians = pretty_int(nlon,xlon,5)
2995        m.drawmeridians(meridians,labels=[True,False,False,True])
2996
2997        parallels = pretty_int(nlat,xlat,5)
2998        m.drawparallels(parallels,labels=[False,True,True,False])
2999
3000        plt.xlabel('W-E')
3001        plt.ylabel('S-N')
3002    else:
3003        print emsg
3004        print '  ' + fname + ': A projection parameter is needed None given !!'
3005        quit(-1)   
3006
3007    figname = 'domain'
3008    graphtit = gtit.replace('_','\_')
3009    cbar.set_label('height ($m$)')
3010
3011    plt.title(graphtit.replace('_','\_').replace('&','\&'))
3012   
3013    output_kind(kfig, figname, closeif)
3014
3015    return
3016
3017def plot_topo_geogrid_boxes(varv, boxesX, boxesY, boxlabels, olon, olat, mint, maxt, \
3018  lonlatLims, gtit, kfig, mapv, closeif):
3019    """ plotting geo_em.d[nn].nc topography from WPS files
3020    plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv)
3021      varv= topography values
3022      boxesX/Y= 4-line sets to draw the boxes
3023      boxlabels= labels for the legend of the boxes
3024      o[lon/lat]= longitude and latitude objects
3025      [min/max]t: minimum and maximum values of topography to draw
3026      lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
3027      gtit= title of the graph
3028      kfig= kind of figure (jpg, pdf, png)
3029      mapv= map characteristics: [proj],[res]
3030        see full documentation: http://matplotlib.org/basemap/
3031        [proj]: projection
3032          * 'cyl', cilindric
3033          * 'lcc', lamvbert conformal
3034        [res]: resolution:
3035          * 'c', crude
3036          * 'l', low
3037          * 'i', intermediate
3038          * 'h', high
3039          * 'f', full
3040      closeif= Boolean value if the figure has to be closed
3041    """
3042    fname = 'plot_topo_geogrid'
3043
3044    if varv == 'h':
3045        print fname + '_____________________________________________________________'
3046        print plot_topo_geogrid.__doc__
3047        quit()
3048
3049    cols = color_lines(len(boxlabels))
3050
3051    dx=varv.shape[1]
3052    dy=varv.shape[0]
3053
3054    plt.rc('text', usetex=True)
3055#    plt.rc('font', family='serif')
3056
3057    if not mapv is None:
3058        if len(olon[:].shape) == 3:
3059            lon0 = olon[0,]
3060            lat0 = olat[0,]
3061        elif len(olon[:].shape) == 2:
3062            lon0 = olon[:]
3063            lat0 = olat[:]
3064        elif len(olon[:].shape) == 1:
3065            lon00 = olon[:]
3066            lat00 = olat[:]
3067            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
3068            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
3069
3070            for iy in range(len(lat00)):
3071                lon0[iy,:] = lon00
3072            for ix in range(len(lon00)):
3073                lat0[:,ix] = lat00
3074
3075        map_proj=mapv.split(',')[0]
3076        map_res=mapv.split(',')[1]
3077        dx = lon0.shape[1]
3078        dy = lon0.shape[0]
3079
3080        if lonlatLims is not None:
3081            print '  ' + fname + ': cutting the domain to plot !!!!'
3082            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
3083            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
3084            nlon =  lonlatLims[0]
3085            xlon =  lonlatLims[2]
3086            nlat =  lonlatLims[1]
3087            xlat =  lonlatLims[3]
3088
3089            if map_proj == 'lcc':
3090                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
3091                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
3092        else:
3093            nlon = np.min(lon0)
3094            xlon = np.max(lon0)
3095            nlat = np.min(lat0)
3096            xlat = np.max(lat0)
3097            lon2 = lon0[dy/2,dx/2]
3098            lat2 = lat0[dy/2,dx/2]
3099
3100        plt.xlim(nlon, xlon)
3101        plt.ylim(nlat, xlat)
3102        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
3103          xlon, ',', xlat
3104
3105        if map_proj == 'cyl':
3106            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
3107              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
3108        elif map_proj == 'lcc':
3109            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
3110              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
3111
3112        if len(olon[:].shape) == 1:
3113            lons, lats = np.meshgrid(olon[:], olat[:])
3114        else:
3115            if len(olon[:].shape) == 3:
3116                lons = olon[0,:,:]
3117                lats = olat[0,:,:]
3118            else:
3119                lons = olon[:]
3120                lats = olat[:]
3121 
3122        x,y = m(lons,lats)
3123
3124        plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt)
3125        cbar = plt.colorbar()
3126
3127        Nboxes = len(boxesX)/4
3128        for ibox in range(Nboxes):
3129            plt.plot(boxesX[ibox*4], boxesY[ibox*4], linestyle='-', linewidth=3,     \
3130              label=boxlabels[ibox], color=cols[ibox])
3131            plt.plot(boxesX[ibox*4+1], boxesY[ibox*4+1], linestyle='-', linewidth=3, \
3132              color=cols[ibox])
3133            plt.plot(boxesX[ibox*4+2], boxesY[ibox*4+2], linestyle='-', linewidth=3, \
3134              color=cols[ibox])
3135            plt.plot(boxesX[ibox*4+3], boxesY[ibox*4+3], linestyle='-', linewidth=3, \
3136              color=cols[ibox])
3137
3138        m.drawcoastlines()
3139
3140        meridians = pretty_int(nlon,xlon,5)
3141        m.drawmeridians(meridians,labels=[True,False,False,True])
3142
3143        parallels = pretty_int(nlat,xlat,5)
3144        m.drawparallels(parallels,labels=[False,True,True,False])
3145
3146        plt.xlabel('W-E')
3147        plt.ylabel('S-N')
3148    else:
3149        print emsg
3150        print '  ' + fname + ': A projection parameter is needed None given !!'
3151        quit(-1)   
3152
3153    figname = 'domain_boxes'
3154    graphtit = gtit.replace('_','\_').replace('&','\&')
3155    cbar.set_label('height ($m$)')
3156
3157    plt.title(graphtit)
3158    plt.legend()
3159   
3160    output_kind(kfig, figname, closeif)
3161
3162    return
3163
3164def plot_2D_shadow(varsv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,          \
3165  colorbar,vs,uts,vtit,kfig,reva,mapv,ifclose):
3166    """ Adding labels and other staff to the graph
3167      varsv= 2D values to plot with shading
3168      vnames= variable names for the figure
3169      dim[x/y]v = values at the axes of x and y
3170      dim[x/y]u = units at the axes of x and y
3171      dimn= dimension names to plot
3172      colorbar= name of the color bar to use
3173      vs= minmum and maximum values to plot in shadow or:
3174        'Srange': for full range
3175        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
3176        'Saroundminmax@val': for min*val,max*val
3177        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
3178          percentile_(100-val)-median)
3179        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
3180        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
3181        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
3182           percentile_(100-val)-median)
3183      uts= units of the variable to shadow
3184      vtit= title of the variable
3185      kfig= kind of figure (jpg, pdf, png)
3186      reva= ('|' for combination)
3187        * 'transpose': reverse the axes (x-->y, y-->x)
3188        * 'flip'@[x/y]: flip the axis x or y
3189      mapv= map characteristics: [proj],[res]
3190        see full documentation: http://matplotlib.org/basemap/
3191        [proj]: projection
3192          * 'cyl', cilindric
3193          * 'lcc', lambert conformal
3194        [res]: resolution:
3195          * 'c', crude
3196          * 'l', low
3197          * 'i', intermediate
3198          * 'h', high
3199          * 'f', full
3200      ifclose= boolean value whether figure should be close (finish) or not
3201    """
3202##    import matplotlib as mpl
3203##    mpl.use('Agg')
3204##    import matplotlib.pyplot as plt
3205    fname = 'plot_2D_shadow'
3206
3207#    print dimyv[73,21]
3208#    dimyv[73,21] = -dimyv[73,21]
3209#    print 'Lluis dimsv: ',np.min(dimxv), np.max(dimxv), ':', np.min(dimyv), np.max(dimyv)
3210
3211    if varsv == 'h':
3212        print fname + '_____________________________________________________________'
3213        print plot_2D_shadow.__doc__
3214        quit()
3215
3216    if len(varsv.shape) != 2:
3217        print errormsg
3218        print '  ' + fname + ': wrong variable shape:',varv.shape,'is has to be 2D!!'
3219        quit(-1)
3220
3221    reva0 = ''
3222    if reva.find('|') != 0:
3223        revas = reva.split('|')
3224    else:
3225        revas = [reva]
3226        reva0 = reva
3227
3228    for rev in revas:
3229        if reva[0:4] == 'flip':
3230            reva0 = 'flip'
3231            if len(reva.split('@')) != 2:
3232                 print errormsg
3233                 print '  ' + fname + ': flip is given', reva, 'but not axis!'
3234                 quit(-1)
3235
3236        if rev == 'transpose':
3237            print '  reversing the axes of the figure (x-->y, y-->x)!!'
3238            varsv = np.transpose(varsv)
3239            dxv = dimyv
3240            dyv = dimxv
3241            dimxv = dxv
3242            dimyv = dyv
3243
3244    if len(dimxv[:].shape) == 3:
3245        xdims = '1,2'
3246    elif len(dimxv[:].shape) == 2:
3247        xdims = '0,1'
3248    elif len(dimxv[:].shape) == 1:
3249        xdims = '0'
3250
3251    if len(dimyv[:].shape) == 3:
3252        ydims = '1,2'
3253    elif len(dimyv[:].shape) == 2:
3254        ydims = '0,1'
3255    elif len(dimyv[:].shape) == 1:
3256        ydims = '0'
3257
3258    lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)
3259
3260    if not mapv is None:
3261        map_proj=mapv.split(',')[0]
3262        map_res=mapv.split(',')[1]
3263
3264        dx = lon0.shape[1]
3265        dy = lat0.shape[0]
3266
3267        nlon = lon0[0,0]
3268        xlon = lon0[dy-1,dx-1]
3269        nlat = lat0[0,0]
3270        xlat = lat0[dy-1,dx-1]
3271
3272# Thats too much! :)
3273#        if lonlatLims is not None:
3274#            print '  ' + fname + ': cutting the domain to plot !!!!'
3275#            plt.xlim(lonlatLims[0], lonlatLims[2])
3276#            plt.ylim(lonlatLims[1], lonlatLims[3])
3277#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
3278#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
3279
3280#            if map_proj == 'cyl':
3281#                nlon = lonlatLims[0]
3282#                nlat = lonlatLims[1]
3283#                xlon = lonlatLims[2]
3284#                xlat = lonlatLims[3]
3285#            elif map_proj == 'lcc':
3286#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
3287#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
3288#                nlon =  lonlatLims[0]
3289#                xlon =  lonlatLims[2]
3290#                nlat =  lonlatLims[1]
3291#                xlat =  lonlatLims[3]
3292
3293        lon2 = lon0[dy/2,dx/2]
3294        lat2 = lat0[dy/2,dx/2]
3295
3296        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
3297          xlon, ',', xlat
3298
3299        if map_proj == 'cyl':
3300            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
3301              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
3302        elif map_proj == 'lcc':
3303            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
3304              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
3305        else:
3306            print errormsg
3307            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
3308            print '    available: cyl, lcc'
3309            quit(-1)
3310 
3311        x,y = m(lon0,lat0)
3312
3313    else:
3314        x = lon0
3315        y = lat0
3316
3317    vsend = np.zeros((2), dtype=np.float)
3318# Changing limits of the colors
3319    if type(vs[0]) != type(np.float(1.)):
3320        if vs[0] == 'Srange':
3321            vsend[0] = np.min(varsv)
3322        elif vs[0][0:11] == 'Saroundmean':
3323            meanv = np.mean(varsv)
3324            permean = np.float(vs[0].split('@')[1])
3325            minv = np.min(varsv)*permean
3326            maxv = np.max(varsv)*permean
3327            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
3328            vsend[0] = meanv-minextrm
3329            vsend[1] = meanv+minextrm
3330        elif vs[0][0:13] == 'Saroundminmax':
3331            permean = np.float(vs[0].split('@')[1])
3332            minv = np.min(varsv)*permean
3333            maxv = np.max(varsv)*permean
3334            vsend[0] = minv
3335            vsend[1] = maxv
3336        elif vs[0][0:17] == 'Saroundpercentile':
3337            medianv = np.median(varsv)
3338            valper = np.float(vs[0].split('@')[1])
3339            minv = np.percentile(varsv, valper)
3340            maxv = np.percentile(varsv, 100.-valper)
3341            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
3342            vsend[0] = medianv-minextrm
3343            vsend[1] = medianv+minextrm
3344        elif vs[0][0:5] == 'Smean':
3345            meanv = np.mean(varsv)
3346            permean = np.float(vs[0].split('@')[1])
3347            minv = np.min(varsv)*permean
3348            maxv = np.max(varsv)*permean
3349            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
3350            vsend[0] = -minextrm
3351            vsend[1] = minextrm
3352        elif vs[0][0:7] == 'Smedian':
3353            medianv = np.median(varsv)
3354            permedian = np.float(vs[0].split('@')[1])
3355            minv = np.min(varsv)*permedian
3356            maxv = np.max(varsv)*permedian
3357            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
3358            vsend[0] = -minextrm
3359            vsend[1] = minextrm
3360        elif vs[0][0:11] == 'Spercentile':
3361            medianv = np.median(varsv)
3362            valper = np.float(vs[0].split('@')[1])
3363            minv = np.percentile(varsv, valper)
3364            maxv = np.percentile(varsv, 100.-valper)
3365            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
3366            vsend[0] = -minextrm
3367            vsend[1] = minextrm
3368        else:
3369            print errormsg
3370            print '  ' + fname + ": range '" + vs[0] + "' not ready!!!"
3371            quit(-1)
3372        print '    ' + fname + ': modified shadow min,max:',vsend
3373    else:
3374        vsend[0] = vs[0]
3375
3376    if type(vs[0]) != type(np.float(1.)):
3377        if vs[1] == 'range':
3378            vsend[1] = np.max(varsv)
3379    else:
3380        vsend[1] = vs[1]
3381
3382    plt.rc('text', usetex=True)
3383
3384    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
3385    cbar = plt.colorbar()
3386
3387    if not mapv is None:
3388        m.drawcoastlines()
3389
3390        meridians = pretty_int(nlon,xlon,5)
3391        m.drawmeridians(meridians,labels=[True,False,False,True])
3392        parallels = pretty_int(nlat,xlat,5)
3393        m.drawparallels(parallels,labels=[False,True,True,False])
3394
3395        plt.xlabel('W-E')
3396        plt.ylabel('S-N')
3397    else:
3398        plt.xlabel(variables_values(dimn[1])[0].replace('_','\_') + ' (' +           \
3399          units_lunits(dimxu) + ')')
3400        plt.ylabel(variables_values(dimn[0])[0].replace('_','\_') + ' (' +           \
3401          units_lunits(dimyu) + ')')
3402
3403    txpos = pretty_int(x.min(),x.max(),5)
3404    typos = pretty_int(y.min(),y.max(),5)
3405    txlabels = list(txpos)
3406    for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
3407    tylabels = list(typos)
3408    for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
3409
3410# set the limits of the plot to the limits of the data
3411
3412    if searchInlist(revas,'transpose'):
3413        x0 = y
3414        y0 = x
3415        x = x0
3416        y = y0
3417#    print 'Lluis reva0:',reva0,'x min,max:',x.min(),x.max(),' y min,max:',y.min(),y.max()
3418
3419    if reva0 == 'flip':
3420        if reva.split('@')[1] == 'x':
3421            plt.axis([x.max(), x.min(), y.min(), y.max()])
3422        elif reva.split('@')[1] == 'y':
3423            plt.axis([x.min(), x.max(), y.max(), y.min()])
3424        else:
3425            plt.axis([x.max(), x.min(), y.max(), y.min()])
3426    else:
3427        plt.axis([x.min(), x.max(), y.min(), y.max()])
3428
3429    if mapv is None:
3430        plt.xticks(txpos, txlabels)
3431        plt.yticks(typos, tylabels)
3432
3433# units labels
3434    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
3435
3436    figname = '2Dfields_shadow'
3437    graphtit = vtit.replace('_','\_').replace('&','\&')
3438
3439    plt.title(graphtit)
3440   
3441    output_kind(kfig, figname, ifclose)
3442
3443    return
3444
3445#Nvals=50
3446#vals1 = np.zeros((Nvals,Nvals), dtype= np.float)
3447#for j in range(Nvals):
3448#    for i in range(Nvals):
3449#      vals1[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.)
3450
3451#plot_2D_shadow(vals1, 'var1', np.arange(50)*1., np.arange(50)*1., 'ms-1',            \
3452#  'm', ['lat','lon'], 'rainbow', [0, Nvals], 'ms-1', 'test var1', 'pdf', 'None',     \
3453#  None, True)
3454#quit()
3455
3456def plot_2D_shadow_time(varsv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,colorbar,vs,uts,   \
3457  vtit,kfig,reva,taxis,tpos,tlabs,ifclose):
3458    """ Plotting a 2D field with one of the axes being time
3459      varsv= 2D values to plot with shading
3460      vnames= shading variable name for the figure
3461      dim[x/y]v= values at the axes of x and y
3462      dim[x/y]u= units at the axes of x and y
3463      dimn= dimension names to plot
3464      colorbar= name of the color bar to use
3465      vs= minmum and maximum values to plot in shadow or:
3466        'Srange': for full range
3467        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
3468        'Saroundminmax@val': for min*val,max*val
3469        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
3470          percentile_(100-val)-median)
3471        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
3472        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
3473        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
3474           percentile_(100-val)-median)
3475      uts= units of the variable to shadow
3476      vtit= title of the variable
3477      kfig= kind of figure (jpg, pdf, png)
3478      reva=
3479        * 'transpose': reverse the axes (x-->y, y-->x)
3480        * 'flip'@[x/y]: flip the axis x or y
3481      taxis= Which is the time-axis
3482      tpos= positions of the time ticks
3483      tlabs= labels of the time ticks
3484      ifclose= boolean value whether figure should be close (finish) or not
3485    """
3486    fname = 'plot_2D_shadow_time'
3487
3488    if varsv == 'h':
3489        print fname + '_____________________________________________________________'
3490        print plot_2D_shadow_time.__doc__
3491        quit()
3492
3493# Definning ticks labels
3494    if taxis == 'x':
3495        txpos = tpos
3496        txlabels = tlabs
3497        plxlabel = dimxu
3498        typos = pretty_int(np.min(dimyv),np.max(dimyv),10)
3499        tylabels = list(typos)
3500        for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
3501        plylabel = variables_values(dimn[0])[0].replace('_','\_') + ' (' +           \
3502          units_lunits(dimyu) + ')'
3503    else:
3504        txpos = pretty_int(np.min(dimxv),np.max(dimxv),10)
3505        txlabels = list(txpos)
3506        plxlabel = variables_values(dimn[1])[0].replace('_','\_') + ' (' +           \
3507          units_lunits(dimxu) + ')'
3508        typos = tpos
3509        tylabels = tlabs
3510        plylabel = dimyu
3511
3512# Transposing/flipping axis
3513    if reva.find('|') != 0:
3514        revas = reva.split('|')
3515        reva0 = ''
3516    else:
3517        revas = [reva]
3518        reva0 = reva
3519
3520    for rev in revas:
3521        if rev[0:4] == 'flip':
3522            reva0 = 'flip'
3523            if len(reva.split('@')) != 2:
3524                 print errormsg
3525                 print '  ' + fname + ': flip is given', reva, 'but not axis!'
3526                 quit(-1)
3527            else:
3528                 print "  flipping '" + rev.split('@')[1] + "' axis !"
3529
3530        if rev == 'transpose':
3531            print '  reversing the axes of the figure (x-->y, y-->x)!!'
3532# Flipping values of variable
3533            varsv = np.transpose(varsv)
3534            dxv = dimyv
3535            dyv = dimxv
3536            dimxv = dxv
3537            dimyv = dyv
3538
3539    if len(dimxv.shape) == 3:
3540        dxget='1,2'
3541    elif len(dimxv.shape) == 2:
3542        dxget='0,1'
3543    elif len(dimxv.shape) == 1:
3544        dxget='0'
3545    else:
3546        print errormsg
3547        print '  ' + fname + ': shape of x-values:',dimxv.shape,'not ready!!'
3548        quit(-1)
3549
3550    if len(dimyv.shape) == 3:
3551        dyget='1,2'
3552    elif len(dimyv.shape) == 2:
3553        dyget='0,1'
3554    elif len(dimyv.shape) == 1:
3555        dyget='0'
3556    else:
3557        print errormsg
3558        print '  ' + fname + ': shape of y-values:',dimyv.shape,'not ready!!'
3559        quit(-1)
3560
3561    x,y = dxdy_lonlat(dimxv,dimyv,dxget,dyget)
3562
3563    plt.rc('text', usetex=True)
3564
3565    vsend = np.zeros((2), dtype=np.float)
3566# Changing limits of the colors
3567    if type(vs[0]) != type(np.float(1.)):
3568        if vs[0] == 'Srange':
3569            vsend[0] = np.min(varsv)
3570        elif vs[0][0:11] == 'Saroundmean':
3571            meanv = np.mean(varsv)
3572            permean = np.float(vs[0].split('@')[1])
3573            minv = np.min(varsv)*permean
3574            maxv = np.max(varsv)*permean
3575            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
3576            vsend[0] = meanv-minextrm
3577            vsend[1] = meanv+minextrm
3578        elif vs[0][0:13] == 'Saroundminmax':
3579            permean = np.float(vs[0].split('@')[1])
3580            minv = np.min(varsv)*permean
3581            maxv = np.max(varsv)*permean
3582            vsend[0] = minv
3583            vsend[1] = maxv
3584        elif vs[0][0:17] == 'Saroundpercentile':
3585            medianv = np.median(varsv)
3586            valper = np.float(vs[0].split('@')[1])
3587            minv = np.percentile(varsv, valper)
3588            maxv = np.percentile(varsv, 100.-valper)
3589            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
3590            vsend[0] = medianv-minextrm
3591            vsend[1] = medianv+minextrm
3592        elif vs[0][0:5] == 'Smean':
3593            meanv = np.mean(varsv)
3594            permean = np.float(vs[0].split('@')[1])
3595            minv = np.min(varsv)*permean
3596            maxv = np.max(varsv)*permean
3597            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
3598            vsend[0] = -minextrm
3599            vsend[1] = minextrm
3600        elif vs[0][0:7] == 'Smedian':
3601            medianv = np.median(varsv)
3602            permedian = np.float(vs[0].split('@')[1])
3603            minv = np.min(varsv)*permedian
3604            maxv = np.max(varsv)*permedian
3605            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
3606            vsend[0] = -minextrm
3607            vsend[1] = minextrm
3608        elif vs[0][0:11] == 'Spercentile':
3609            medianv = np.median(varsv)
3610            valper = np.float(vs[0].split('@')[1])
3611            minv = np.percentile(varsv, valper)
3612            maxv = np.percentile(varsv, 100.-valper)
3613            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
3614            vsend[0] = -minextrm
3615            vsend[1] = minextrm
3616        else:
3617            print errormsg
3618            print '  ' + fname + ": range '" + vs[0] + "' not ready!!!"
3619            quit(-1)
3620        print '    ' + fname + ': modified shadow min,max:',vsend
3621    else:
3622        vsend[0] = vs[0]
3623
3624    if type(vs[0]) != type(np.float(1.)):
3625        if vs[1] == 'range':
3626            vsend[1] = np.max(varsv)
3627    else:
3628        vsend[1] = vs[1]
3629
3630    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
3631    cbar = plt.colorbar()
3632
3633#    print 'Lluis reva0:',reva0,'x min,max:',x.min(),x.max(),' y min,max:',y.min(),y.max()
3634
3635# set the limits of the plot to the limits of the data
3636    if reva0 == 'flip':
3637        if reva.split('@')[1] == 'x':
3638            plt.axis([x.max(), x.min(), y.min(), y.max()])
3639        elif reva.split('@')[1] == 'y':
3640            plt.axis([x.min(), x.max(), y.max(), y.min()])
3641        else:
3642            plt.axis([x.max(), x.min(), y.max(), y.min()])
3643    else:
3644        plt.axis([x.min(), x.max(), y.min(), y.max()])
3645
3646    if searchInlist(revas, 'transpose'):
3647        plt.xticks(typos, tylabels)
3648        plt.yticks(txpos, txlabels)
3649        plt.xlabel(plylabel)
3650        plt.ylabel(plxlabel)
3651    else:
3652        plt.xticks(txpos, txlabels)
3653        plt.yticks(typos, tylabels)
3654        plt.xlabel(plxlabel)
3655        plt.ylabel(plylabel)
3656
3657# units labels
3658    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
3659
3660    figname = '2Dfields_shadow_time'
3661    graphtit = vtit.replace('_','\_').replace('&','\&')
3662
3663    plt.title(graphtit)
3664   
3665    output_kind(kfig, figname, ifclose)
3666
3667    return
3668
3669def plot_2D_shadow_contour(varsv,varcv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,          \
3670  colorbar,ckind,clabfmt,vs,vc,uts,vtit,kfig,reva,mapv):
3671    """ Adding labels and other staff to the graph
3672      varsv= 2D values to plot with shading
3673      varcv= 2D values to plot with contours
3674      vnames= variable names for the figure
3675      dim[x/y]v = values at the axes of x and y
3676      dim[x/y]u = units at the axes of x and y
3677      dimn= dimension names to plot
3678      colorbar= name of the color bar to use
3679      ckind= contour kind
3680        'cmap': as it gets from colorbar
3681        'fixc,[colname]': fixed color [colname], all stright lines
3682        'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
3683      clabfmt= format of the labels in the contour plot (None, no labels)
3684      vs= minmum and maximum values to plot in shadow
3685      vc= vector with the levels for the contour
3686      uts= units of the variable [u-shadow, u-contour]
3687      vtit= title of the variable
3688      kfig= kind of figure (jpg, pdf, png)
3689      reva=
3690        * 'transpose': reverse the axes (x-->y, y-->x)
3691        * 'flip'@[x/y]: flip the axis x or y
3692      mapv= map characteristics: [proj],[res]
3693        see full documentation: http://matplotlib.org/basemap/
3694        [proj]: projection
3695          * 'cyl', cilindric
3696          * 'lcc', lamvbert conformal
3697        [res]: resolution:
3698          * 'c', crude
3699          * 'l', low
3700          * 'i', intermediate
3701          * 'h', high
3702          * 'f', full
3703    """
3704##    import matplotlib as mpl
3705##    mpl.use('Agg')
3706##    import matplotlib.pyplot as plt
3707    fname = 'plot_2D_shadow_contour'
3708
3709
3710    if varsv == 'h':
3711        print fname + '_____________________________________________________________'
3712        print plot_2D_shadow_contour.__doc__
3713        quit()
3714
3715    if reva[0:4] == 'flip':
3716        reva0 = 'flip'
3717        if len(reva.split('@')) != 2:
3718             print errormsg
3719             print '  ' + fname + ': flip is given', reva, 'but not axis!'
3720             quit(-1)
3721    else:
3722        reva0 = reva
3723
3724    if reva0 == 'transpose':
3725        print '  reversing the axes of the figure (x-->y, y-->x)!!'
3726        varsv = np.transpose(varsv)
3727        varcv = np.transpose(varcv)
3728        dxv = dimyv
3729        dyv = dimxv
3730        dimxv = dxv
3731        dimyv = dyv
3732
3733    if not mapv is None:
3734        if len(dimxv[:].shape) == 3:
3735            lon0 = dimxv[0,]
3736            lat0 = dimyv[0,]
3737        elif len(dimxv[:].shape) == 2:
3738            lon0 = dimxv[:]
3739            lat0 = dimyv[:]
3740        elif len(dimxv[:].shape) == 1:
3741            lon00 = dimxv[:]
3742            lat00 = dimyv[:]
3743            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
3744            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
3745
3746            for iy in range(len(lat00)):
3747                lon0[iy,:] = lon00
3748            for ix in range(len(lon00)):
3749                lat0[:,ix] = lat00
3750
3751        map_proj=mapv.split(',')[0]
3752        map_res=mapv.split(',')[1]
3753
3754        dx = lon0.shape[1]
3755        dy = lon0.shape[0]
3756
3757        nlon = lon0[0,0]
3758        xlon = lon0[dy-1,dx-1]
3759        nlat = lat0[0,0]
3760        xlat = lat0[dy-1,dx-1]
3761
3762# Thats too much! :)
3763#        if lonlatLims is not None:
3764#            print '  ' + fname + ': cutting the domain to plot !!!!'
3765#            plt.xlim(lonlatLims[0], lonlatLims[2])
3766#            plt.ylim(lonlatLims[1], lonlatLims[3])
3767#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
3768#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
3769
3770#            if map_proj == 'cyl':
3771#                nlon = lonlatLims[0]
3772#                nlat = lonlatLims[1]
3773#                xlon = lonlatLims[2]
3774#                xlat = lonlatLims[3]
3775#            elif map_proj == 'lcc':
3776#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
3777#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
3778#                nlon =  lonlatLims[0]
3779#                xlon =  lonlatLims[2]
3780#                nlat =  lonlatLims[1]
3781#                xlat =  lonlatLims[3]
3782
3783        lon2 = lon0[dy/2,dx/2]
3784        lat2 = lat0[dy/2,dx/2]
3785
3786        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
3787          xlon, ',', xlat
3788
3789        if map_proj == 'cyl':
3790            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
3791              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
3792        elif map_proj == 'lcc':
3793            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
3794              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
3795
3796        if len(dimxv.shape) == 1:
3797            lons, lats = np.meshgrid(dimxv, dimyv)
3798        else:
3799            if len(dimxv.shape) == 3:
3800                lons = dimxv[0,:,:]
3801                lats = dimyv[0,:,:]
3802            else:
3803                lons = dimxv[:]
3804                lats = dimyv[:]
3805 
3806        x,y = m(lons,lats)
3807
3808    else:
3809        if len(dimxv.shape) == 2:
3810            x = dimxv
3811        else:
3812            if len(dimyv.shape) == 1:
3813                x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
3814                for j in range(len(dimyv)):
3815                    x[j,:] = dimxv
3816            else:
3817                x = np.zeros((dimyv.shape), dtype=np.float)
3818                if x.shape[0] == dimxv.shape[0]:
3819                    for j in range(x.shape[1]):
3820                        x[:,j] = dimxv
3821                else:
3822                    for j in range(x.shape[0]):
3823                        x[j,:] = dimxv
3824
3825        if len(dimyv.shape) == 2:
3826            y = dimyv
3827        else:
3828            if len(dimxv.shape) == 1:
3829                y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
3830                for i in range(len(dimxv)):
3831                    y[:,i] = dimyv
3832            else:
3833                y = np.zeros((dimxv.shape), dtype=np.float)
3834
3835                if y.shape[0] == dimyv.shape[0]:
3836                    for i in range(y.shape[1]):
3837                        y[i,:] = dimyv
3838                else:
3839                    for i in range(y.shape[0]):
3840                        y[i,:] = dimyv
3841
3842    plt.rc('text', usetex=True)
3843
3844    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
3845    cbar = plt.colorbar()
3846
3847# contour
3848##
3849    contkind = ckind.split(',')[0]
3850    if contkind == 'cmap':
3851        cplot = plt.contour(x, y, varcv, levels=vc)
3852    elif  contkind == 'fixc':
3853        plt.rcParams['contour.negative_linestyle'] = 'solid'
3854        coln = ckind.split(',')[1]
3855        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
3856    elif  contkind == 'fixsigc':
3857        coln = ckind.split(',')[1]
3858        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
3859    else:
3860        print errormsg
3861        print '  ' + fname + ': contour kind "' + contkind + '" not defined !!!!!'
3862        quit(-1)
3863
3864    if clabfmt is not None:
3865        plt.clabel(cplot, fmt=clabfmt)
3866        mincntS = format(vc[0], clabfmt[1:len(clabfmt)])
3867        maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)])
3868    else:
3869        mincntS = '{:g}'.format(vc[0])
3870        maxcntS = '{:g}'.format(vc[len(vc)-1])       
3871
3872    if not mapv is None:
3873        m.drawcoastlines()
3874
3875        meridians = pretty_int(nlon,xlon,5)
3876        m.drawmeridians(meridians,labels=[True,False,False,True])
3877        parallels = pretty_int(nlat,xlat,5)
3878        m.drawparallels(parallels,labels=[False,True,True,False])
3879
3880        plt.xlabel('W-E')
3881        plt.ylabel('S-N')
3882    else:
3883        plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(dimxu) + ')')
3884        plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(dimyu) + ')')
3885
3886    txpos = pretty_int(x.min(),x.max(),10)
3887    typos = pretty_int(y.min(),y.max(),10)
3888    txlabels = list(txpos)
3889    for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
3890    tylabels = list(typos)
3891    for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
3892
3893# set the limits of the plot to the limits of the data
3894    if reva0 == 'flip':
3895        if reva.split('@')[1] == 'x':
3896            plt.axis([x.max(), x.min(), y.min(), y.max()])
3897        else:
3898            plt.axis([x.min(), x.max(), y.max(), y.min()])
3899    else:
3900        plt.axis([x.min(), x.max(), y.min(), y.max()])
3901
3902    plt.xticks(txpos, txlabels)
3903    plt.yticks(typos, tylabels)
3904
3905# units labels
3906    cbar.set_label(vnames[0].replace('_','\_') + ' (' + units_lunits(uts[0]) + ')')
3907    plt.annotate(vnames[1].replace('_','\_') +' (' + units_lunits(uts[1]) + ') [' +  \
3908      mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction',    \
3909      color=coln)
3910
3911    figname = '2Dfields_shadow-contour'
3912    graphtit = vtit.replace('_','\_').replace('&','\&')
3913
3914    plt.title(graphtit)
3915   
3916    output_kind(kfig, figname, True)
3917
3918    return
3919
3920#Nvals=50
3921#vals1 = np.zeros((Nvals,Nvals), dtype= np.float)
3922#vals2 = np.zeros((Nvals,Nvals), dtype= np.float)
3923#for j in range(Nvals):
3924#    for i in range(Nvals):
3925#      vals1[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.)
3926#      vals2[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.) - Nvals/2
3927
3928#prettylev=pretty_int(-Nvals/2,Nvals/2,10)
3929
3930#plot_2D_shadow_contour(vals1, vals2, ['var1', 'var2'], np.arange(50)*1.,             \
3931#  np.arange(50)*1., ['x-axis','y-axis'], 'rainbow', 'fixc,b', "%.2f", [0, Nvals],    \
3932#  prettylev, ['$ms^{-1}$','$kJm^{-1}s^{-1}$'], 'test var1 & var2', 'pdf', False)
3933
3934def plot_2D_shadow_contour_time(varsv,varcv,vnames,valv,timv,timpos,timlab,valu,     \
3935  timeu,axist,dimn,colorbar,ckind,clabfmt,vs,vc,uts,vtit,kfig,reva,mapv):
3936    """ Adding labels and other staff to the graph
3937      varsv= 2D values to plot with shading
3938      varcv= 2D values to plot with contours
3939      vnames= variable names for the figure
3940      valv = values at the axes which is not time
3941      timv = values for the axis time
3942      timpos = positions at the axis time
3943      timlab = labes at the axis time
3944      valu = units at the axes which is not time
3945      timeu = units at the axes which is not time
3946      axist = which is the axis time
3947      dimn= dimension names to plot
3948      colorbar= name of the color bar to use
3949      ckind= contour kind
3950        'cmap': as it gets from colorbar
3951        'fixc,[colname]': fixed color [colname], all stright lines
3952        'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
3953      clabfmt= format of the labels in the contour plot (None, no labels)
3954      vs= minmum and maximum values to plot in shadow
3955      vc= vector with the levels for the contour
3956      uts= units of the variable [u-shadow, u-contour]
3957      vtit= title of the variable
3958      kfig= kind of figure (jpg, pdf, png)
3959      reva=
3960        * 'transpose': reverse the axes (x-->y, y-->x)
3961        * 'flip'@[x/y]: flip the axis x or y
3962      mapv= map characteristics: [proj],[res]
3963        see full documentation: http://matplotlib.org/basemap/
3964        [proj]: projection
3965          * 'cyl', cilindric
3966          * 'lcc', lamvbert conformal
3967        [res]: resolution:
3968          * 'c', crude
3969          * 'l', low
3970          * 'i', intermediate
3971          * 'h', high
3972          * 'f', full
3973    """
3974##    import matplotlib as mpl
3975##    mpl.use('Agg')
3976##    import matplotlib.pyplot as plt
3977    fname = 'plot_2D_shadow_contour'
3978
3979    if varsv == 'h':
3980        print fname + '_____________________________________________________________'
3981        print plot_2D_shadow_contour.__doc__
3982        quit()
3983
3984    if axist == 'x':
3985        dimxv = timv.copy()
3986        dimyv = valv.copy()
3987    else:
3988        dimxv = valv.copy()
3989        dimyv = timv.copy()
3990
3991    if reva[0:4] == 'flip':
3992        reva0 = 'flip'
3993        if len(reva.split('@')) != 2:
3994             print errormsg
3995             print '  ' + fname + ': flip is given', reva, 'but not axis!'
3996             quit(-1)
3997    else:
3998        reva0 = reva
3999
4000    if reva0 == 'transpose':
4001        if axist == 'x': 
4002            axist = 'y'
4003        else:
4004            axist = 'x'
4005
4006    if not mapv is None:
4007        if len(dimxv[:].shape) == 3:
4008            lon0 = dimxv[0,]
4009            lat0 = dimyv[0,]
4010        elif len(dimxv[:].shape) == 2:
4011            lon0 = dimxv[:]
4012            lat0 = dimyv[:]
4013        elif len(dimxv[:].shape) == 1:
4014            lon00 = dimxv[:]
4015            lat00 = dimyv[:]
4016            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
4017            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
4018
4019            for iy in range(len(lat00)):
4020                lon0[iy,:] = lon00
4021            for ix in range(len(lon00)):
4022                lat0[:,ix] = lat00
4023        if reva0 == 'transpose':
4024            print '  reversing the axes of the figure (x-->y, y-->x)!!'
4025            varsv = np.transpose(varsv)
4026            varcv = np.transpose(varcv)
4027            lon0 = np.transpose(lon0)
4028            lat0 = np.transpose(lat0)
4029
4030        map_proj=mapv.split(',')[0]
4031        map_res=mapv.split(',')[1]
4032
4033        dx = lon0.shape[1]
4034        dy = lon0.shape[0]
4035
4036        nlon = lon0[0,0]
4037        xlon = lon0[dy-1,dx-1]
4038        nlat = lat0[0,0]
4039        xlat = lat0[dy-1,dx-1]
4040
4041# Thats too much! :)
4042#        if lonlatLims is not None:
4043#            print '  ' + fname + ': cutting the domain to plot !!!!'
4044#            plt.xlim(lonlatLims[0], lonlatLims[2])
4045#            plt.ylim(lonlatLims[1], lonlatLims[3])
4046#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
4047#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
4048
4049#            if map_proj == 'cyl':
4050#                nlon = lonlatLims[0]
4051#                nlat = lonlatLims[1]
4052#                xlon = lonlatLims[2]
4053#                xlat = lonlatLims[3]
4054#            elif map_proj == 'lcc':
4055#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
4056#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
4057#                nlon =  lonlatLims[0]
4058#                xlon =  lonlatLims[2]
4059#                nlat =  lonlatLims[1]
4060#                xlat =  lonlatLims[3]
4061
4062        lon2 = lon0[dy/2,dx/2]
4063        lat2 = lat0[dy/2,dx/2]
4064
4065        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
4066          xlon, ',', xlat
4067
4068        if map_proj == 'cyl':
4069            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
4070              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
4071        elif map_proj == 'lcc':
4072            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
4073              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
4074
4075        if len(dimxv.shape) == 1:
4076            lons, lats = np.meshgrid(dimxv, dimyv)
4077        else:
4078            if len(dimxv.shape) == 3:
4079                lons = dimxv[0,:,:]
4080                lats = dimyv[0,:,:]
4081            else:
4082                lons = dimxv[:]
4083                lats = dimyv[:]
4084 
4085        x,y = m(lons,lats)
4086
4087    else:
4088        if reva0  == 'transpose':
4089            print '  reversing the axes of the figure (x-->y, y-->x)!!'
4090            varsv = np.transpose(varsv)
4091            varcv = np.transpose(varcv)
4092            dimn0 = []
4093            dimn0.append(dimn[1] + '')
4094            dimn0.append(dimn[0] + '')
4095            dimn = dimn0
4096            if len(dimyv.shape) == 2:
4097                x = np.transpose(dimyv)
4098            else:
4099                if len(dimxv.shape) == 2:
4100                    ddx = len(dimyv)
4101                    ddy = dimxv.shape[1]
4102                else:
4103                    ddx = len(dimyv)
4104                    ddy = len(dimxv)
4105   
4106                x = np.zeros((ddy,ddx), dtype=np.float)
4107                for j in range(ddy):
4108                    x[j,:] = dimyv
4109
4110            if len(dimxv.shape) == 2:
4111                y = np.transpose(dimxv)
4112            else:
4113                if len(dimyv.shape) == 2:
4114                    ddx = dimyv.shape[0]
4115                    ddy = len(dimxv)
4116                else:
4117                    ddx = len(dimyv)
4118                    ddy = len(dimxv)
4119
4120                y = np.zeros((ddy,ddx), dtype=np.float)
4121                for i in range(ddx):
4122                    y[:,i] = dimxv
4123        else:
4124            if len(dimxv.shape) == 2:
4125                x = dimxv
4126            else:
4127                if len(dimyv.shape) == 1:
4128                    x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
4129                    for j in range(len(dimyv)):
4130                        x[j,:] = dimxv
4131                else:
4132                    x = np.zeros((dimyv.shape), dtype=np.float)
4133                    if x.shape[0] == dimxv.shape[0]:
4134                        for j in range(x.shape[1]):
4135                            x[:,j] = dimxv
4136                    else:
4137                        for j in range(x.shape[0]):
4138                            x[j,:] = dimxv
4139
4140            if len(dimyv.shape) == 2:
4141                y = dimyv
4142            else:
4143                if len(dimxv.shape) == 1:
4144                    y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
4145                    for i in range(len(dimxv)):
4146                        y[:,i] = dimyv
4147                else:
4148                    y = np.zeros((dimxv.shape), dtype=np.float)
4149                    if y.shape[0] == dimyv.shape[0]:
4150                        for i in range(y.shape[1]):
4151                            y[:,i] = dimyv
4152                    else:
4153                        for i in range(y.shape[0]):
4154                            y[i,:] = dimyv
4155
4156    dx=varsv.shape[1]
4157    dy=varsv.shape[0]
4158   
4159    plt.rc('text', usetex=True)
4160
4161    if axist == 'x':
4162        valpos = pretty_int(y.min(),y.max(),10)
4163        vallabels = list(valpos)
4164        for i in range(len(vallabels)): vallabels[i] = str(vallabels[i])
4165    else:
4166        valpos = pretty_int(x.min(),x.max(),10)
4167        vallabels = list(valpos)
4168        for i in range(len(vallabels)): vallabels[i] = str(vallabels[i])
4169
4170    if reva0 == 'flip':
4171        if reva.split('@')[1] == 'x':
4172            varsv[:,0:dx-1] = varsv[:,dx-1:0:-1]
4173            varcv[:,0:dx-1] = varcv[:,dx-1:0:-1]
4174            plt.xticks(valpos, vallabels[::-1])
4175        else:
4176            varsv[0:dy-1,:] = varsv[dy-1:0:-1,:]
4177            varcv[0:dy-1,:] = varcv[dy-1:0:-1,:]
4178            plt.yticks(valpos, vallabels[::-1])
4179    else:
4180        plt.xlim(0,dx-1)
4181        plt.ylim(0,dy-1)
4182
4183    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
4184    cbar = plt.colorbar()
4185   
4186# contour
4187##
4188    contkind = ckind.split(',')[0]
4189    if contkind == 'cmap':
4190        cplot = plt.contour(x, y, varcv, levels=vc)
4191    elif  contkind == 'fixc':
4192        plt.rcParams['contour.negative_linestyle'] = 'solid'
4193        coln = ckind.split(',')[1]
4194        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
4195    elif  contkind == 'fixsigc':
4196        coln = ckind.split(',')[1]
4197        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
4198    else:
4199        print errormsg
4200        print '  ' + fname + ': contour kind "' + contkind + '" not defined !!!!!'
4201        quit(-1)
4202
4203    if clabfmt is not None:
4204        plt.clabel(cplot, fmt=clabfmt)
4205        mincntS = format(vc[0], clabfmt[1:len(clabfmt)])
4206        maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)])
4207    else:
4208        mincntS = '{:g}'.format(vc[0])
4209        maxcntS = '{:g}'.format(vc[len(vc)-1])       
4210
4211    if not mapv is None:
4212        m.drawcoastlines()
4213
4214        meridians = pretty_int(nlon,xlon,5)
4215        m.drawmeridians(meridians,labels=[True,False,False,True])
4216        parallels = pretty_int(nlat,xlat,5)
4217        m.drawparallels(parallels,labels=[False,True,True,False])
4218
4219        plt.xlabel('W-E')
4220        plt.ylabel('S-N')
4221    else:
4222        if axist == 'x':
4223            plt.xlabel(timeu)
4224            plt.xticks(timpos, timlab)
4225            plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(valu) + ')')
4226            plt.yticks(valpos, vallabels)
4227        else:
4228            plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(valu) + ')')
4229            plt.xticks(valpos, vallabels)
4230            plt.ylabel(timeu)
4231            plt.yticks(timpos, timlab)
4232
4233# set the limits of the plot to the limits of the data
4234    plt.axis([x.min(), x.max(), y.min(), y.max()])
4235
4236# units labels
4237    cbar.set_label(vnames[0].replace('_','\_') + ' (' + units_lunits(uts[0]) + ')')
4238    plt.annotate(vnames[1].replace('_','\_') +' (' + units_lunits(uts[1]) + ') [' +  \
4239      mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction',    \
4240      color=coln)
4241
4242    figname = '2Dfields_shadow-contour'
4243    graphtit = vtit.replace('_','\_').replace('&','\&')
4244
4245    plt.title(graphtit)
4246   
4247    output_kind(kfig, figname, True)
4248
4249    return
4250
4251def dxdy_lonlat(dxv,dyv,ddx,ddy):
4252    """ Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values
4253    dxdy_lonlat(dxv,dyv,Lv,lv)
4254      dx: values for the x
4255      dy: values for the y
4256      ddx: ',' list of which dimensions to use from values along x
4257      ddy: ',' list of which dimensions to use from values along y
4258    """
4259
4260    fname = 'dxdy_lonlat'
4261
4262    if ddx.find(',') > -1:
4263        dxk = 2
4264        ddxv = ddx.split(',')
4265        ddxy = int(ddxv[0])
4266        ddxx = int(ddxv[1])
4267    else:
4268        dxk = 1
4269        ddxy = int(ddx)
4270        ddxx = int(ddx)
4271
4272    if ddy.find(',') > -1:
4273        dyk = 2
4274        ddyv = ddy.split(',')
4275        ddyy = int(ddyv[0])
4276        ddyx = int(ddyv[1])
4277    else:
4278        dyk = 1
4279        ddyy = int(ddy)
4280        ddyx = int(ddy)
4281
4282    ddxxv = dxv.shape[ddxx]
4283    ddxyv = dxv.shape[ddxy]
4284    ddyxv = dyv.shape[ddyx]
4285    ddyyv = dyv.shape[ddyy]
4286
4287    slicex = []
4288    if len(dxv.shape) > 1:
4289        for idim in range(len(dxv.shape)):
4290            if idim == ddxx or idim == ddxy:
4291                slicex.append(slice(0,dxv.shape[idim]))
4292            else:
4293                slicex.append(0)
4294    else:
4295        slicex.append(slice(0,len(dxv)))
4296
4297    slicey = []
4298    if len(dyv.shape) > 1:
4299        for idim in range(len(dyv.shape)):
4300            if idim == ddyx or idim == ddyy:
4301                slicey.append(slice(0,dyv.shape[idim]))
4302            else:
4303                slicey.append(0)
4304    else:
4305        slicey.append(slice(0,len(dyv)))
4306
4307#    print '    ' + fname + ' Lluis shapes dxv:',dxv.shape,'dyv:',dyv.shape
4308#    print '    ' + fname + ' Lluis slicex:',slicex,'slicey:',slicey
4309
4310    if dxk == 2 and dyk == 2:
4311        if ddxxv != ddyxv:
4312            print errormsg
4313            print '  ' + fname + ': wrong dx dimensions! ddxx=',ddxxv,'ddyx=',ddyxv
4314            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
4315            quit(-1)
4316        if ddxyv != ddyyv:
4317            print errormsg
4318            print '  ' + fname + ': wrong dy dimensions! ddxy=',ddxyv,'ddyy=',ddyv
4319            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
4320            quit(-1)
4321        dx = ddxxv
4322        dy = ddxyv
4323
4324        print '  ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
4325        lonv = np.zeros((dy,dx), dtype=np.float)
4326        latv = np.zeros((dy,dx), dtype=np.float)
4327
4328
4329        lonv = dxv[tuple(slicex)]
4330        latv = dyv[tuple(slicey)]
4331
4332    elif dxk == 2 and dyk == 1:
4333        if not ddxxv == ddyxv and not ddxyv == ddyyv:
4334            print errormsg
4335            print '  ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv,    \
4336              'ddyx=',ddyxv,'ddyy=',ddyyv
4337            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
4338            quit(-1)
4339        dx = ddxvv
4340        dy = ddxyv
4341
4342        print '  ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
4343        lonv = np.zeros((dy,dx), dtype=np.float)
4344        latv = np.zeros((dy,dx), dtype=np.float)
4345        lonv = dxv[tuple(slicex)]
4346
4347        if ddxxv == ddyxv: 
4348            for iy in range(dy):
4349                latv[iy,:] = dyv[tuple(slicey)]
4350        else:
4351            for ix in range(dx):
4352                latv[:,ix] = dyv[tuple(slicey)]
4353
4354    elif dxk == 1 and dyk == 2:
4355        if not ddxxv == ddyxv and not ddxyv == ddyyv:
4356            print errormsg
4357            print '  ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv,    \
4358              'ddyx=',ddyxv,'ddyy=',ddyyv
4359            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
4360            quit(-1)
4361        dx = ddyxv
4362        dy = ddyyv
4363 
4364        print '  ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
4365        lonv = np.zeros((dy,dx), dtype=np.float)
4366        latv = np.zeros((dy,dx), dtype=np.float)
4367
4368        latv = dyv[tuple(slicey)]
4369
4370        if ddyxv == ddxxv: 
4371            for iy in range(dy):
4372                lonv[iy,:] = dxv[tuple(slicex)]
4373        else:
4374            for ix in range(dx):
4375                lonv[:,ix] = dxv[tuple(slicex)]
4376
4377
4378    elif dxk == 1 and dyk == 1:
4379        dx = ddxxv
4380        dy = ddyyv
4381 
4382#        print 'dx:',dx,'dy:',dy
4383
4384        lonv = np.zeros((dy,dx), dtype=np.float)
4385        latv = np.zeros((dy,dx), dtype=np.float)
4386
4387        for iy in range(dy):
4388            lonv[iy,:] = dxv[tuple(slicex)]
4389        for ix in range(dx):
4390            latv[:,ix] = dyv[tuple(slicey)]
4391
4392    return lonv,latv
4393
4394def plot_2D_shadow_line(varsv,varlv,vnames,vnamel,dimxv,dimyv,dimxu,dimyu,dimn,             \
4395  colorbar,colln,vs,uts,utl,vtit,kfig,reva,mapv,ifclose):
4396    """ Plotting a 2D field with shadows and another one with a line
4397      varsv= 2D values to plot with shading
4398      varlv= 1D values to plot with line
4399      vnames= variable names for the shadow variable in the figure
4400      vnamel= variable names for the line varibale in the figure
4401      dim[x/y]v = values at the axes of x and y
4402      dim[x/y]u = units at the axes of x and y
4403      dimn= dimension names to plot
4404      colorbar= name of the color bar to use
4405      colln= color for the line
4406      vs= minmum and maximum values to plot in shadow
4407      uts= units of the variable to shadow
4408      utl= units of the variable to line
4409      vtit= title of the variable
4410      kfig= kind of figure (jpg, pdf, png)
4411      reva=
4412        * 'transpose': reverse the axes (x-->y, y-->x)
4413        * 'flip'@[x/y]: flip the axis x or y
4414      mapv= map characteristics: [proj],[res]
4415        see full documentation: http://matplotlib.org/basemap/
4416        [proj]: projection
4417          * 'cyl', cilindric
4418          * 'lcc', lambert conformal
4419        [res]: resolution:
4420          * 'c', crude
4421          * 'l', low
4422          * 'i', intermediate
4423          * 'h', high
4424          * 'f', full
4425      ifclose= boolean value whether figure should be close (finish) or not
4426    """
4427##    import matplotlib as mpl
4428##    mpl.use('Agg')
4429##    import matplotlib.pyplot as plt
4430    fname = 'plot_2D_shadow_line'
4431
4432    if varsv == 'h':
4433        print fname + '_____________________________________________________________'
4434        print plot_2D_shadow_line.__doc__
4435        quit()
4436
4437    if reva[0:4] == 'flip':
4438        reva0 = 'flip'
4439        if len(reva.split('@')) != 2:
4440             print errormsg
4441             print '  ' + fname + ': flip is given', reva, 'but not axis!'
4442             quit(-1)
4443    else:
4444        reva0 = reva
4445
4446    if reva0 == 'transpose':
4447        print '  reversing the axes of the figure (x-->y, y-->x)!!'
4448        varsv = np.transpose(varsv)
4449        dxv = dimyv
4450        dyv = dimxv
4451        dimxv = dxv
4452        dimyv = dyv
4453
4454    if len(dimxv[:].shape) == 3:
4455        lon0 = dimxv[0,]
4456    elif len(dimxv[:].shape) == 2:
4457        lon0 = dimxv[:]
4458
4459    if len(dimyv[:].shape) == 3:
4460        lat0 = dimyv[0,]
4461    elif len(dimyv[:].shape) == 2:
4462        lat0 = dimyv[:]
4463
4464    if len(dimxv[:].shape) == 1 and len(dimyv[:].shape) == 1:
4465        lon00 = dimxv[:]
4466        lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
4467
4468        for iy in range(len(lat00)):
4469            lon0[iy,:] = lon00
4470        for ix in range(len(lon00)):
4471            lat0[:,ix] = lat00
4472
4473    if not mapv is None:
4474        map_proj=mapv.split(',')[0]
4475        map_res=mapv.split(',')[1]
4476
4477        dx = lon0.shape[1]
4478        dy = lat0.shape[0]
4479
4480        nlon = lon0[0,0]
4481        xlon = lon0[dy-1,dx-1]
4482        nlat = lat0[0,0]
4483        xlat = lat0[dy-1,dx-1]
4484
4485# Thats too much! :)
4486#        if lonlatLims is not None:
4487#            print '  ' + fname + ': cutting the domain to plot !!!!'
4488#            plt.xlim(lonlatLims[0], lonlatLims[2])
4489#            plt.ylim(lonlatLims[1], lonlatLims[3])
4490#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
4491#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
4492
4493#            if map_proj == 'cyl':
4494#                nlon = lonlatLims[0]
4495#                nlat = lonlatLims[1]
4496#                xlon = lonlatLims[2]
4497#                xlat = lonlatLims[3]
4498#            elif map_proj == 'lcc':
4499#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
4500#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
4501#                nlon =  lonlatLims[0]
4502#                xlon =  lonlatLims[2]
4503#                nlat =  lonlatLims[1]
4504#                xlat =  lonlatLims[3]
4505
4506        lon2 = lon0[dy/2,dx/2]
4507        lat2 = lat0[dy/2,dx/2]
4508
4509        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
4510          xlon, ',', xlat
4511
4512        if map_proj == 'cyl':
4513            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
4514              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
4515        elif map_proj == 'lcc':
4516            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
4517              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
4518        else:
4519            print errormsg
4520            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
4521            print '    available: cyl, lcc'
4522            quit(-1)
4523
4524        if len(dimxv.shape) == 1:
4525            lons, lats = np.meshgrid(dimxv, dimyv)
4526        else:
4527            if len(dimxv.shape) == 3:
4528                lons = dimxv[0,:,:]
4529            else:
4530                lons = dimxv[:]
4531
4532            if len(dimyv.shape) == 3:
4533                lats = dimyv[0,:,:]
4534            else:
4535                lats = dimyv[:]
4536 
4537        x,y = m(lons,lats)
4538
4539    else:
4540        if len(dimxv.shape) == 3:
4541            x = dimxv[0,:,:]
4542        elif len(dimxv.shape) == 2:
4543            x = dimxv
4544        else:
4545# Attempt of simplier way...
4546#            x = np.zeros((lon0.shape), dtype=np.float)
4547#            for j in range(lon0.shape[0]):
4548#                x[j,:] = dimxv
4549
4550## This way is too complicated and maybe not necessary ? (assuming dimxv.shape == dimyv.shape)
4551            if len(dimyv.shape) == 1:
4552                x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
4553                for j in range(len(dimxv)):
4554                    x[j,:] = dimxv
4555            else:
4556                x = np.zeros((dimyv.shape), dtype=np.float)
4557                if x.shape[0] == dimxv.shape[0]:
4558                    for j in range(x.shape[1]):
4559                        x[:,j] = dimxv
4560                else:
4561                    for j in range(x.shape[0]):
4562                        x[j,:] = dimxv
4563
4564        if len(dimyv.shape) == 3:
4565            y = dimyv[0,:,:]
4566        elif len(dimyv.shape) == 2:
4567            y = dimyv
4568        else:
4569#            y = np.zeros((lat0.shape), dtype=np.float)
4570#            for i in range(lat0.shape[1]):
4571#                x[:,i] = dimyv
4572
4573# Idem
4574            if len(dimxv.shape) == 1:
4575                y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
4576                for i in range(len(dimxv)):
4577                    y[:,i] = dimyv
4578            else:
4579                y = np.zeros((dimxv.shape), dtype=np.float)
4580                if y.shape[0] == dimyv.shape[0]:
4581                    for i in range(y.shape[1]):
4582                        y[:,i] = dimyv
4583                else:
4584                    for j in range(y.shape[0]):
4585                        y[j,:] = dimyv
4586
4587    plt.rc('text', usetex=True)
4588
4589    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
4590    cbar = plt.colorbar()
4591
4592    if not mapv is None:
4593        m.drawcoastlines()
4594
4595        meridians = pretty_int(nlon,xlon,5)
4596        m.drawmeridians(meridians,labels=[True,False,False,True])
4597        parallels = pretty_int(nlat,xlat,5)
4598        m.drawparallels(parallels,labels=[False,True,True,False])
4599
4600        plt.xlabel('W-E')
4601        plt.ylabel('S-N')
4602    else:
4603        plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(dimxu) + ')')
4604        plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(dimyu) + ')')
4605
4606# Line
4607##
4608
4609    if reva0 == 'flip' and reva.split('@')[1] == 'y':
4610        b=-np.max(y[0,:])/np.max(varlv)
4611        a=np.max(y[0,:])
4612    else:
4613        b=np.max(y[0,:])/np.max(varlv)
4614        a=0.
4615
4616    newlinv = varlv*b+a
4617    if reva0 == 'transpose':
4618        plt.plot(newlinv, x[0,:], '-', color=colln, linewidth=2)
4619    else:
4620        plt.plot(x[0,:], newlinv, '-', color=colln, linewidth=2)
4621
4622    txpos = pretty_int(x.min(),x.max(),10)
4623    typos = pretty_int(y.min(),y.max(),10)
4624    txlabels = list(txpos)
4625    for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
4626    tylabels = list(typos)
4627    for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
4628
4629    tllabels = pretty_int(np.min(varlv),np.max(varlv),len(txlabels))
4630    for it in range(len(tllabels)):
4631        yval = (tllabels[it]*b+a)
4632        plt.plot([x.max()*0.97, x.max()], [yval, yval], '-', color='k')
4633        plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)),             \
4634          xycoords='axes fraction')
4635
4636# set the limits of the plot to the limits of the data
4637    if reva0 == 'flip':
4638        if reva.split('@')[1] == 'x':
4639            plt.axis([x.max(), x.min(), y.min(), y.max()])
4640        else:
4641            plt.axis([x.min(), x.max(), y.max(), y.min()])
4642    else:
4643        plt.axis([x.min(), x.max(), y.min(), y.max()])
4644
4645    plt.tick_params(axis='y',right='off')
4646    if mapv is None:
4647        plt.xticks(txpos, txlabels)
4648        plt.yticks(typos, tylabels)
4649
4650    tllabels = pretty_int(np.min(varlv),np.max(varlv),len(txlabels))
4651    for it in range(len(tllabels)):
4652        plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)), xycoords='axes fraction')
4653
4654# units labels
4655    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
4656
4657    plt.annotate(vnamel +' (' + units_lunits(utl) + ')', xy=(0.75,0.04), 
4658      xycoords='figure fraction', color=colln)
4659    figname = '2Dfields_shadow_line'
4660    graphtit = vtit.replace('_','\_').replace('&','\&')
4661
4662    plt.title(graphtit)
4663   
4664    output_kind(kfig, figname, ifclose)
4665
4666    return
4667
4668def plot_Neighbourghood_evol(varsv, dxv, dyv, vnames, ttits, tpos, tlabels, colorbar, \
4669  Nng, vs, uts, gtit, kfig, ifclose):
4670    """ Plotting neighbourghood evolution
4671      varsv= 2D values to plot with shading
4672      vnames= shading variable name for the figure
4673      d[x/y]v= values at the axes of x and y
4674      ttits= titles of both time axis
4675      tpos= positions of the time ticks
4676      tlabels= labels of the time ticks
4677      colorbar= name of the color bar to use
4678      Nng= Number of grid points of the full side of the box (odd value)
4679      vs= minmum and maximum values to plot in shadow or:
4680        'Srange': for full range
4681        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
4682        'Saroundminmax@val': for min*val,max*val
4683        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
4684          percentile_(100-val)-median)
4685        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
4686        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
4687        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
4688           percentile_(100-val)-median)
4689      uts= units of the variable to shadow
4690      gtit= title of the graph
4691      kfig= kind of figure (jpg, pdf, png)
4692      ifclose= boolean value whether figure should be close (finish) or not
4693    """
4694    import numpy.ma as ma
4695
4696    fname = 'plot_Neighbourghood_evol'
4697
4698    if varsv == 'h':
4699        print fname + '_____________________________________________________________'
4700        print plot_Neighbourghood_evol.__doc__
4701        quit()
4702
4703    if len(varsv.shape) != 2:
4704        print errormsg
4705        print '  ' + fname + ': wrong number of dimensions of the values: ',         \
4706          varsv.shape
4707        quit(-1)
4708
4709    varsvmask = ma.masked_equal(varsv,fillValue)
4710
4711    vsend = np.zeros((2), dtype=np.float)
4712# Changing limits of the colors
4713    if type(vs[0]) != type(np.float(1.)):
4714        if vs[0] == 'Srange':
4715            vsend[0] = np.min(varsvmask)
4716        elif vs[0][0:11] == 'Saroundmean':
4717            meanv = np.mean(varsvmask)
4718            permean = np.float(vs[0].split('@')[1])
4719            minv = np.min(varsvmask)*permean
4720            maxv = np.max(varsvmask)*permean
4721            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
4722            vsend[0] = meanv-minextrm
4723            vsend[1] = meanv+minextrm
4724        elif vs[0][0:13] == 'Saroundminmax':
4725            permean = np.float(vs[0].split('@')[1])
4726            minv = np.min(varsvmask)*permean
4727            maxv = np.max(varsvmask)*permean
4728            vsend[0] = minv
4729            vsend[1] = maxv
4730        elif vs[0][0:17] == 'Saroundpercentile':
4731            medianv = np.median(varsvmask)
4732            valper = np.float(vs[0].split('@')[1])
4733            minv = np.percentile(varsvmask, valper)
4734            maxv = np.percentile(varsvmask, 100.-valper)
4735            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
4736            vsend[0] = medianv-minextrm
4737            vsend[1] = medianv+minextrm
4738        elif vs[0][0:5] == 'Smean':
4739            meanv = np.mean(varsvmask)
4740            permean = np.float(vs[0].split('@')[1])
4741            minv = np.min(varsvmask)*permean
4742            maxv = np.max(varsvmask)*permean
4743            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
4744            vsend[0] = -minextrm
4745            vsend[1] = minextrm
4746        elif vs[0][0:7] == 'Smedian':
4747            medianv = np.median(varsvmask)
4748            permedian = np.float(vs[0].split('@')[1])
4749            minv = np.min(varsvmask)*permedian
4750            maxv = np.max(varsvmask)*permedian
4751            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
4752            vsend[0] = -minextrm
4753            vsend[1] = minextrm
4754        elif vs[0][0:11] == 'Spercentile':
4755            medianv = np.median(varsvmask)
4756            valper = np.float(vs[0].split('@')[1])
4757            minv = np.percentile(varsvmask, valper)
4758            maxv = np.percentile(varsvmask, 100.-valper)
4759            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
4760            vsend[0] = -minextrm
4761            vsend[1] = minextrm
4762        else:
4763            print errormsg
4764            print '  ' + fname + ": range '" + vs[0] + "' not ready!!!"
4765            quit(-1)
4766        print '    ' + fname + ': modified shadow min,max:',vsend
4767    else:
4768        vsend[0] = vs[0]
4769
4770    if type(vs[0]) != type(np.float(1.)):
4771        if vs[1] == 'range':
4772            vsend[1] = np.max(varsv)
4773    else:
4774        vsend[1] = vs[1]
4775
4776    plt.rc('text', usetex=True)
4777
4778#    plt.pcolormesh(dxv, dyv, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
4779    plt.pcolormesh(varsvmask, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
4780    cbar = plt.colorbar()
4781
4782    newtposx = (tpos[0][:] - np.min(dxv)) * len(dxv) * Nng / (np.max(dxv) - np.min(dxv))
4783    newtposy = (tpos[1][:] - np.min(dyv)) * len(dyv) * Nng / (np.max(dyv) - np.min(dyv))
4784
4785    plt.xticks(newtposx, tlabels[0])
4786    plt.yticks(newtposy, tlabels[1])
4787    plt.xlabel(ttits[0])
4788    plt.ylabel(ttits[1])
4789
4790    plt.axes().set_aspect('equal')
4791# From: http://stackoverflow.com/questions/14406214/moving-x-axis-to-the-top-of-a-plot-in-matplotlib
4792    plt.axes().xaxis.tick_top
4793    plt.axes().xaxis.set_ticks_position('top')
4794
4795# units labels
4796    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
4797
4798    figname = 'Neighbourghood_evol'
4799    graphtit = gtit.replace('_','\_').replace('&','\&')
4800
4801    plt.title(graphtit, position=(0.5,1.05))
4802   
4803    output_kind(kfig, figname, ifclose)
4804
4805    return
4806
4807def plot_lines(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, gtit, gloc, kfig):
4808    """ Function to plot a collection of lines
4809      vardv= list of set of dimension values
4810      varvv= list of set of values
4811      vaxis= which axis will be used for the values ('x', or 'y')
4812      dtit= title for the common dimension
4813      linesn= names for the legend
4814      vtit= title for the vaxis
4815      vunit= units of the vaxis
4816      gtit= main title
4817      gloc= location of the legend (-1, autmoatic)
4818        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4819        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4820        9: 'upper center', 10: 'center'
4821      kfig= kind of figure
4822      plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)',      \
4823  ['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf')
4824    """
4825    fname = 'plot_lines'
4826
4827    if vardv == 'h':
4828        print fname + '_____________________________________________________________'
4829        print plot_lines.__doc__
4830        quit()
4831
4832# Canging line kinds every 7 lines (end of standard colors)
4833    linekinds=['.-','x-','o-']
4834
4835    Ntraj = len(vardv)
4836
4837    N7lines = 0
4838
4839    plt.rc('text', usetex=True)
4840
4841    if vaxis == 'x':
4842        for il in range(Ntraj):
4843            plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il])
4844            if il == 6: N7lines = N7lines + 1
4845
4846        plt.xlabel(vtit + ' (' + vunit + ')')
4847        plt.ylabel(dtit)
4848        plt.xlim(np.min(varvv[:]),np.max(varvv[:]))
4849        plt.ylim(np.min(vardv[:]),np.max(vardv[:]))
4850
4851    else:
4852        for il in range(Ntraj):
4853            plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il])
4854            if il == 6: N7lines = N7lines + 1
4855
4856        plt.xlabel(dtit)
4857        plt.ylabel(vtit + ' (' + vunit + ')')
4858        plt.xlim(np.min(vardv[:]),np.max(vardv[:]))
4859        plt.ylim(np.min(varvv[:]),np.max(varvv[:]))
4860
4861    figname = 'lines'
4862    graphtit = gtit.replace('_','\_').replace('&','\&')
4863
4864    plt.title(graphtit)
4865    plt.legend(loc=gloc)
4866   
4867    output_kind(kfig, figname, True)
4868
4869    return
4870
4871def plot_lines_time(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, tpos, tlabs,     \
4872  gtit, gloc, kfig):
4873    """ Function to plot a collection of lines with a time axis
4874      vardv= list of set of dimension values
4875      varvv= list of set of values
4876      vaxis= which axis will be used for the time values ('x', or 'y')
4877      dtit= title for the common dimension
4878      linesn= names for the legend
4879      vtit= title for the vaxis
4880      vunit= units of the vaxis
4881      tpos= positions of the time ticks
4882      tlabs= labels of the time ticks
4883      gtit= main title
4884      gloc= location of the legend (-1, autmoatic)
4885        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4886        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4887        9: 'upper center', 10: 'center'
4888      kfig= kind of figure
4889      plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)',      \
4890  ['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf')
4891    """
4892    fname = 'plot_lines'
4893
4894    if vardv == 'h':
4895        print fname + '_____________________________________________________________'
4896        print plot_lines.__doc__
4897        quit()
4898
4899# Canging line kinds every 7 lines (end of standard colors)
4900    linekinds=['.-','x-','o-']
4901
4902    Ntraj = len(vardv)
4903
4904    N7lines = 0
4905
4906    plt.rc('text', usetex=True)
4907    varTvv = []
4908    varTdv = []
4909
4910    if vaxis == 'x':
4911        for il in range(Ntraj):
4912            plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il])
4913            varTvv = varTvv + list(varvv[il])
4914            varTdv = varTdv + list(vardv[il])
4915            if il == 6: N7lines = N7lines + 1
4916
4917        plt.xlabel(vtit + ' (' + vunit + ')')
4918        plt.ylabel(dtit)
4919        plt.xlim(np.min(varTvv),np.max(varTvv))
4920        plt.ylim(np.min(varTdv),np.max(varTdv))
4921        plt.yticks(tpos, tlabs)
4922    else:
4923        for il in range(Ntraj):
4924            plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il])
4925            varTvv = varTvv + list(varvv[il])
4926            varTdv = varTdv + list(vardv[il])
4927            if il == 6: N7lines = N7lines + 1
4928
4929        plt.xlabel(dtit)
4930        plt.ylabel(vtit + ' (' + vunit + ')')
4931
4932        plt.xlim(np.min(varTdv),np.max(varTdv))
4933        plt.ylim(np.min(varTvv),np.max(varTvv))
4934        plt.xticks(tpos, tlabs)
4935
4936    figname = 'lines_time'
4937    graphtit = gtit.replace('_','\_').replace('&','\&')
4938
4939    plt.title(graphtit)
4940    plt.legend(loc=gloc)
4941   
4942    output_kind(kfig, figname, True)
4943
4944    return
Note: See TracBrowser for help on using the repository browser.