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

Last change on this file since 658 was 655, checked in by lfita, 10 years ago

Adding legend on `draw_topo_geogrid_boxes'
Fixing typo on legend position in all figures

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