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

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

Fixing the extra time-ticks issue on 'CFtime_plot'

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