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

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

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