Changeset 746 in lmdz_wrf


Ignore:
Timestamp:
May 4, 2016, 9:06:46 AM (9 years ago)
Author:
lfita
Message:

In progress version of the new `python script'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var_tools.py

    r744 r746  
    44import re
    55import numpy.ma as ma
     6import generic_tools as gen
    67
    78main = 'nc_var_tools.py'
    8 
    9 errormsg='ERROR -- error -- ERROR -- error'
    10 warnmsg='WARNING -- warning -- WARNING -- warning'
    11 
    12 fillValue = 1.e20
    13 fillValueF = 1.e20
    14 fillValueC = '-'
    15 fillValueI = -99999
    16 fillValueI16 = -99999
    17 fillValueF64 = 1.e20
    18 fillValueI32 = -99999
    19 
    20 def values_line(line, splitv, chars):
    21     """ Function to retrieve values from a line of ASCII text removing a series of characters
    22       line= line of ASCII text to work with
    23       splitv= character to use to split with the line
    24       chars= list of characters/strings to remove
    25     >>> values_line('0.0 0.0 1 [ 0 0 ]', ' ', ['[', ']'])
    26     ['0.0', '0.0', '1', '0', '0']
    27     """
    28     fname = 'values_line'
    29 
    30     vals = line.replace('\n', '')
    31     for ic in chars:
    32         vals = vals.replace(ic, '')
    33 
    34     vals0 = vals.split(splitv)
    35     vals = []
    36     for val in vals0:
    37         if not len(val) == 0: vals.append(val)
    38 
    39     return vals
    40 
    41 def filexist(filen, emsg, fmsg):
    42     """ Check existence of a file
    43     filen= file name
    44     emsg= general error message
    45     fmsg= message before generic text of inexistence of file
    46     """
    47     if not os.path.isfile(filen):
    48         print emsg
    49         print '  ' + fmsg + ' file "' + filen + '" does not exist !!'
    50         print emsg
    51         quit(-1)   
    529
    5310def varinfile(ncf, filen, emsg, vmsg, varn):
    … …  
    6623        quit(-1)   
    6724
     25    return
     26
    6827class testvarinfile(object):
    6928    """ Check existence of a variable inside a netCDF file
    … …  
    8847            else:
    8948                self.exist = True
     49        return
    9050
    9151def attrinvar(varobj, emsg, vmsg, attrn):
    … …  
    10565        quit(-1)   
    10666
    107 def reduce_spaces(string):
    108     """ Function to give words of a line of text removing any extra space
    109     """
    110     values = string.replace('\n','').split(' ')
    111     vals = []
    112     for val in values:
    113          if len(val) > 0:
    114              vals.append(val)
    115 
    116     return vals
    117 
    118 def reduce_tabulars(string):
    119     """ Function to give words of a line of text removing any extra space and tabular separation
    120     """
    121     values = string.replace('\n','').split('\t')
    122     vals = []
    123     for val in values:
    124          if len(val) > 0:
    125              vals.append(val)
    126 
    127     return vals
    128 
    129 def check_arguments(funcname,args,expectargs,char):
    130     """ Function to check the number of arguments if they are coincident
    131     check_arguments(funcname,Nargs,args,char)
    132       funcname= name of the function/program to check
    133       args= passed arguments
    134       expectargs= expected arguments
    135       char= character used to split the arguments
    136     """
    137 
    138     fname = 'check_arguments'
    139 
    140     Nvals = len(args.split(char))
    141     Nargs = len(expectargs.split(char))
    142 
    143     if Nvals != Nargs:
    144         print errormsg
    145         print '  ' + fname + ': wrong number of arguments:',Nvals," passed to  '",   \
    146           funcname, "' which requires:",Nargs,'!!'
    147         print '     # given expected _______'
    148         Nmax = np.max([Nvals, Nargs])
    149         for ia in range(Nmax):
    150             if ia > Nvals-1:
    151                 aval = ' '
    152             else:
    153                 aval = args.split(char)[ia]
    154             if ia > Nargs-1:
    155                 aexp = ' '
    156             else:
    157                 aexp = expectargs.split(char)[ia]
    158 
    159             print '    ', ia, aval, aexp
    160         quit(-1)
    161 
    16267    return
    163 
    164 def Str_Bool(val):
    165     """ Function to transform from a String value to a boolean one
    166     >>> Str_Bool('True')
    167     True
    168     >>> Str_Bool('0')
    169     False
    170     >>> Str_Bool('no')
    171     False
    172     """
    173 
    174     fname = 'Str_Bool'
    175 
    176     if val == 'True' or val == 'true' or val == '1' or val == 'yes' or val == 'T' or val == 't':
    177         boolv = True
    178         boolv = True
    179     elif val == 'False' or val == 'false' or val == '0' or val== 'no' or val == 'F' or val == 'f':
    180         boolv = False
    181     else:
    182         print errormsg
    183         print '  ' + fname + ": value '" + val + "' not ready!!"
    184         quit(-1)
    185 
    186     return boolv
    187 
    188 def Nchar(N,string):
    189     """ Function to provide 'N' times concatenated 'string'
    190       N= number of copies
    191       strin= string to reproduce
    192     >>> Nchar(4,'#')
    193     ####
    194     """
    195     fname = 'Nchar'
    196 
    197     newstring=''
    198     for it in range(N):
    199         newstring=newstring + string
    200 
    201     return newstring
    202 
    203 def string_dicvals_char(dictionary, string, char):
    204     """ Function to provide a string taking values from a given [string] where each single character
    205       coincide with the key of the dictionary [dictionary], separated by a given set of characters [char]
    206       dictionary= dictionary of the values dim([c]) = {[value]}
    207       string = [c1][...[cN]] String with the combination of [c]s
    208       >>> dictvals = {'i': 'dimx', 'j': 'dimy', 'k': 'dimz', 't': 'dimt'}
    209       >>> string_dicvals_char(dictvals, 'ijkt', ', ')
    210       dimx, dimy, dimz, dimt
    211     """
    212     fname = 'string_dicvals_char'
    213 
    214     Lstring = len(string)
    215     newstring = ''
    216     for ic in range(Lstring):
    217         cc = string[ic:ic+1]
    218         if not dictionary.has_key(cc):
    219             print errormsg
    220             print '  ' + fname + ": dictinoary of values does not have character '" +\
    221               cc + "' !!"
    222             print '    it has:',dictionary
    223 
    224         if ic == 0:
    225             newstring = dictionary[cc]
    226         else:
    227             newstring = newstring + char + dictionary[cc]
    228 
    229     return newstring
    230 
    231 def substring(stringv,pos,char):
    232     """ Function to substitute a character of a given string
    233       stringv= string to change
    234       pos= position to change
    235       char= characters to use as substitution
    236     >>> substring('0123456',2,'ii')
    237     01ii3456
    238     """
    239     fname = 'substring'
    240 
    241 # From: http://stackoverflow.com/questions/1228299/change-one-character-in-a-string-in-python
    242     newstring = stringv[:pos] + char + stringv[pos+1:]
    243 
    244     return newstring
    245 
    246 def period_information(idate, edate, totunits):
    247     """ Function to provide the information of a given period idate, edate (dates in [YYYY][MM][DD][HH][MI][SS] format)
    248       [idate]= initial date of the period (in [YYYY][MM][DD][HH][MI][SS] format)
    249       [edate]= end date of the period (in [YYYY][MM][DD][HH][MI][SS] format)
    250       [totunits]= latest number of entire temporal units to compute: 'year', 'month', 'day', 'hour', 'minute', 'second'
    251 
    252     resultant information given as ',' list of:
    253       [dT]: sign of the temporal increment ('pos', 'neg')
    254       [idate]: initial date as [YYYY]:[MM]:[DD]:[HH]:[MI]:[SS]
    255       [edate]: end date as [YYYY]:[MM]:[DD]:[HH]:[MI]:[SS]
    256       [Nsecs]: total seconds between dates
    257       [Nyears]: Number of years taken as entire units ([iYYYY]0101000000 to [eYYYY+1]0101000000)
    258       [ExactYears]: Exact year-dates of the period as a '@' separated list of  [YYYY]0101000000
    259       [Nmonths]: Number of months taken as entire units ([iYYYY][iMM]01000000 to [eYYYY][eMM+1]01000000)
    260       [ExactMonths]: Exact mon-dates of the period as  a '@' separated list of [YYYY][MM]01000000
    261       [Ndays]: Number of days taken as entire units ([iYYYY][iMM][iDD]000000 to [eYYYY][eMM][eDD+1]000000)
    262       [ExactDays]: Exact day-dates of the period as  a '@' separated list of [YYYY][MM][DD]000000
    263       [Nhours]: Number of hours taken as entire units ([iYYYY][iMM][iDD][iHH]0000 to [eYYYY][eMM][eDD][eHH+1]0000)
    264       [ExactHours]: Exact hour-dates of the period as a '@' separated list of [YYYY][MM][DD][HH]0000
    265       [Nminutes]: Number of minutes taken as entire units ([iYYYY][iMM][iDD][iHH][iMI]00 to [eYYYY][eMM][eDD][eHH][eMI+1]00)
    266       [ExactMinutes]: Exact minutes-dates of the period as a '@' separated list of [YYYY][MM][DD][HH][MI]00
    267       [Nsecs]: Number of minutes taken as entire units ([iYYYY][iMM][iDD][iHH][iMI][iS] to [eYYYY][eMM][eDD][eHH][eMI][eSE+1])
    268       [ExactSeconds]: Exact seconds-dates of the period as a '@' separated list of [YYYY][MM][DD][HH][MI][SS]
    269 
    270     >>> period_information( '19760217083110', '19770828000000', 'month')
    271     pos,1976:02:17:08:31:10,1977:08:28:00:00:00,48180530.0,19,19760201000000@19760301000000@19760401000000@
    272     19760501000000@19760601000000@19760701000000@19760801000000@19760901000000@19761001000000@19761101000000@
    273     19761201000000@19770101000000@19770201000000@19770301000000@19770401000000@19770501000000@19770601000000@
    274     19770701000000@19770801000000@19770901000000
    275     """
    276     import datetime as dt
    277 
    278     fname = 'period_information'
    279 
    280     if idate == 'h':
    281         print fname + '_____________________________________________________________'
    282         print variables_values.__doc__
    283         quit()
    284 
    285     expectargs = '[idate],[edate],[totunits]'
    286  
    287     check_arguments(fname,str(idate)+','+str(edate)+','+str(totunits),expectargs,',')
    288 
    289 # Checking length
    290     Lidate = len(idate)
    291     Ledate = len(edate)
    292 
    293     if Lidate != Ledate:
    294         print errormsg
    295         print '   ' + fname + ': string dates of the period have different length !!'
    296         print "    idate: '" + idate + "' (L=", Lidate, ") edate: '" + edate +       \
    297           "' (L=", Ledate, ')'
    298         quit(-1)
    299     else:
    300         if Lidate != 14:
    301             print errormsg
    302             print '  ' + fname + ": wrong initial date: '" + idate + "' must " +     \
    303               'have 14 characters!!'
    304             print '    and it has:', Lidate
    305             quit(-1)
    306         if Ledate != 14:
    307             print errormsg
    308             print '  ' + fname + ": wrong end date: '" + edate + "' must " +         \
    309               'have 14 characters!!'
    310             print '    and it has:', Ledate
    311             quit(-1)
    312 
    313 # Checking for right order of dates
    314     if int(idate) > int(edate):
    315         print warnmsg
    316         print '  ' + fname + ": initial date '" + idate + "' after end date '" +     \
    317           edate + "' !!"
    318         print "  re-sorting them!"
    319         i1 = edate
    320         e1 = idate
    321         idate = i1
    322         edate = e1
    323         oinf = 'neg'
    324     else:
    325         oinf = 'pos'
    326 
    327 # Year, month, day, hour, minute, second initial date
    328     iyrS = idate[0:4]
    329     imoS = idate[4:6]
    330     idaS = idate[6:8]
    331     ihoS = idate[8:10]
    332     imiS = idate[10:12]
    333     iseS = idate[12:14]
    334 
    335     oinf = oinf + ',' + iyrS+':'+imoS+':'+idaS+':'+ihoS+':'+imiS+':'+iseS
    336 
    337 # Year, month, day, hour, minute, second end date
    338     eyrS = edate[0:4]
    339     emoS = edate[4:6]
    340     edaS = edate[6:8]
    341     ehoS = edate[8:10]
    342     emiS = edate[10:12]
    343     eseS = edate[12:14]
    344 
    345     oinf = oinf + ',' + eyrS+':'+emoS+':'+edaS+':'+ehoS+':'+emiS+':'+eseS
    346 
    347     idateT = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),int(iseS))
    348     edateT = dt.datetime(int(eyrS),int(emoS),int(edaS),int(ehoS),int(emiS),int(eseS))
    349 
    350 # Seconds between dates
    351     DT = edateT - idateT
    352     Nsecs = DT.total_seconds()
    353 
    354     oinf = oinf + ',' + str(Nsecs)
    355 
    356 # Looking for number of years tacking exact beginning of the units [iYYYY]0101000000, [eYYYY+1]0101000000
    357 ##
    358     if totunits == 'year':
    359         itime = int(iyrS)
    360         if edate[4:15] != '0101000000':
    361             etime = int(eyrS)+1
    362         else:
    363             etime = int(eyrS)
    364 
    365         itimeT = dt.datetime(itime,1,1,0,0,0)
    366 
    367         Nyears = 1
    368         ExactYears = itimeT.strftime("%Y%m%d%H%M%S")
    369         it = int(iyrS)
    370         while it + 1 <= etime:
    371             it = it + 1
    372             Nyears = Nyears + 1
    373             itT = dt.datetime(it,1,1,0,0,0)
    374             ExactYears = ExactYears + '@' + itT.strftime("%Y%m%d%H%M%S")
    375 
    376         oinf = oinf + ',' + str(Nyears) + ',' + ExactYears
    377 
    378 # Looking for number of months tacking exact beginning of the units [iYYYY][iMM]01000000, [eYYYY][eMM+1]01000000
    379 ##
    380     if totunits == 'month':
    381         itime = int(iyrS)*100 + int(imoS)
    382         if edate[6:15] != '01000000':
    383             emo1 = int(emoS)+1
    384             if emo1 > 13:
    385                 eyr1 = str(int(eyrS) + 1)
    386                 emo1 = 01
    387             else:
    388                 eyr1 = eyrS
    389            
    390         else:
    391             eyr1 = eyrS
    392             emo1 = emoS
    393    
    394         etime = int(eyr1)*100 + int(emo1)
    395         itimeT = dt.datetime(int(iyrS),int(imoS),1,0,0,0)
    396        
    397         Nmonths = 1
    398         ExactMonths = itimeT.strftime("%Y%m%d%H%M%S")
    399         it = itime
    400         while it + 1 <= etime:
    401             it = it + 1
    402             Nmonths = Nmonths + 1
    403             ityr = it/100
    404             itmo = it - ityr*100
    405             if itmo > 12:
    406                 ityr = ityr + 1
    407                 itmo = 1
    408                 it = ityr * 100 + itmo
    409    
    410             itT = dt.datetime(ityr,itmo,1,0,0,0)
    411             ExactMonths = ExactMonths + '@' + itT.strftime("%Y%m%d%H%M%S")
    412    
    413         oinf = oinf + ',' + str(Nmonths) + ',' + ExactMonths
    414 
    415 # Looking for number of days tacking exact beginning of the units [iYYYY][iMM][iDD]000000, [eYYYY][eMM][eDD+1]000000
    416 ##
    417     if totunits == 'day':
    418         itime = dt.datetime(int(iyrS),int(imoS),int(idaS),0,0,0)
    419         if edate[8:15] != '000000':
    420             etime = dt.datetime(int(eyrS), int(emoS), int(edaS)+1,0,0,0)
    421         else:
    422             etime = edateT
    423    
    424         Ndays = 1
    425         ExactDays = itime.strftime("%Y%m%d%H%M%S")
    426         it = itime
    427         while it + dt.timedelta(days=1) <= etime:
    428             it = it + dt.timedelta(days=1)
    429             Ndays = Ndays + 1
    430             ExactDays = ExactDays + '@' + it.strftime("%Y%m%d%H%M%S")
    431    
    432         oinf = oinf + ',' + str(Ndays) + ',' + ExactDays
    433 
    434 # Looking for number of hours tacking exact beginning of the units [iYYYY][iMM][iDD][iHH]0000, [eYYYY][eMM][eDD][iHH+1]0000
    435 ##
    436     if totunits == 'hour':
    437         itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),0,0)
    438         if edate[10:15] != '0000':
    439             etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS)+1,0,0)
    440         else:
    441             etime = edateT
    442    
    443         Nhours = 1
    444         ExactHours = itime.strftime("%Y%m%d%H%M%S")
    445         it = itime
    446         while it  + dt.timedelta(hours=1) <= etime:
    447             it = it + dt.timedelta(hours=1)
    448             Nhours = Nhours + 1
    449             ExactHours = ExactHours + '@' + it.strftime("%Y%m%d%H%M%S")
    450    
    451         oinf = oinf + ',' + str(Nhours) + ',' + ExactHours
    452 
    453 # Looking for number of minutes tacking exact beginning of the units [iYYYY][iMM][iDD][iHH][iMI]00, [eYYYY][eMM][eDD][iHH][iMI+1]00
    454 ##
    455     if totunits == 'minute':
    456         itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),0)
    457         if edate[12:15] != '00':
    458             etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS), int(emiS)+1,0)
    459         else:
    460             etime = edateT
    461 
    462         Nminutes = 1
    463         ExactMinutes = itime.strftime("%Y%m%d%H%M%S")
    464         it = itime
    465         while it + dt.timedelta(minutes=1) <= etime:
    466             it = it + dt.timedelta(minutes=1)
    467             Nminutes = Nminutes + 1
    468             ExactMinutes = ExactMinutes + '@' + it.strftime("%Y%m%d%H%M%S")
    469    
    470         oinf = oinf + ',' + str(Nminutes) + ',' + ExactMinutes
    471 
    472 # Looking for number of seconds tacking exact beginning of the units [iYYYY][iMM][iDD][iHH][iMI][iSE],
    473 #   [eYYYY][eMM][eDD][iHH][iMI][iSE+1]
    474 ##
    475     if totunits == 'second':
    476         itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),int(iseS))
    477         print 'Lluis ', edate[12:15], '00'
    478         if edate[12:15] != '00':
    479             etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS), int(emiS), int(eseS)+1)
    480         else:
    481             etime = edateT
    482 
    483         Nseconds = 1
    484         ExactSeconds = itime.strftime("%Y%m%d%H%M%S")
    485         it = itime
    486         while it + dt.timedelta(seconds=1) < etime:
    487             it = it + dt.timedelta(seconds=1)
    488             Nseconds = Nseconds + 1
    489             ExactSeconds = ExactSeconds + '@' + it.strftime("%Y%m%d%H%M%S")
    490    
    491         oinf = oinf + ',' + str(Nseconds) + ',' + ExactSeconds
    492 
    493     return oinf
    494 
    495 # LluisWORKING
    496 
    497 ####### ###### ##### #### ### ## #
    498 
    499 
    500 def variables_values(varName):
    501     """ Function to provide values to plot the different variables values from ASCII file
    502       'variables_values.dat'
    503     variables_values(varName)
    504       [varName]= name of the variable
    505         return: [var name], [std name], [minimum], [maximum],
    506           [long name]('|' for spaces), [units], [color palette] (following:
    507           http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
    508      [varn]: original name of the variable
    509        NOTE: It might be better doing it with an external ASII file. But then we
    510          got an extra dependency...
    511     >>> variables_values('WRFght')
    512     ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
    513     """
    514     import subprocess as sub
    515 
    516     fname='variables_values'
    517 
    518     if varName == 'h':
    519         print fname + '_____________________________________________________________'
    520         print variables_values.__doc__
    521         quit()
    522 
    523 # This does not work....
    524 #    folderins = sub.Popen(["pwd"], stdout=sub.PIPE)
    525 #    folder = list(folderins.communicate())[0].replace('\n','')
    526 # From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
    527     folder = os.path.dirname(os.path.realpath(__file__))
    528 
    529     infile = folder + '/variables_values.dat'
    530 
    531     if not os.path.isfile(infile):
    532         print errormsg
    533         print '  ' + fname + ": File '" + infile + "' does not exist !!"
    534         quit(-1)
    535 
    536 # Variable name might come with a statistical surname...
    537     stats=['min','max','mean','stdv', 'sum']
    538 
    539 # Variables with a statistical section on their name...
    540     NOstatsvars = ['zmaxth', 'zmax_th', 'lmax_th', 'lmaxth']
    541 
    542     ifst = False
    543     if not searchInlist(NOstatsvars, varName.lower()):
    544         for st in stats:
    545             if varName.find(st) > -1:
    546                 print '    '+ fname + ": varibale '" + varName + "' with a " +       \
    547                   "statistical surname: '",st,"' !!"
    548                 Lst = len(st)
    549                 LvarName = len(varName)
    550                 varn = varName[0:LvarName - Lst]
    551                 ifst = True
    552                 break
    553     if not ifst:
    554         varn = varName
    555 
    556     ncf = open(infile, 'r')
    557 
    558     for line in ncf:
    559         if line[0:1] != '#':
    560             values = line.replace('\n','').split(',')
    561             if len(values) != 8:
    562                 print errormsg
    563                 print "problem in varibale:'", values[0],                            \
    564                   'it should have 8 values and it has',len(values)
    565                 quit(-1)
    566 
    567             if varn[0:6] == 'varDIM':
    568 # Variable from a dimension (all with 'varDIM' prefix)
    569                 Lvarn = len(varn)
    570                 varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                 \
    571                   "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1',  \
    572                    'rainbow']
    573             else:
    574                 varvals = [values[1].replace(' ',''), values[2].replace(' ',''),     \
    575                   np.float(values[3]), np.float(values[4]),values[5].replace(' ',''),\
    576                   values[6].replace(' ',''), values[7].replace(' ','')]
    577             if values[0] == varn:
    578                 ncf.close()
    579                 return varvals
    580                 break
    581 
    582     print errormsg
    583     print '  ' + fname + ": variable '" + varn + "' not defined !!!"
    584     ncf.close()
    585     quit(-1)
    586 
    587     return
    588 
    589 def variables_values_old(varName):
    590     """ Function to provide values to plot the different variables
    591     variables_values(varName)
    592       [varName]= name of the variable
    593         return: [var name], [std name], [minimum], [maximum],
    594           [long name]('|' for spaces), [units], [color palette] (following:
    595           http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
    596      [varn]: original name of the variable
    597        NOTE: It might be better doing it with an external ASII file. But then we
    598          got an extra dependency...
    599     >>> variables_values('WRFght')
    600     ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
    601     """
    602     fname='variables_values'
    603 
    604     if varName == 'h':
    605         print fname + '_____________________________________________________________'
    606         print variables_values.__doc__
    607         quit()
    608 
    609 # Variable name might come with a statistical surname...
    610     stats=['min','max','mean','stdv', 'sum']
    611 
    612     ifst = False
    613     for st in stats:
    614         if varName.find(st) > -1:
    615             print '    '+ fname + ": varibale '" + varName + "' with a statistical "+\
    616               " surname: '",st,"' !!"
    617             Lst = len(st)
    618             LvarName = len(varName)
    619             varn = varName[0:LvarName - Lst]
    620             ifst = True
    621             break
    622     if not ifst:
    623         varn = varName
    624 
    625     if varn[0:6] == 'varDIM':
    626 # Variable from a dimension (all with 'varDIM' prefix)
    627         Lvarn = len(varn)
    628         varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                         \
    629           "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1', 'rainbox']
    630     elif varn == 'a_tht' or varn == 'LA_THT':
    631         varvals = ['ath', 'total_thermal_plume_cover', 0., 1.,                       \
    632         'total|column|thermal|plume|cover', '1', 'YlGnBu']
    633     elif varn == 'acprc' or varn == 'RAINC':
    634         varvals = ['acprc', 'accumulated_cmulus_precipitation', 0., 3.e4,            \
    635           'accumulated|cmulus|precipitation', 'mm', 'Blues']
    636     elif varn == 'acprnc' or varn == 'RAINNC':
    637         varvals = ['acprnc', 'accumulated_non-cmulus_precipitation', 0., 3.e4,       \
    638           'accumulated|non-cmulus|precipitation', 'mm', 'Blues']
    639     elif varn == 'bils' or varn == 'LBILS':
    640         varvals = ['bils', 'surface_total_heat_flux', -100., 100.,                   \
    641           'surface|total|heat|flux', 'Wm-2', 'seismic']
    642     elif varn == 'landcat' or varn == 'category':
    643         varvals = ['landcat', 'land_categories', 0., 22., 'land|categories', '1',    \
    644           'rainbow']
    645     elif varn == 'c' or varn == 'QCLOUD' or varn == 'oliq' or varn == 'OLIQ':
    646         varvals = ['c', 'condensed_water_mixing_ratio', 0., 3.e-4,                   \
    647           'condensed|water|mixing|ratio', 'kgkg-1', 'BuPu']
    648     elif varn == 'ci' or varn == 'iwcon' or varn == 'LIWCON':
    649         varvals = ['ci', 'cloud_iced_water_mixing_ratio', 0., 0.0003,                \
    650          'cloud|iced|water|mixing|ratio', 'kgkg-1', 'Purples']
    651     elif varn == 'cl' or varn == 'lwcon' or varn == 'LLWCON':
    652         varvals = ['cl', 'cloud_liquidwater_mixing_ratio', 0., 0.0003,               \
    653          'cloud|liquid|water|mixing|ratio', 'kgkg-1', 'Blues']
    654     elif varn == 'cld' or varn == 'CLDFRA' or varn == 'rneb' or varn == 'lrneb' or   \
    655       varn == 'LRNEB':
    656         varvals = ['cld', 'cloud_area_fraction', 0., 1., 'cloud|fraction', '1',      \
    657           'gist_gray']
    658     elif varn == 'cldc' or varn == 'rnebcon' or varn == 'lrnebcon' or                \
    659       varn == 'LRNEBCON':
    660         varvals = ['cldc', 'convective_cloud_area_fraction', 0., 1.,                 \
    661           'convective|cloud|fraction', '1', 'gist_gray']
    662     elif varn == 'cldl' or varn == 'rnebls' or varn == 'lrnebls' or varn == 'LRNEBLS':
    663         varvals = ['cldl', 'large_scale_cloud_area_fraction', 0., 1.,                \
    664           'large|scale|cloud|fraction', '1', 'gist_gray']
    665     elif varn == 'clt' or varn == 'CLT' or varn == 'cldt' or                         \
    666       varn == 'Total cloudiness':
    667         varvals = ['clt', 'cloud_area_fraction', 0., 1., 'total|cloud|cover', '1',   \
    668           'gist_gray']
    669     elif varn == 'cll' or varn == 'cldl' or varn == 'LCLDL' or                       \
    670       varn == 'Low-level cloudiness':
    671         varvals = ['cll', 'low_level_cloud_area_fraction', 0., 1.,                   \
    672           'low|level|(p|>|680|hPa)|cloud|fraction', '1', 'gist_gray']
    673     elif varn == 'clm' or varn == 'cldm' or varn == 'LCLDM' or                       \
    674       varn == 'Mid-level cloudiness':
    675         varvals = ['clm', 'mid_level_cloud_area_fraction', 0., 1.,                   \
    676           'medium|level|(440|<|p|<|680|hPa)|cloud|fraction', '1', 'gist_gray']
    677     elif varn == 'clh' or varn == 'cldh' or varn == 'LCLDH' or                       \
    678       varn == 'High-level cloudiness':
    679         varvals = ['clh', 'high_level_cloud_area_fraction', 0., 1.,                  \
    680           'high|level|(p|<|440|hPa)|cloud|fraction', '1', 'gist_gray']
    681     elif varn == 'clmf' or varn == 'fbase' or varn == 'LFBASE':
    682         varvals = ['clmf', 'cloud_base_max_flux', -0.3, 0.3, 'cloud|base|max|flux',  \
    683           'kgm-2s-1', 'seismic']
    684     elif varn == 'clp' or varn == 'pbase' or varn == 'LPBASE':
    685         varvals = ['clp', 'cloud_base_pressure', -0.3, 0.3, 'cloud|base|pressure',   \
    686           'Pa', 'Reds']
    687     elif varn == 'cpt' or varn == 'ptconv' or varn == 'LPTCONV':
    688         varvals = ['cpt', 'convective_point', 0., 1., 'convective|point', '1',       \
    689           'seismic']
    690     elif varn == 'dqajs' or varn == 'LDQAJS':
    691         varvals = ['dqajs', 'dry_adjustment_water_vapor_tendency', -0.0003, 0.0003,  \
    692         'dry|adjustment|water|vapor|tendency', 'kg/kg/s', 'seismic']
    693     elif varn == 'dqcon' or varn == 'LDQCON':
    694         varvals = ['dqcon', 'convective_water_vapor_tendency', -3e-8, 3.e-8,         \
    695         'convective|water|vapor|tendency', 'kg/kg/s', 'seismic']
    696     elif varn == 'dqdyn' or varn == 'LDQDYN':
    697         varvals = ['dqdyn', 'dynamics_water_vapor_tendency', -3.e-7, 3.e-7,          \
    698         'dynamics|water|vapor|tendency', 'kg/kg/s', 'seismic']
    699     elif varn == 'dqeva' or varn == 'LDQEVA':
    700         varvals = ['dqeva', 'evaporation_water_vapor_tendency', -3.e-6, 3.e-6,       \
    701         'evaporation|water|vapor|tendency', 'kg/kg/s', 'seismic']
    702     elif varn == 'dqlscst' or varn == 'LDQLSCST':
    703         varvals = ['dqlscst', 'stratocumulus_water_vapor_tendency', -3.e-7, 3.e-7,   \
    704         'stratocumulus|water|vapor|tendency', 'kg/kg/s', 'seismic']
    705     elif varn == 'dqlscth' or varn == 'LDQLSCTH':
    706         varvals = ['dqlscth', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7,        \
    707         'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic']
    708     elif varn == 'dqlsc' or varn == 'LDQLSC':
    709         varvals = ['dqlsc', 'condensation_water_vapor_tendency', -3.e-6, 3.e-6,      \
    710         'condensation|water|vapor|tendency', 'kg/kg/s', 'seismic']
    711     elif varn == 'dqphy' or varn == 'LDQPHY':
    712         varvals = ['dqphy', 'physics_water_vapor_tendency', -3.e-7, 3.e-7,           \
    713         'physics|water|vapor|tendency', 'kg/kg/s', 'seismic']
    714     elif varn == 'dqthe' or varn == 'LDQTHE':
    715         varvals = ['dqthe', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7,          \
    716         'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic']
    717     elif varn == 'dqvdf' or varn == 'LDQVDF':
    718         varvals = ['dqvdf', 'vertical_difussion_water_vapor_tendency', -3.e-8, 3.e-8,\
    719         'vertical|difussion|water|vapor|tendency', 'kg/kg/s', 'seismic']
    720     elif varn == 'dqwak' or varn == 'LDQWAK':
    721         varvals = ['dqwak', 'wake_water_vapor_tendency', -3.e-7, 3.e-7,              \
    722         'wake|water|vapor|tendency', 'kg/kg/s', 'seismic']
    723     elif varn == 'dta' or varn == 'tnt' or varn == 'LTNT':
    724         varvals = ['dta', 'tendency_air_temperature', -3.e-3, 3.e-3,                 \
    725         'tendency|of|air|temperature', 'K/s', 'seismic']
    726     elif varn == 'dtac' or varn == 'tntc' or varn == 'LTNTC':
    727         varvals = ['dtac', 'moist_convection_tendency_air_temperature', -3.e-3,      \
    728         3.e-3, 'moist|convection|tendency|of|air|temperature', 'K/s', 'seismic']
    729     elif varn == 'dtar' or varn == 'tntr' or varn == 'LTNTR':
    730         varvals = ['dtar', 'radiative_heating_tendency_air_temperature', -3.e-3,     \
    731           3.e-3, 'radiative|heating|tendency|of|air|temperature', 'K/s', 'seismic']
    732     elif varn == 'dtascpbl' or varn == 'tntscpbl' or varn == 'LTNTSCPBL':
    733         varvals = ['dtascpbl',                                                       \
    734           'stratiform_cloud_precipitation_BL_mixing_tendency_air_temperature',       \
    735           -3.e-6, 3.e-6,                                                             \
    736           'stratiform|cloud|precipitation|Boundary|Layer|mixing|tendency|air|'       +
    737           'temperature', 'K/s', 'seismic']
    738     elif varn == 'dtajs' or varn == 'LDTAJS':
    739         varvals = ['dtajs', 'dry_adjustment_thermal_tendency', -3.e-5, 3.e-5,        \
    740         'dry|adjustment|thermal|tendency', 'K/s', 'seismic']
    741     elif varn == 'dtcon' or varn == 'LDTCON':
    742         varvals = ['dtcon', 'convective_thermal_tendency', -3.e-5, 3.e-5,            \
    743         'convective|thermal|tendency', 'K/s', 'seismic']
    744     elif varn == 'dtdyn' or varn == 'LDTDYN':
    745         varvals = ['dtdyn', 'dynamics_thermal_tendency', -3.e-4, 3.e-4,              \
    746         'dynamics|thermal|tendency', 'K/s', 'seismic']
    747     elif varn == 'dteva' or varn == 'LDTEVA':
    748         varvals = ['dteva', 'evaporation_thermal_tendency', -3.e-3, 3.e-3,           \
    749         'evaporation|thermal|tendency', 'K/s', 'seismic']
    750     elif varn == 'dtlscst' or varn == 'LDTLSCST':
    751         varvals = ['dtlscst', 'stratocumulus_thermal_tendency', -3.e-4, 3.e-4,       \
    752         'stratocumulus|thermal|tendency', 'K/s', 'seismic']
    753     elif varn == 'dtlscth' or varn == 'LDTLSCTH':
    754         varvals = ['dtlscth', 'thermals_thermal_tendency', -3.e-4, 3.e-4,            \
    755         'thermal|plumes|thermal|tendency', 'K/s', 'seismic']
    756     elif varn == 'dtlsc' or varn == 'LDTLSC':
    757         varvals = ['dtlsc', 'condensation_thermal_tendency', -3.e-3, 3.e-3,          \
    758         'condensation|thermal|tendency', 'K/s', 'seismic']
    759     elif varn == 'dtlwr' or varn == 'LDTLWR':
    760         varvals = ['dtlwr', 'long_wave_thermal_tendency', -3.e-3, 3.e-3, \
    761         'long|wave|radiation|thermal|tendency', 'K/s', 'seismic']
    762     elif varn == 'dtphy' or varn == 'LDTPHY':
    763         varvals = ['dtphy', 'physics_thermal_tendency', -3.e-4, 3.e-4,               \
    764         'physics|thermal|tendency', 'K/s', 'seismic']
    765     elif varn == 'dtsw0' or varn == 'LDTSW0':
    766         varvals = ['dtsw0', 'cloudy_sky_short_wave_thermal_tendency', -3.e-4, 3.e-4, \
    767         'cloudy|sky|short|wave|radiation|thermal|tendency', 'K/s', 'seismic']
    768     elif varn == 'dtthe' or varn == 'LDTTHE':
    769         varvals = ['dtthe', 'thermals_thermal_tendency', -3.e-4, 3.e-4,              \
    770         'thermal|plumes|thermal|tendency', 'K/s', 'seismic']
    771     elif varn == 'dtvdf' or varn == 'LDTVDF':
    772         varvals = ['dtvdf', 'vertical_difussion_thermal_tendency', -3.e-5, 3.e-5,    \
    773         'vertical|difussion|thermal|tendency', 'K/s', 'seismic']
    774     elif varn == 'dtwak' or varn == 'LDTWAK':
    775         varvals = ['dtwak', 'wake_thermal_tendency', -3.e-4, 3.e-4,                  \
    776         'wake|thermal|tendency', 'K/s', 'seismic']
    777     elif varn == 'ducon' or varn == 'LDUCON':
    778         varvals = ['ducon', 'convective_eastward_wind_tendency', -3.e-3, 3.e-3,      \
    779         'convective|eastward|wind|tendency', 'ms-2', 'seismic']
    780     elif varn == 'dudyn' or varn == 'LDUDYN':
    781         varvals = ['dudyn', 'dynamics_eastward_wind_tendency', -3.e-3, 3.e-3,        \
    782         'dynamics|eastward|wind|tendency', 'ms-2', 'seismic']
    783     elif varn == 'duvdf' or varn == 'LDUVDF':
    784         varvals = ['duvdf', 'vertical_difussion_eastward_wind_tendency', -3.e-3,     \
    785          3.e-3, 'vertical|difussion|eastward|wind|tendency', 'ms-2', 'seismic']
    786     elif varn == 'dvcon' or varn == 'LDVCON':
    787         varvals = ['dvcon', 'convective_difussion_northward_wind_tendency', -3.e-3,  \
    788          3.e-3, 'convective|northward|wind|tendency', 'ms-2', 'seismic']
    789     elif varn == 'dvdyn' or varn == 'LDVDYN':
    790         varvals = ['dvdyn', 'dynamics_northward_wind_tendency', -3.e-3,              \
    791          3.e-3, 'dynamics|difussion|northward|wind|tendency', 'ms-2', 'seismic']
    792     elif varn == 'dvvdf' or varn == 'LDVVDF':
    793         varvals = ['dvvdf', 'vertical_difussion_northward_wind_tendency', -3.e-3,    \
    794          3.e-3, 'vertical|difussion|northward|wind|tendency', 'ms-2', 'seismic']
    795     elif varn == 'etau' or varn == 'ZNU':
    796         varvals = ['etau', 'etau', 0., 1, 'eta values on half (mass) levels', '-',   \
    797         'reds']
    798     elif varn == 'evspsbl' or varn == 'LEVAP' or varn == 'evap' or varn == 'SFCEVPde':
    799         varvals = ['evspsbl', 'water_evaporation_flux', 0., 1.5e-4,                  \
    800           'water|evaporation|flux', 'kgm-2s-1', 'Blues']
    801     elif varn == 'evspsbl' or varn == 'SFCEVPde':
    802         varvals = ['evspsblac', 'water_evaporation_flux_ac', 0., 1.5e-4,             \
    803           'accumulated|water|evaporation|flux', 'kgm-2', 'Blues']
    804     elif varn == 'g' or varn == 'QGRAUPEL':
    805         varvals = ['g', 'grauepl_mixing_ratio', 0., 0.0003, 'graupel|mixing|ratio',  \
    806           'kgkg-1', 'Purples']
    807     elif varn == 'h2o' or varn == 'LH2O':
    808         varvals = ['h2o', 'water_mass_fraction', 0., 3.e-2,                          \
    809           'mass|fraction|of|water', '1', 'Blues']
    810     elif varn == 'h' or varn == 'QHAIL':
    811         varvals = ['h', 'hail_mixing_ratio', 0., 0.0003, 'hail|mixing|ratio',        \
    812           'kgkg-1', 'Purples']
    813     elif varn == 'hfls' or varn == 'LH' or varn == 'LFLAT' or varn == 'flat':
    814         varvals = ['hfls', 'surface_upward_latent_heat_flux', -400., 400.,           \
    815           'upward|latnt|heat|flux|at|the|surface', 'Wm-2', 'seismic']
    816     elif varn == 'hfss' or varn == 'LSENS' or varn == 'sens':
    817         varvals = ['hfss', 'surface_upward_sensible_heat_flux', -150., 150.,         \
    818           'upward|sensible|heat|flux|at|the|surface', 'Wm-2', 'seismic']
    819     elif varn == 'hus' or varn == 'WRFrh' or varn == 'LMDZrh' or varn == 'rhum' or   \
    820       varn == 'LRHUM':
    821         varvals = ['hus', 'specific_humidity', 0., 1., 'specific|humidty', '1',      \
    822           'BuPu']
    823     elif varn == 'huss' or varn == 'WRFrhs' or varn == 'LMDZrhs' or varn == 'rh2m' or\
    824       varn == 'LRH2M':
    825         varvals = ['huss', 'specific_humidity', 0., 1., 'specific|humidty|at|2m',    \
    826           '1', 'BuPu']
    827     elif varn == 'i' or varn == 'QICE':
    828         varvals = ['i', 'iced_water_mixing_ratio', 0., 0.0003,                       \
    829          'iced|water|mixing|ratio', 'kgkg-1', 'Purples']
    830     elif varn == 'lat' or varn == 'XLAT' or varn == 'XLAT_M' or varn == 'latitude':
    831         varvals = ['lat', 'latitude', -90., 90., 'latitude', 'degrees North',        \
    832           'seismic']
    833     elif varn == 'lcl' or varn == 's_lcl' or varn == 'ls_lcl' or varn == 'LS_LCL':
    834         varvals = ['lcl', 'condensation_level', 0., 2500., 'level|of|condensation',  \
    835           'm', 'Greens']
    836     elif varn == 'lambdath' or varn == 'lambda_th' or varn == 'LLAMBDA_TH':
    837         varvals = ['lambdath', 'thermal_plume_vertical_velocity', -30., 30.,         \
    838           'thermal|plume|vertical|velocity', 'm/s', 'seismic']
    839     elif varn == 'lmaxth' or varn == 'LLMAXTH':
    840         varvals = ['lmaxth', 'upper_level_thermals', 0., 100., 'upper|level|thermals'\
    841           , '1', 'Greens']
    842     elif varn == 'lon' or varn == 'XLONG' or varn == 'XLONG_M':
    843         varvals = ['lon', 'longitude', -180., 180., 'longitude', 'degrees East',     \
    844           'seismic']
    845     elif varn == 'longitude':
    846         varvals = ['lon', 'longitude', 0., 360., 'longitude', 'degrees East',        \
    847           'seismic']
    848     elif varn == 'orog' or varn == 'HGT' or varn == 'HGT_M':
    849         varvals = ['orog', 'orography',  0., 3000., 'surface|altitude', 'm','terrain']
    850     elif varn == 'pfc' or varn == 'plfc' or varn == 'LPLFC':
    851         varvals = ['pfc', 'pressure_free_convection', 100., 1100.,                   \
    852           'pressure|free|convection', 'hPa', 'BuPu']
    853     elif varn == 'plcl' or varn == 'LPLCL':
    854         varvals = ['plcl', 'pressure_lifting_condensation_level', 700., 1100.,       \
    855           'pressure|lifting|condensation|level', 'hPa', 'BuPu']
    856     elif varn == 'pr' or varn == 'RAINTOT' or varn == 'precip' or                    \
    857       varn == 'LPRECIP' or varn == 'Precip Totale liq+sol':
    858         varvals = ['pr', 'precipitation_flux', 0., 1.e-4, 'precipitation|flux',      \
    859           'kgm-2s-1', 'BuPu']
    860     elif varn == 'prprof' or varn == 'vprecip' or varn == 'LVPRECIP':
    861         varvals = ['prprof', 'precipitation_profile', 0., 1.e-3,                     \
    862           'precipitation|profile', 'kg/m2/s', 'BuPu']
    863     elif varn == 'prprofci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
    864         varvals = ['prprofci', 'precipitation_profile_convective_i', 0., 1.e-3,      \
    865           'precipitation|profile|convective|i', 'kg/m2/s', 'BuPu']
    866     elif varn == 'prprofcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
    867         varvals = ['prprofcl', 'precipitation_profile_convective_l', 0., 1.e-3,      \
    868           'precipitation|profile|convective|l', 'kg/m2/s', 'BuPu']
    869     elif varn == 'prprofli' or varn == 'pr_lsc_i' or varn == 'LPR_LSC_I':
    870         varvals = ['prprofli', 'precipitation_profile_large_scale_i', 0., 1.e-3,     \
    871           'precipitation|profile|large|scale|i', 'kg/m2/s', 'BuPu']
    872     elif varn == 'prprofll' or varn == 'pr_lsc_l' or varn == 'LPR_LSC_L':
    873         varvals = ['prprofll', 'precipitation_profile_large_scale_l', 0., 1.e-3,     \
    874           'precipitation|profile|large|scale|l', 'kg/m2/s', 'BuPu']
    875     elif varn == 'pracc' or varn == 'ACRAINTOT':
    876         varvals = ['pracc', 'precipitation_amount', 0., 100.,                        \
    877           'accumulated|precipitation', 'kgm-2', 'BuPu']
    878     elif varn == 'prc' or varn == 'LPLUC' or varn == 'pluc' or varn == 'WRFprc' or   \
    879       varn == 'RAINCde':
    880         varvals = ['prc', 'convective_precipitation_flux', 0., 2.e-4,                \
    881           'convective|precipitation|flux', 'kgm-2s-1', 'Blues']
    882     elif varn == 'prci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
    883         varvals = ['prci', 'convective_ice_precipitation_flux', 0., 0.003,           \
    884           'convective|ice|precipitation|flux', 'kgm-2s-1', 'Purples']
    885     elif varn == 'prcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
    886         varvals = ['prcl', 'convective_liquid_precipitation_flux', 0., 0.003,        \
    887           'convective|liquid|precipitation|flux', 'kgm-2s-1', 'Blues']
    888     elif varn == 'pres' or varn == 'presnivs' or varn == 'pressure' or               \
    889       varn == 'lpres' or varn == 'LPRES':
    890         varvals = ['pres', 'air_pressure', 0., 103000., 'air|pressure', 'Pa',        \
    891           'Blues']
    892     elif varn == 'prls' or varn == 'WRFprls' or varn == 'LPLUL' or varn == 'plul' or \
    893        varn == 'RAINNCde':
    894         varvals = ['prls', 'large_scale_precipitation_flux', 0., 2.e-4,              \
    895           'large|scale|precipitation|flux', 'kgm-2s-1', 'Blues']
    896     elif varn == 'prsn' or varn == 'SNOW' or varn == 'snow' or varn == 'LSNOW':
    897         varvals = ['prsn', 'snowfall', 0., 1.e-4, 'snowfall|flux', 'kgm-2s-1', 'BuPu']
    898     elif varn == 'prw' or varn == 'WRFprh':
    899         varvals = ['prw', 'atmosphere_water_vapor_content', 0., 10.,                 \
    900           'water|vapor"path', 'kgm-2', 'Blues']
    901     elif varn == 'ps' or varn == 'psfc' or varn =='PSFC' or varn == 'psol' or        \
    902       varn == 'Surface Pressure':
    903         varvals=['ps', 'surface_air_pressure', 85000., 105400., 'surface|pressure',  \
    904           'hPa', 'cool']
    905     elif varn == 'psl' or varn == 'mslp' or varn =='WRFmslp':
    906         varvals=['psl', 'air_pressure_at_sea_level', 85000., 104000.,                \
    907           'mean|sea|level|pressure', 'Pa', 'Greens']
    908     elif varn == 'qth' or varn == 'q_th' or varn == 'LQ_TH':
    909         varvals = ['qth', 'thermal_plume_total_water_content', 0., 25.,              \
    910           'total|water|cotent|in|thermal|plume', 'mm', 'YlOrRd']
    911     elif varn == 'r' or varn == 'QVAPOR' or varn == 'ovap' or varn == 'LOVAP':
    912         varvals = ['r', 'water_mixing_ratio', 0., 0.03, 'water|mixing|ratio',        \
    913           'kgkg-1', 'BuPu']
    914     elif varn == 'r2' or varn == 'Q2':
    915         varvals = ['r2', 'water_mixing_ratio_at_2m', 0., 0.03, 'water|mixing|' +     \
    916           'ratio|at|2|m','kgkg-1', 'BuPu']
    917     elif varn == 'rsds' or varn == 'SWdnSFC' or varn == 'SWdn at surface' or         \
    918       varn == 'SWDOWN':
    919         varvals=['rsds', 'surface_downwelling_shortwave_flux_in_air',  0., 1200.,    \
    920           'downward|SW|surface|radiation', 'Wm-2' ,'Reds']
    921     elif varn == 'rsdsacc':
    922         varvals=['rsdsacc', 'accumulated_surface_downwelling_shortwave_flux_in_air', \
    923           0., 1200., 'accumulated|downward|SW|surface|radiation', 'Wm-2' ,'Reds']
    924     elif varn == 'rvor' or varn == 'WRFrvor':
    925         varvals = ['rvor', 'air_relative_vorticity', -2.5E-3, 2.5E-3,                \
    926           'air|relative|vorticity', 's-1', 'seismic']
    927     elif varn == 'rvors' or varn == 'WRFrvors':
    928         varvals = ['rvors', 'surface_air_relative_vorticity', -2.5E-3, 2.5E-3,       \
    929           'surface|air|relative|vorticity', 's-1', 'seismic']
    930     elif varn == 's' or varn == 'QSNOW':
    931         varvals = ['s', 'snow_mixing_ratio', 0., 0.0003, 'snow|mixing|ratio',        \
    932           'kgkg-1', 'Purples']
    933     elif varn == 'stherm' or varn == 'LS_THERM':
    934         varvals = ['stherm', 'thermals_excess', 0., 0.8, 'thermals|excess', 'K',     \
    935           'Reds']
    936     elif varn == 'ta' or varn == 'WRFt' or varn == 'temp' or varn == 'LTEMP' or      \
    937       varn == 'Air temperature':
    938         varvals = ['ta', 'air_temperature', 195., 320., 'air|temperature', 'K',      \
    939           'YlOrRd']
    940     elif varn == 'tah' or varn == 'theta' or varn == 'LTHETA':
    941         varvals = ['tah', 'potential_air_temperature', 195., 320.,                   \
    942           'potential|air|temperature', 'K', 'YlOrRd']
    943     elif varn == 'tas' or varn == 'T2' or varn == 't2m' or varn == 'T2M' or          \
    944       varn == 'Temperature 2m':
    945         varvals = ['tas', 'air_temperature', 240., 310., 'air|temperature|at|2m', '  \
    946           K', 'YlOrRd']
    947     elif varn == 'tds' or varn == 'TH2':
    948         varvals = ['tds', 'air_dew_point_temperature', 240., 310.,                   \
    949           'air|dew|point|temperature|at|2m', 'K', 'YlGnBu']
    950     elif varn == 'tke' or varn == 'TKE' or varn == 'tke' or varn == 'LTKE':
    951         varvals = ['tke', 'turbulent_kinetic_energy', 0., 0.003,                     \
    952           'turbulent|kinetic|energy', 'm2/s2', 'Reds']
    953     elif varn == 'time'or varn == 'time_counter':
    954         varvals = ['time', 'time', 0., 1000., 'time',                                \
    955           'hours|since|1949/12/01|00:00:00', 'Reds']
    956     elif varn == 'tmla' or varn == 's_pblt' or varn == 'LS_PBLT':
    957         varvals = ['tmla', 'atmosphere_top_boundary_layer_temperature', 250., 330.,  \
    958           'atmosphere|top|boundary|layer|temperature', 'K', 'Reds']
    959     elif varn == 'ua' or varn == 'vitu' or varn == 'U' or varn == 'Zonal wind' or    \
    960       varn == 'LVITU':
    961         varvals = ['ua', 'eastward_wind', -30., 30., 'eastward|wind', 'ms-1',        \
    962           'seismic']
    963     elif varn == 'uas' or varn == 'u10m' or varn == 'U10' or varn =='Vent zonal 10m':
    964         varvals = ['uas', 'eastward_wind', -30., 30., 'eastward|2m|wind',    \
    965           'ms-1', 'seismic']
    966     elif varn == 'va' or varn == 'vitv' or varn == 'V' or varn == 'Meridional wind'  \
    967       or varn == 'LVITV':
    968         varvals = ['va', 'northward_wind', -30., 30., 'northward|wind', 'ms-1',      \
    969           'seismic']
    970     elif varn == 'vas' or varn == 'v10m' or varn == 'V10' or                         \
    971       varn =='Vent meridien 10m':
    972         varvals = ['vas', 'northward_wind', -30., 30., 'northward|2m|wind', 'ms-1',  \
    973           'seismic']
    974     elif varn == 'wakedeltaq' or varn == 'wake_deltaq' or varn == 'lwake_deltaq' or  \
    975       varn == 'LWAKE_DELTAQ':
    976         varvals = ['wakedeltaq', 'wake_delta_vapor', -0.003, 0.003,                  \
    977           'wake|delta|mixing|ratio', '-', 'seismic']
    978     elif varn == 'wakedeltat' or varn == 'wake_deltat' or varn == 'lwake_deltat' or  \
    979       varn == 'LWAKE_DELTAT':
    980         varvals = ['wakedeltat', 'wake_delta_temp', -0.003, 0.003,                   \
    981           'wake|delta|temperature', '-', 'seismic']
    982     elif varn == 'wakeh' or varn == 'wake_h' or varn == 'LWAKE_H':
    983         varvals = ['wakeh', 'wake_height', 0., 1000., 'height|of|the|wakes', 'm',    \
    984           'YlOrRd']
    985     elif varn == 'wakeomg' or varn == 'wake_omg' or varn == 'lwake_omg' or           \
    986       varn == 'LWAKE_OMG':
    987         varvals = ['wakeomg', 'wake_omega', 0., 3., 'wake|omega', \
    988           '-', 'BuGn']
    989     elif varn == 'wakes' or varn == 'wake_s' or varn == 'LWAKE_S':
    990         varvals = ['wakes', 'wake_area_fraction', 0., 0.5, 'wake|spatial|fraction',  \
    991           '1', 'BuGn']
    992     elif varn == 'wa' or varn == 'W' or varn == 'Vertical wind':
    993         varvals = ['wa', 'upward_wind', -10., 10., 'upward|wind', 'ms-1',            \
    994           'seismic']
    995     elif varn == 'wap' or varn == 'vitw' or varn == 'LVITW':
    996         varvals = ['wap', 'upward_wind', -3.e-10, 3.e-10, 'upward|wind', 'mPa-1',    \
    997           'seismic']
    998     elif varn == 'wss' or varn == 'SPDUV':
    999         varvals = ['wss', 'air_velocity',  0., 30., 'surface|horizontal|wind|speed', \
    1000           'ms-1', 'Reds']
    1001 # Water budget
    1002 # Water budget de-accumulated
    1003     elif varn == 'ccond' or varn == 'CCOND' or varn == 'ACCCONDde':
    1004         varvals = ['ccond', 'cw_cond',  0., 30.,                                     \
    1005           'cloud|water|condensation', 'mm', 'Reds']
    1006     elif varn == 'wbr' or varn == 'ACQVAPORde':
    1007         varvals = ['wbr', 'wbr',  0., 30., 'Water|Budget|water|wapor', 'mm', 'Blues']
    1008     elif varn == 'diabh' or varn == 'DIABH' or varn == 'ACDIABHde':
    1009         varvals = ['diabh', 'diabh',  0., 30., 'diabatic|heating', 'K', 'Reds']
    1010     elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde':
    1011         varvals = ['wbpw', 'water_budget_pw',  0., 30., 'Water|Budget|water|content',\
    1012            'mms-1', 'Reds']
    1013     elif varn == 'wbf' or varn == 'WBACF' or varn == 'WBACFde':
    1014         varvals = ['wbf', 'water_budget_hfcqv',  0., 30.,                       \
    1015           'Water|Budget|horizontal|convergence|of|water|vapour|(+,|' +   \
    1016           'conv.;|-,|div.)', 'mms-1', 'Reds']
    1017     elif varn == 'wbfc' or varn == 'WBFC' or varn == 'WBACFCde':
    1018         varvals = ['wbfc', 'water_budget_fc',  0., 30.,                         \
    1019           'Water|Budget|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\
    1020           'div.)', 'mms-1', 'Reds']
    1021     elif varn == 'wbfp' or varn == 'WBFP' or varn == 'WBACFPde':
    1022         varvals = ['wbfp', 'water_budget_cfp',  0., 30.,                        \
    1023           'Water|Budget|horizontal|convergence|of|precipitation|(+,|' +  \
    1024           'conv.;|-,|div.)', 'mms-1', 'Reds']
    1025     elif varn == 'wbz' or varn == 'WBZ' or varn == 'WBACZde':
    1026         varvals = ['wbz', 'water_budget_z',  0., 30.,                           \
    1027           'Water|Budget|vertical|convergence|of|water|vapour|(+,|conv.' +\
    1028           ';|-,|div.)', 'mms-1', 'Reds']
    1029     elif varn == 'wbc' or varn == 'WBC' or varn == 'WBACCde':
    1030         varvals = ['wbc', 'water_budget_c',  0., 30.,                           \
    1031           'Water|Budget|Cloud|water|species','mms-1', 'Reds']
    1032     elif varn == 'wbqvd' or varn == 'WBQVD' or varn == 'WBACQVDde':
    1033         varvals = ['wbqvd', 'water_budget_qvd',  0., 30.,                       \
    1034           'Water|Budget|water|vapour|divergence', 'mms-1', 'Reds']
    1035     elif varn == 'wbqvblten' or varn == 'WBQVBLTEN' or varn == 'WBACQVBLTENde':
    1036         varvals = ['wbqvblten', 'water_budget_qv_blten',  0., 30.,              \
    1037           'Water|Budget|QV|tendency|due|to|pbl|parameterization',        \
    1038           'kg kg-1 s-1', 'Reds']
    1039     elif varn == 'wbqvcuten' or varn == 'WBQVCUTEN' or varn == 'WBACQVCUTENde':
    1040         varvals = ['wbqvcuten', 'water_budget_qv_cuten',  0., 30.,              \
    1041           'Water|Budget|QV|tendency|due|to|cu|parameterization',         \
    1042           'kg kg-1 s-1', 'Reds']
    1043     elif varn == 'wbqvshten' or varn == 'WBQVSHTEN' or varn == 'WBACQVSHTENde':
    1044         varvals = ['wbqvshten', 'water_budget_qv_shten',  0., 30.,              \
    1045           'Water|Budget|QV|tendency|due|to|shallow|cu|parameterization', \
    1046           'kg kg-1 s-1', 'Reds']
    1047     elif varn == 'wbpr' or varn == 'WBP' or varn == 'WBACPde':
    1048         varvals = ['wbpr', 'water_budget_pr',  0., 30.,                         \
    1049           'Water|Budget|recipitation', 'mms-1', 'Reds']
    1050     elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde':
    1051         varvals = ['wbpw', 'water_budget_pw',  0., 30.,                         \
    1052           'Water|Budget|water|content', 'mms-1', 'Reds']
    1053     elif varn == 'wbcondt' or varn == 'WBCONDT' or varn == 'WBACCONDTde':
    1054         varvals = ['wbcondt', 'water_budget_condt',  0., 30.,                   \
    1055           'Water|Budget|condensation|and|deposition', 'mms-1', 'Reds']
    1056     elif varn == 'wbqcm' or varn == 'WBQCM' or varn == 'WBACQCMde':
    1057         varvals = ['wbqcm', 'water_budget_qcm',  0., 30.,                       \
    1058           'Water|Budget|hydrometeor|change|and|convergence', 'mms-1', 'Reds']
    1059     elif varn == 'wbsi' or varn == 'WBSI' or varn == 'WBACSIde':
    1060         varvals = ['wbsi', 'water_budget_si',  0., 30.,                         \
    1061           'Water|Budget|hydrometeor|sink', 'mms-1', 'Reds']
    1062     elif varn == 'wbso' or varn == 'WBSO' or varn == 'WBACSOde':
    1063         varvals = ['wbso', 'water_budget_so',  0., 30.,                         \
    1064           'Water|Budget|hydrometeor|source', 'mms-1', 'Reds']
    1065 # Water Budget accumulated
    1066     elif varn == 'ccondac' or varn == 'ACCCOND':
    1067         varvals = ['ccondac', 'cw_cond_ac',  0., 30.,                                \
    1068           'accumulated|cloud|water|condensation', 'mm', 'Reds']
    1069     elif varn == 'rac' or varn == 'ACQVAPOR':
    1070         varvals = ['rac', 'ac_r',  0., 30., 'accumualted|water|wapor', 'mm', 'Blues']
    1071     elif varn == 'diabhac' or varn == 'ACDIABH':
    1072         varvals = ['diabhac', 'diabh_ac',  0., 30., 'accumualted|diabatic|heating',  \
    1073           'K', 'Reds']
    1074     elif varn == 'wbpwac' or varn == 'WBACPW':
    1075         varvals = ['wbpwac', 'water_budget_pw_ac',  0., 30.,                         \
    1076           'Water|Budget|accumulated|water|content', 'mm', 'Reds']
    1077     elif varn == 'wbfac' or varn == 'WBACF':
    1078         varvals = ['wbfac', 'water_budget_hfcqv_ac',  0., 30.,                       \
    1079           'Water|Budget|accumulated|horizontal|convergence|of|water|vapour|(+,|' +   \
    1080           'conv.;|-,|div.)', 'mm', 'Reds']
    1081     elif varn == 'wbfcac' or varn == 'WBACFC':
    1082         varvals = ['wbfcac', 'water_budget_fc_ac',  0., 30.,                         \
    1083           'Water|Budget|accumulated|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\
    1084           'div.)', 'mm', 'Reds']
    1085     elif varn == 'wbfpac' or varn == 'WBACFP':
    1086         varvals = ['wbfpac', 'water_budget_cfp_ac',  0., 30.,                        \
    1087           'Water|Budget|accumulated|horizontal|convergence|of|precipitation|(+,|' +  \
    1088           'conv.;|-,|div.)', 'mm', 'Reds']
    1089     elif varn == 'wbzac' or varn == 'WBACZ':
    1090         varvals = ['wbzac', 'water_budget_z_ac',  0., 30.,                           \
    1091           'Water|Budget|accumulated|vertical|convergence|of|water|vapour|(+,|conv.' +\
    1092           ';|-,|div.)', 'mm', 'Reds']
    1093     elif varn == 'wbcac' or varn == 'WBACC':
    1094         varvals = ['wbcac', 'water_budget_c_ac',  0., 30.,                           \
    1095           'Water|Budget|accumulated|Cloud|water|species','mm', 'Reds']
    1096     elif varn == 'wbqvdac' or varn == 'WBACQVD':
    1097         varvals = ['wbqvdac', 'water_budget_qvd_ac',  0., 30.,                       \
    1098           'Water|Budget|accumulated|water|vapour|divergence', 'mm', 'Reds']
    1099     elif varn == 'wbqvbltenac' or varn == 'WBACQVBLTEN':
    1100         varvals = ['wbqvbltenac', 'water_budget_qv_blten_ac',  0., 30.,              \
    1101           'Water|Budget|accumulated|QV|tendency|due|to|pbl|parameterization',        \
    1102           'kg kg-1 s-1', 'Reds']
    1103     elif varn == 'wbqvcutenac' or varn == 'WBACQVCUTEN':
    1104         varvals = ['wbqvcutenac', 'water_budget_qv_cuten_ac',  0., 30.,              \
    1105           'Water|Budget|accumulated|QV|tendency|due|to|cu|parameterization',         \
    1106           'kg kg-1 s-1', 'Reds']
    1107     elif varn == 'wbqvshtenac' or varn == 'WBACQVSHTEN':
    1108         varvals = ['wbqvshtenac', 'water_budget_qv_shten_ac',  0., 30.,              \
    1109           'Water|Budget|accumulated|QV|tendency|due|to|shallow|cu|parameterization', \
    1110           'kg kg-1 s-1', 'Reds']
    1111     elif varn == 'wbprac' or varn == 'WBACP':
    1112         varvals = ['wbprac', 'water_budget_pr_ac',  0., 30.,                         \
    1113           'Water|Budget|accumulated|precipitation', 'mm', 'Reds']
    1114     elif varn == 'wbpwac' or varn == 'WBACPW':
    1115         varvals = ['wbpwac', 'water_budget_pw_ac',  0., 30.,                         \
    1116           'Water|Budget|accumulated|water|content', 'mm', 'Reds']
    1117     elif varn == 'wbcondtac' or varn == 'WBACCONDT':
    1118         varvals = ['wbcondtac', 'water_budget_condt_ac',  0., 30.,                   \
    1119           'Water|Budget|accumulated|condensation|and|deposition', 'mm', 'Reds']
    1120     elif varn == 'wbqcmac' or varn == 'WBACQCM':
    1121         varvals = ['wbqcmac', 'water_budget_qcm_ac',  0., 30.,                       \
    1122           'Water|Budget|accumulated|hydrometeor|change|and|convergence', 'mm', 'Reds']
    1123     elif varn == 'wbsiac' or varn == 'WBACSI':
    1124         varvals = ['wbsiac', 'water_budget_si_ac',  0., 30.,                         \
    1125           'Water|Budget|accumulated|hydrometeor|sink', 'mm', 'Reds']
    1126     elif varn == 'wbsoac' or varn == 'WBACSO':
    1127         varvals = ['wbsoac', 'water_budget_so_ac',  0., 30.,                         \
    1128           'Water|Budget|accumulated|hydrometeor|source', 'mm', 'Reds']
    1129 
    1130     elif varn == 'xtime' or varn == 'XTIME':
    1131         varvals = ['xtime', 'time',  0., 1.e5, 'time',                               \
    1132           'minutes|since|simulation|start', 'Reds']
    1133     elif varn == 'x' or varn == 'X':
    1134         varvals = ['x', 'x',  0., 100., 'x', '-', 'Reds']
    1135     elif varn == 'y' or varn == 'Y':
    1136         varvals = ['y', 'y',  0., 100., 'y', '-', 'Blues']
    1137     elif varn == 'z' or varn == 'Z':
    1138         varvals = ['z', 'z',  0., 100., 'z', '-', 'Greens']
    1139     elif varn == 'zg' or varn == 'WRFght' or varn == 'Geopotential height' or        \
    1140       varn == 'geop' or varn == 'LGEOP':
    1141         varvals = ['zg', 'geopotential_height', 0., 80000., 'geopotential|height',   \
    1142           'm2s-2', 'rainbow']
    1143     elif varn == 'zmaxth' or varn == 'zmax_th'  or varn == 'LZMAX_TH':
    1144         varvals = ['zmaxth', 'thermal_plume_height', 0., 4000.,                     \
    1145           'maximum|thermals|plume|height', 'm', 'YlOrRd']
    1146     elif varn == 'zmla' or varn == 's_pblh' or varn == 'LS_PBLH':
    1147         varvals = ['zmla', 'atmosphere_boundary_layer_thickness', 0., 2500.,         \
    1148           'atmosphere|boundary|layer|thickness', 'm', 'Blues']
    1149     else:
    1150         print errormsg
    1151         print '  ' + fname + ": variable '" + varn + "' not defined !!!"
    1152         quit(-1)
    1153 
    1154     return varvals
    1155 
    1156 def numVector_String(vec,char):
    1157     """ Function to transform a vector of numbers to a single string [char] separated
    1158     numVector_String(vec,char)
    1159       vec= vector with the numerical values
    1160       char= single character to split the values
    1161     >>> print numVector_String(np.arange(10),' ')
    1162     0 1 2 3 4 5 6 7 8 9
    1163     """
    1164     fname = 'numVector_String'
    1165 
    1166 #    if vec == 'h':
    1167 #        print fname + '_____________________________________________________________'
    1168 #        print numVector_String.__doc__
    1169 #        quit()
    1170 
    1171     Nvals = len(vec)
    1172 
    1173     string=''
    1174     for i in range(Nvals):
    1175         if i == 0:
    1176             string = str(vec[i])
    1177         else:
    1178             string = string + char + str(vec[i])
    1179 
    1180     return string
    1181 
    1182 def index_mat(mat,val):
    1183     """ Function to provide the coordinates of a given value inside a matrix
    1184     index_mat(mat,val)
    1185       mat= matrix with values
    1186       val= value to search
    1187     >>> index_mat(np.arange(27).reshape(3,3,3),22)
    1188     [2 1 1]
    1189     """
    1190 
    1191     fname = 'index_mat'
    1192 
    1193     matshape = mat.shape
    1194 
    1195     matlist = list(mat.flatten())
    1196     ifound = matlist.index(val)
    1197 
    1198     Ndims = len(matshape)
    1199     valpos = np.zeros((Ndims), dtype=int)
    1200     baseprevdims = np.zeros((Ndims), dtype=int)
    1201 
    1202     for dimid in range(Ndims):
    1203         baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
    1204         if dimid == 0:
    1205             alreadyplaced = 0
    1206         else:
    1207             alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
    1208         valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
    1209 
    1210     return valpos
    1211 
    1212 def multi_index_mat(mat,val):
    1213     """ Function to provide the multiple coordinates of a given value inside a matrix
    1214     index_mat(mat,val)
    1215       mat= matrix with values
    1216       val= value to search
    1217     >>> vals = np.ones((24), dtype=np.float).reshape(2,3,4)
    1218     vals[:,:,2] = 0.
    1219     vals[1,:,:] = np.pi
    1220     vals[:,2,:] = -1.
    1221     multi_index_mat(vals,1.)
    1222     [array([0, 0, 0]), array([0, 0, 1]), array([0, 0, 3]), array([0, 1, 0]), array([0, 1, 1]), array([0, 1, 3])]
    1223     """
    1224     fname = 'multi_index_mat'
    1225 
    1226     matshape = mat.shape
    1227 
    1228     ivalpos = []
    1229     matlist = list(mat.flatten())
    1230     Lmatlist = len(matlist)
    1231    
    1232     val0 = val - val
    1233     if val != val0:
    1234         valdiff = val0
    1235     else:
    1236         valdiff = np.ones((1), dtype = type(val))
    1237    
    1238     ifound = 0
    1239     while ifound < Lmatlist:
    1240         if matlist.count(val) == 0:
    1241             ifound = Lmatlist + 1
    1242         else:
    1243             ifound = matlist.index(val)
    1244 
    1245             Ndims = len(matshape)
    1246             valpos = np.zeros((Ndims), dtype=int)
    1247             baseprevdims = np.zeros((Ndims), dtype=int)
    1248 
    1249             for dimid in range(Ndims):
    1250                 baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
    1251                 if dimid == 0:
    1252                     alreadyplaced = 0
    1253                 else:
    1254                     alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
    1255                 valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
    1256             matlist[ifound] = valdiff
    1257             ivalpos.append(valpos)
    1258 
    1259     return ivalpos
    1260 
    1261 def addfileInfile(origfile,destfile,addfile,addsign):
    1262     """ Function to add the content of a file [addfile] to another one [origfile] at
    1263       a given position [addsign] creating a new file called [destfile]
    1264       origfile= original file
    1265       destfile= destination file
    1266       addfile= content of the file to add
    1267       addsign= sign where to add the file
    1268     """
    1269     fname = 'addfileInfile'
    1270 
    1271     if not os.path.isfile(origfile):
    1272         print errormsg
    1273         print '  ' + fname + ": original file '" + origfile + "' does not exist !!"
    1274         quit()
    1275 
    1276     if not os.path.isfile(addfile):
    1277         print errormsg
    1278         print '  ' + fname + ": adding file '" + addfile + "' does not exist !!"
    1279         quit()
    1280    
    1281     ofile = open(origfile, 'r')
    1282 
    1283 # Inspecting flag
    1284 ##
    1285     Nfound = 0
    1286     for line in ofile:
    1287         if line == addsign + '\n': Nfound = Nfound + 1
    1288 
    1289     if Nfound == 0:
    1290         print errormsg
    1291         print '  ' + fname + ": adding sign '" + addsign + "' not found !!"
    1292         quit()
    1293     print '  '+ fname + ': found', Nfound," times sign '" + addsign + "' "
    1294 
    1295     ofile.seek(0)
    1296     dfile = open(destfile, 'w')
    1297 
    1298     for line in ofile:
    1299         if line != addsign + '\n':
    1300             dfile.write(line)
    1301         else:
    1302             afile = open(addfile,'r')
    1303             for aline in afile:
    1304                 print aline
    1305                 dfile.write(aline)
    1306 
    1307             afile.close()
    1308 
    1309     ofile.close()
    1310     dfile.close()
    1311 
    1312     print main + " successful writting of file '" + destfile + "' !!"
    1313     return
    1314 
    1315 
    1316 ####### ###### ##### #### ### ## #
    1317 
    1318 def valmodoper(varVal, valuesS):
    1319     """ Function to run the modification of a variable
    1320       varVAl= matrix with the values
    1321       valuesS= [modins],[[modval1],...,[modvalN]] modification instruction, value with which modify
    1322         [modins]: Modifiers:
    1323           sumc: add [modval1]
    1324           subc: remove [modval1]
    1325           mulc: multiply by [modval1]
    1326           divc: divide by [modval1]
    1327           lowthres: if [val] < [modval1]; val = [modval2]
    1328           upthres: if [val] > [modval1]; val = [modval2]
    1329           lowthres@oper: if [val] < [modval1]; val = [oper] (operation as [modval2],[modval3])
    1330           upthres@oper: if [val] > [modval1]; val = [oper] (operation as [modval2],[modval3])
    1331           potc: [val] ** [modval1]
    1332     """
    1333     fname='valmodoper'
    1334 
    1335     if valuesS == 'h':
    1336         print fname + '_____________________________________________________________'
    1337         print valmodoper.__doc__
    1338         quit()
    1339 
    1340     valsS = valuesS.split(',')
    1341     modins = valsS[0]
    1342     modval = float(valsS[1])
    1343 
    1344     if modins == 'sumc':
    1345         varVal[:] = varVal[:] + modval
    1346     elif modins == 'subc':
    1347         varVal[:] = varVal[:] - modval
    1348     elif modins == 'mulc':
    1349         varVal[:] = varVal[:] * modval
    1350     elif modins == 'divc':
    1351         varVal[:] = varVal[:] / modval
    1352     elif modins == 'lowthres':
    1353         varVal2 = np.where(varVal[:] < float(valsS[1]), float(valsS[2]), varVal[:])
    1354         varVal[:] = varVal2
    1355     elif modins == 'upthres':
    1356         varVal2 = np.where(varVal[:] > float(valsS[1]), float(valsS[2]), varVal[:])
    1357         varVal[:] = varVal2
    1358     elif modins == 'lowthres@oper':
    1359         if valsS[2] == 'sumc':
    1360            varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
    1361              varVal[:] + float(valsS[3]), varVal[:])
    1362            varVal[:] = varVal2
    1363         elif valsS[2] == 'subc':
    1364            varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
    1365              varVal[:] - float(valsS[3]), varVal[:])
    1366            varVal[:] = varVal2
    1367         elif valsS[2] == 'mulc':
    1368            varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
    1369              varVal[:] * float(valsS[3]), varVal[:])
    1370            varVal[:] = varVal2
    1371         elif valsS[2] == 'divc':
    1372            varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
    1373              varVal[:] / float(valsS[3]), varVal[:])
    1374            varVal[:] = varVal2
    1375         elif valsS[2] == 'potc':
    1376            varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
    1377              varVal[:] ** float(valsS[3]), varVal[:])
    1378            varVal[:] = varVal2
    1379         else:
    1380             print errormsg
    1381             print '  ' + fname + ": Operation to modify values '" + modins +         \
    1382               "' is not defined !!"
    1383             quit(-1)
    1384     elif modins == 'upthres@oper':
    1385         if valsS[2] == 'sumc':
    1386            varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
    1387              varVal[:] + float(valsS[3]), varVal[:])
    1388            varVal[:] = varVal2
    1389         elif valsS[2] == 'subc':
    1390            varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
    1391              varVal[:] - float(valsS[3]), varVal[:])
    1392            varVal[:] = varVal2
    1393         elif valsS[2] == 'mulc':
    1394            varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
    1395              varVal[:] * float(valsS[3]), varVal[:])
    1396            varVal[:] = varVal2
    1397         elif valsS[2] == 'divc':
    1398            varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
    1399              varVal[:] / float(valsS[3]), varVal[:])
    1400            varVal[:] = varVal2
    1401         elif valsS[2] == 'potc':
    1402            varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
    1403              varVal[:] ** float(valsS[3]), varVal[:])
    1404            varVal[:] = varVal2
    1405         else:
    1406             print errormsg
    1407             print '  ' + fname + ": Operation to modify values '" + modins +         \
    1408               "' is not defined !!"
    1409             quit(-1)
    1410     elif modins == 'potc':
    1411         varVal[:] = varVal[:] ** modval
    1412     else:
    1413         print errormsg
    1414         print '  ' + fname + ": Operation to modify values '" + modins +             \
    1415           "' is not defined !!"
    1416         quit(-1)
    1417 
    1418     return varVal
    141968
    142069def valmod(values, ncfile, varn):
    … …  
    1519168    ncf.close()
    1520169
    1521 def rangedim(end, shape):
    1522     """Gives the instruction to retrieve values from a dimension of a variable
    1523     >>> print rangedim(-1, 15)
    1524     15
    1525     """
    1526     if end == -1:
    1527       return shape
    1528     else:
    1529       return end
     170    return
    1530171
    1531172def varaddref(values, ncfile, varn):
    … …  
    1723364  ncf.sync()
    1724365  ncf.close()
     366
     367  return
    1725368
    1726369def varoutold(values, ncfile, varn):
    … …  
    1819462      print '%2s %f' % ( 'NC', varValrange[i] )
    1820463
     464  return
    1821465
    1822466def varout(values, ncfile, varn):
    … …  
    1957601                  ncf = NetCDFFile(ncfile,'a')
    1958602
     603    return
     604
    1959605def chvarname(values, ncfile, varn):
    1960606  """Changing the name of the variable
    … …  
    1994640  ncf.close()
    1995641
    1996 def searchInlist(listname, nameFind):
    1997     """ Function to search a value within a list
    1998     listname = list
    1999     nameFind = value to find
    2000     >>> searInlist(['1', '2', '3', '5'], '5')
    2001     True
    2002     """
    2003     for x in listname:
    2004       if x == nameFind:
    2005         return True
    2006     return False
     642  return
    2007643
    2008644def set_attribute(ncv, attrname, attrvalue):
    … …  
    2124760  ncf.close()
    2125761
     762  return
     763
    2126764def gaddattrk(values, ncfile):
    2127765  """ Add a global attribute to a netCDF caring about the type. Removes previous attribute if it exist
    … …  
    2166804  ncf.close()
    2167805
     806  return
     807
    2168808def varaddattr(values, ncfile, varn):
    2169809  """ Add an attribute to a variable. Removes previous attribute if it exists
    … …  
    2207847  ncf.sync()
    2208848  ncf.close()
     849
     850  return
    2209851
    2210852def varaddattrk(values, ncfile, varn):
    … …  
    2273915  ncf.close()
    2274916
     917  return
     918
    2275919def varrmattr(values, ncfile, varn):
    2276920  """ Removing an attribute from a variable
    … …  
    2338982  ncf.sync()
    2339983  ncf.close()
     984
     985  return
    2340986
    2341987def fvaradd(values, ncfile):
    … …  
    25051151  ncref.close()
    25061152
     1153  return
     1154
    25071155def fattradd(var, values, ncfile):
    25081156  """ Adding attributes from a reference file
    … …  
    25621210  ncref.close()
    25631211
     1212  return
     1213
    25641214def fgaddattr(values, ncfile):
    25651215  """ Adding global attributes from a reference file
    … …  
    26401290  shu.copyfile('tmp_py.nc', ncfile)
    26411291  os.remove('tmp_py.nc')
     1292
     1293  return
    26421294 
    26431295def ivars(ncfile):
    … …  
    26611313  print '  # allvars= ' + allvars
    26621314  ncf.close()
     1315
     1316  return
    26631317
    26641318def igattrs(ncfile):
    … …  
    29011555  ncf.close()
    29021556
     1557  return
     1558
    29031559def vrattr(values, ncfile, varn):
    29041560    """ Function to remove an atribute from a variable
    … …  
    29331589
    29341590    return
    2935 
    2936 def datetimeStr_datetime(StringDT):
    2937     """ Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
    2938     >>> datetimeStr_datetime('1976-02-17_00:00:00')
    2939     1976-02-17 00:00:00
    2940     """
    2941     import datetime as dt
    2942 
    2943     fname = 'datetimeStr_datetime'
    2944 
    2945     dateD = np.zeros((3), dtype=int)
    2946     timeT = np.zeros((3), dtype=int)
    2947 
    2948     dateD[0] = int(StringDT[0:4])
    2949     dateD[1] = int(StringDT[5:7])
    2950     dateD[2] = int(StringDT[8:10])
    2951 
    2952     trefT = StringDT.find(':')
    2953     if not trefT == -1:
    2954 #        print '  ' + fname + ': refdate with time!'
    2955         timeT[0] = int(StringDT[11:13])
    2956         timeT[1] = int(StringDT[14:16])
    2957         timeT[2] = int(StringDT[17:19])
    2958 
    2959     if int(dateD[0]) == 0:
    2960         print warnmsg
    2961         print '    ' + fname + ': 0 reference year!! changing to 1'
    2962         dateD[0] = 1
    2963  
    2964     newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])
    2965 
    2966     return newdatetime
    2967 
    2968 def dateStr_date(StringDate):
    2969   """ Function to transform a string date ([YYYY]-[MM]-[DD] format) to a date object
    2970   >>> dateStr_date('1976-02-17')
    2971   1976-02-17
    2972   """
    2973   import datetime as dt
    2974 
    2975   dateD = StringDate.split('-')
    2976   if int(dateD[0]) == 0:
    2977     print warnmsg
    2978     print '    dateStr_date: 0 reference year!! changing to 1'
    2979     dateD[0] = 1
    2980   newdate = dt.date(int(dateD[0]), int(dateD[1]), int(dateD[2]))
    2981   return newdate
    2982 
    2983 def timeStr_time(StringDate):
    2984   """ Function to transform a string date ([HH]:[MI]:[SS] format) to a time object
    2985   >>> datetimeStr_datetime('04:32:54')
    2986   04:32:54
    2987   """
    2988   import datetime as dt
    2989 
    2990   timeT = StringDate.split(':')
    2991   if len(timeT) == 3:
    2992     newtime = dt.time(int(timeT[0]), int(timeT[1]), int(timeT[2]))
    2993   else:
    2994     newtime = dt.time(int(timeT[0]), int(timeT[1]), 0)
    2995 
    2996   return newtime
    2997 
    2998 def timeref_datetime(refd, timeval, tu):
    2999     """ Function to transform from a [timeval] in [tu] units from the time referece [tref] to datetime object
    3000     refd: time of reference (as datetime object)
    3001     timeval: time value (as [tu] from [tref])
    3002     tu: time units
    3003     >>> timeref = date(1949,12,1,0,0,0)
    3004     >>> timeref_datetime(timeref, 229784.36, hours)
    3005     1976-02-17 08:21:36
    3006     """
    3007     import datetime as dt
    3008     import numpy as np
    3009 
    3010 ## Not in timedelta
    3011 #    if tu == 'years':
    3012 #        realdate = refdate + dt.timedelta(years=float(timeval))
    3013 #    elif tu == 'months':
    3014 #        realdate = refdate + dt.timedelta(months=float(timeval))
    3015     if tu == 'weeks':
    3016         realdate = refd + dt.timedelta(weeks=float(timeval))
    3017     elif tu == 'days':
    3018         realdate = refd + dt.timedelta(days=float(timeval))
    3019     elif tu == 'hours':
    3020         realdate = refd + dt.timedelta(hours=float(timeval))
    3021     elif tu == 'minutes':
    3022         realdate = refd + dt.timedelta(minutes=float(timeval))
    3023     elif tu == 'seconds':
    3024         realdate = refd + dt.timedelta(seconds=float(timeval))
    3025     elif tu == 'milliseconds':
    3026         realdate = refd + dt.timedelta(milliseconds=float(timeval))
    3027     else:
    3028           print errormsg
    3029           print '    timeref_datetime: time units "' + tu + '" not ready!!!!'
    3030           quit(-1)
    3031 
    3032     return realdate
    3033 
    3034 def timeref_datetime_mat(refd, timeval, tu):
    3035     """ Function to transform from a [timeval] in [tu] units from the time referece [tref] to matrix with: year, day, month, hour, minute, second
    3036     refd: time of reference (as datetime object)
    3037     timeval: time value (as [tu] from [tref])
    3038     tu: time units
    3039     >>> timeref = date(1949,12,1,0,0,0)
    3040     >>> timeref_datetime(timeref, 229784.36, hours)
    3041     [1976    2   17    8   36   21]
    3042     """
    3043     import datetime as dt
    3044     import numpy as np
    3045 
    3046 
    3047     realdates = np.zeros(6, dtype=int)
    3048 ## Not in timedelta
    3049 #    if tu == 'years':
    3050 #        realdate = refdate + dt.timedelta(years=float(timeval))
    3051 #    elif tu == 'months':
    3052 #        realdate = refdate + dt.timedelta(months=float(timeval))
    3053     if tu == 'weeks':
    3054         realdate = refd + dt.timedelta(weeks=float(timeval))
    3055     elif tu == 'days':
    3056         realdate = refd + dt.timedelta(days=float(timeval))
    3057     elif tu == 'hours':
    3058         realdate = refd + dt.timedelta(hours=float(timeval))
    3059     elif tu == 'minutes':
    3060         realdate = refd + dt.timedelta(minutes=float(timeval))
    3061     elif tunits == 'seconds':
    3062         realdate = refd + dt.timedelta(seconds=float(timeval))
    3063     elif tunits == 'milliseconds':
    3064         realdate = refd + dt.timedelta(milliseconds=float(timeval))
    3065     else:
    3066           print errormsg
    3067           print '    timeref_datetime: time units "' + tu + '" not ready!!!!'
    3068           quit(-1)
    3069 
    3070     realdates[0] = int(realdate.year)
    3071     realdates[1] = int(realdate.month)
    3072     realdates[2] = int(realdate.day)
    3073     realdates[3] = int(realdate.hour)
    3074     realdates[4] = int(realdate.second)
    3075     realdates[5] = int(realdate.minute)
    3076 
    3077     return realdates
    3078 
    3079 def realdatetime_CFcompilant(times, Srefdate, tunits):
    3080     """ Function to transform a matrix with real time values ([year, month, day, hour, minute, second]) to a netCDF one
    3081     times= matrix with times
    3082     Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
    3083     tunits= units of time respect to Srefdate
    3084     >>> realdatetime_CFcompilant(np.array([ [1976, 2, 17, 8, 20, 0], [1976, 2, 18, 8, 20, 0]], dtype=int), '19491201000000', 'hours')
    3085     [ 229784.33333333  229808.33333333]
    3086     """
    3087 
    3088     import datetime as dt
    3089     yrref=int(Srefdate[0:4])
    3090     monref=int(Srefdate[4:6])
    3091     dayref=int(Srefdate[6:8])
    3092     horref=int(Srefdate[8:10])
    3093     minref=int(Srefdate[10:12])
    3094     secref=int(Srefdate[12:14])
    3095  
    3096     refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref)
    3097 
    3098     dimt=times.shape[0]
    3099        
    3100     cfdates = np.zeros((dimt), dtype=np.float64)
    3101     if tunits == 'weeks':
    3102         for it in range(dimt):
    3103             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3104             cfdates[it] = (cfdate.days + cfdate.seconds/(3600.*24.))/7.
    3105     elif tunits == 'days':
    3106         for it in range(dimt):
    3107             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3108             cfdates[it] = cfdate.days + cfdate.seconds/(3600.*24.)
    3109     elif tunits == 'hours':
    3110         for it in range(dimt):
    3111             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3112             cfdates[it] = cfdate.days*24. + cfdate.seconds/3600.
    3113     elif tunits == 'minutes':
    3114         for it in range(dimt):
    3115             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3116             cfdates[it] = cfdate.days*24.*60. + cfdate.seconds/60.
    3117     elif tunits == 'seconds':
    3118         for it in range(dimt):
    3119             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3120             cfdates[it] = cfdate.days*24.*3600. + cfdate.seconds
    3121     elif tunits == 'milliseconds':
    3122         for it in range(dimt):
    3123             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3124             cfdates[it] = cfdate.days*1000.*24.*3600. + cfdate.seconds*1000.
    3125     elif tunits == 'microseconds':
    3126         for it in range(dimt):
    3127             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    3128             cfdates[it] = cfdate.days*1000000.*24.*3600. + cfdate.seconds*1000000.
    3129     else:
    3130         print errormsg
    3131         print '  ' + fname + ': time units "' + tunits + '" is not ready!!!'
    3132         quit(-1)
    3133 
    3134     return cfdates
    3135 
    3136 def netCDFdatetime_realdatetime(units, tcalendar, times):
    3137     """ Function to transfrom from netCDF CF-compilant times to real time
    3138     """
    3139     import datetime as dt
    3140 
    3141     txtunits = units.split(' ')
    3142     tunits = txtunits[0]
    3143     Srefdate = txtunits[len(txtunits) - 1]
    3144 
    3145 # Calendar type
    3146 ##
    3147     is360 = False
    3148     if tcalendar is not None:
    3149       print '  netCDFdatetime_realdatetime: There is a calendar attribute'
    3150       if tcalendar == '365_day' or tcalendar == 'noleap':
    3151           print '    netCDFdatetime_realdatetime: No leap years!'
    3152           isleapcal = False
    3153       elif tcalendar == 'proleptic_gregorian' or tcalendar == 'standard' or tcalendar == 'gregorian':
    3154           isleapcal = True
    3155       elif tcalendar == '360_day':
    3156           is360 = True
    3157           isleapcal = False
    3158       else:
    3159           print errormsg
    3160           print '    netCDFdatetime_realdatetime: Calendar "' + tcalendar + '" not prepared!'
    3161           quit(-1)
    3162 
    3163 # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
    3164 ##
    3165     timeval = Srefdate.find(':')
    3166 
    3167     if not timeval == -1:
    3168         print '  netCDFdatetime_realdatetime: refdate with time!'
    3169         refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    3170     else:
    3171         refdate = dateStr_date(Srefdate)
    3172 
    3173     dimt = len(times)
    3174 #    datetype = type(dt.datetime(1972,02,01))
    3175 #    realdates = np.array(dimt, datetype)
    3176 #    print realdates
    3177 
    3178 ## Not in timedelta
    3179 #  if tunits == 'years':
    3180 #    for it in range(dimt):
    3181 #      realdate = refdate + dt.timedelta(years=float(times[it]))
    3182 #      realdates[it] = int(realdate.year)
    3183 #  elif tunits == 'months':
    3184 #    for it in range(dimt):
    3185 #      realdate = refdate + dt.timedelta(months=float(times[it]))
    3186 #      realdates[it] = int(realdate.year)
    3187 #    realdates = []
    3188     realdates = np.zeros((dimt, 6), dtype=int)
    3189     if tunits == 'weeks':
    3190         for it in range(dimt):
    3191             realdate = refdate + dt.timedelta(weeks=float(times[it]))
    3192             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3193     elif tunits == 'days':
    3194         for it in range(dimt):
    3195             realdate = refdate + dt.timedelta(days=float(times[it]))
    3196             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3197     elif tunits == 'hours':
    3198         for it in range(dimt):
    3199             realdate = refdate + dt.timedelta(hours=float(times[it]))
    3200 #            if not isleapcal:
    3201 #                Nleapdays = cal.leapdays(int(refdate.year), int(realdate.year))
    3202 #                realdate = realdate - dt.timedelta(days=Nleapdays)
    3203 #            if is360:
    3204 #                Nyears360 = int(realdate.year) - int(refdate.year) + 1
    3205 #                realdate = realdate -dt.timedelta(days=Nyears360*5)
    3206 #            realdates[it] = realdate
    3207 #        realdates = refdate + dt.timedelta(hours=float(times))
    3208             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3209     elif tunits == 'minutes':
    3210         for it in range(dimt):
    3211             realdate = refdate + dt.timedelta(minutes=float(times[it]))
    3212             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3213     elif tunits == 'seconds':
    3214         for it in range(dimt):
    3215             realdate = refdate + dt.timedelta(seconds=float(times[it]))
    3216             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3217     elif tunits == 'milliseconds':
    3218         for it in range(dimt):
    3219             realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
    3220             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3221     elif tunits == 'microseconds':
    3222         for it in range(dimt):
    3223             realdate = refdate + dt.timedelta(microseconds=float(times[it]))
    3224             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    3225     else:
    3226         print errormsg
    3227         print '  netCDFdatetime_realdatetime: time units "' + tunits + '" is not ready!!!'
    3228         quit(-1)
    3229 
    3230     return realdates
    32311591
    32321592class cls_time_information(object):
    … …  
    34331793    return timeinf
    34341794
    3435 def date_juliandate(refyr, TOTsecs):
    3436     """ Function to transform from a total quantity of seconds since the
    3437       beginning of a year to 360 days / year calendar
    3438       TOTsecs= total number of seconds
    3439     """
    3440     fname = 'date_juliandate'
    3441 
    3442     secsYear = 3600*24*30*12.
    3443     secsMonth = 3600*24*30.
    3444     secsDay = 3600*24.
    3445     secsHour = 3600.
    3446     secsMinute = 60.
    3447 
    3448     rldate = np.zeros((6), dtype=int)
    3449     rldate[0] = refyr
    3450     if TOTsecs > secsYear:
    3451         rldate[0] = refyr + int(TOTsecs/secsYear) + 1
    3452         TOTsecs = TOTsecs - int(TOTsecs/secsYear)*secsYear
    3453 
    3454     if np.mod(TOTsecs,secsMonth) == 0:
    3455         rldate[1] = int(TOTsecs/secsMonth)
    3456     else:
    3457         rldate[1] = int(TOTsecs/secsMonth) + 1
    3458     TOTsecs = TOTsecs - int(TOTsecs/secsMonth)*secsMonth
    3459 
    3460     if np.mod(TOTsecs,secsDay) == 0:
    3461         rldate[2] = int(TOTsecs/secsDay)
    3462     else:
    3463         rldate[2] = int(TOTsecs/secsDay) + 1
    3464     TOTsecs = TOTsecs - int(TOTsecs/secsDay)*secsDay
    3465 
    3466     if np.mod(TOTsecs,secsHour) == 0:
    3467         rldate[3] = int(TOTsecs/secsHour)
    3468     else:
    3469         rldate[3] = int(TOTsecs/secsHour) + 1
    3470     TOTsecs = TOTsecs - int(TOTsecs/secsHour)*secsHour
    3471 
    3472     if np.mod(TOTsecs,secsMinute) == 0:
    3473         rldate[4] = int(TOTsecs/secsMinute)
    3474     else:
    3475         rldate[4] = int(TOTsecs/secsMinute) + 1
    3476     TOTsecs = TOTsecs - int(TOTsecs/secsMinute)*secsMinute
    3477 
    3478     rldate[5] = TOTsecs
    3479 
    3480 #    print refyr,TOTsecs,':',rldate
    3481 #    quit()
    3482 
    3483     return rldate
    3484 
    3485 def CFtimes_datetime(ncfile, tname):
    3486     """ Provide date/time array from a file with a series of netCDF CF-compilant time variable
    3487     ncfile = netCDF file name
    3488     tname = name of the variable time in [ncfile]
    3489     output:
    3490       array(dimt, 0) = year
    3491       array(dimt, 1) = month
    3492       array(dimt, 2) = day
    3493       array(dimt, 3) = hour
    3494       array(dimt, 4) = minute
    3495       array(dimt, 5) = second
    3496     """
    3497     import datetime as dt
    3498     fname = 'CFtimes_datetime'
    3499 
    3500     times = ncfile.variables[tname]
    3501     timevals = times[:]
    3502 
    3503     attvar = times.ncattrs()
    3504     if not searchInlist(attvar, 'units'):
    3505         print errormsg
    3506         print '  ' + fname + ": '" + tname + "' does not have attribute: 'units'"
    3507         quit(-1)
    3508     else:
    3509         units = times.getncattr('units')
    3510  
    3511     txtunits = units.split(' ')
    3512     tunits = txtunits[0]
    3513     Srefdate = txtunits[len(txtunits) - 1]
    3514 # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
    3515 ##
    3516     timeval = Srefdate.find(':')
    3517 
    3518     if not timeval == -1:
    3519 #        print '  refdate with time!'
    3520         refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    3521     else:
    3522         refdate = dateStr_date(Srefdate)
    3523 
    3524     dimt = len(timevals)
    3525     realdates = np.zeros((dimt, 6), dtype=int)
    3526 
    3527     secsDay=3600*24.
    3528 
    3529 # Checking calendar!
    3530 ##
    3531     y360 = False
    3532     daycal360 = ['earth_360d', '360d', '360days', '360_day']
    3533     if searchInlist(attvar, 'calendar'):
    3534         calendar = times.getncattr('calendar')
    3535         if searchInlist(daycal360,calendar):
    3536             print warnmsg
    3537             print '  ' + fname + ': calendar of 12 months of 30 days !!'
    3538             y360 = True
    3539 
    3540 ## Not in timedelta
    3541 #    if tunits == 'years':
    3542 #        for it in range(dimt):
    3543 #            realdate = refdate + dt.timedelta(years=float(times[it]))
    3544 #            realdates[it] = int(realdate.year)
    3545 #    elif tunits == 'months':
    3546 #        for it in range(dimt):
    3547 #            realdate = refdate + dt.timedelta(months=float(times[it]))
    3548 #            realdates[it] = int(realdate.year)
    3549     if y360:
    3550         if tunits == 'weeks':
    3551             for it in range(dimt):
    3552                 deltat = dt.timedelta(weeks=float(times[it]))
    3553                 Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
    3554                 realdates[it,:] = date_juliandate(refdate.year,Tsecs)
    3555         elif tunits == 'days':
    3556             for it in range(dimt):
    3557                 deltat = dt.timedelta(days=float(times[it]))
    3558                 Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
    3559                 realdates[it,:] = date_juliandate(refdate.year,Tsecs)
    3560         elif tunits == 'hours':
    3561            for it in range(dimt):
    3562                 realdate = dt.timedelta(hours=float(times[it]))
    3563                 Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
    3564                 realdates[it,:] = date_juliandate(refdate.year,Tsecs)
    3565         elif tunits == 'minutes':
    3566            for it in range(dimt):
    3567                 realdate = dt.timedelta(minutes=float(times[it]))
    3568                 Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
    3569                 realdates[it,:] = date_juliandate(refdate.year,Tsecs)
    3570         elif tunits == 'seconds':
    3571            for it in range(dimt):
    3572                 realdate = dt.timedelta(seconds=float(times[it]))
    3573                 Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
    3574                 realdates[it,:] = date_juliandate(refdate.year,Tsecs)
    3575         elif tunits == 'miliseconds':
    3576            for it in range(dimt):
    3577                 realdate = dt.timedelta(miliseconds=float(times[it]))
    3578                 Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
    3579                 realdates[it,:] = date_juliandate(refdate.year,Tsecs)
    3580         else:
    3581               print errormsg
    3582               print '    CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
    3583               quit(-1)
    3584     else:
    3585         if tunits == 'weeks':
    3586             for it in range(dimt):
    3587                 realdate = refdate + dt.timedelta(weeks=float(times[it]))
    3588                 realdates[it,0] = int(realdate.year)
    3589                 realdates[it,1] = int(realdate.month)
    3590                 realdates[it,2] = int(realdate.day)
    3591                 realdates[it,3] = int(realdate.hour)
    3592                 realdates[it,4] = int(realdate.second)
    3593                 realdates[it,5] = int(realdate.minute)
    3594         elif tunits == 'days':
    3595             for it in range(dimt):
    3596                 realdate = refdate + dt.timedelta(days=float(times[it]))
    3597                 realdates[it,0] = int(realdate.year)
    3598                 realdates[it,1] = int(realdate.month)
    3599                 realdates[it,2] = int(realdate.day)
    3600                 realdates[it,3] = int(realdate.hour)
    3601                 realdates[it,4] = int(realdate.second)
    3602                 realdates[it,5] = int(realdate.minute)
    3603         elif tunits == 'hours':
    3604            for it in range(dimt):
    3605                 realdate = refdate + dt.timedelta(hours=float(times[it]))
    3606                 realdates[it,0] = int(realdate.year)
    3607                 realdates[it,1] = int(realdate.month)
    3608                 realdates[it,2] = int(realdate.day)
    3609                 realdates[it,3] = int(realdate.hour)
    3610                 realdates[it,4] = int(realdate.second)
    3611                 realdates[it,5] = int(realdate.minute)
    3612         elif tunits == 'minutes':
    3613            for it in range(dimt):
    3614                 realdate = refdate + dt.timedelta(minutes=float(times[it]))
    3615                 realdates[it,0] = int(realdate.year)
    3616                 realdates[it,1] = int(realdate.month)
    3617                 realdates[it,2] = int(realdate.day)
    3618                 realdates[it,3] = int(realdate.hour)
    3619                 realdates[it,4] = int(realdate.second)
    3620                 realdates[it,5] = int(realdate.minute)
    3621         elif tunits == 'seconds':
    3622            for it in range(dimt):
    3623                 realdate = refdate + dt.timedelta(seconds=float(times[it]))
    3624                 realdates[it,0] = int(realdate.year)
    3625                 realdates[it,1] = int(realdate.month)
    3626                 realdates[it,2] = int(realdate.day)
    3627                 realdates[it,3] = int(realdate.hour)
    3628                 realdates[it,4] = int(realdate.second)
    3629                 realdates[it,5] = int(realdate.minute)
    3630         elif tunits == 'milliseconds':
    3631            for it in range(dimt):
    3632                 realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
    3633                 realdates[it,0] = int(realdate.year)
    3634                 realdates[it,1] = int(realdate.month)
    3635                 realdates[it,2] = int(realdate.day)
    3636                 realdates[it,3] = int(realdate.hour)
    3637                 realdates[it,4] = int(realdate.second)
    3638                 realdates[it,5] = int(realdate.minute)
    3639         else:
    3640               print errormsg
    3641               print '    CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
    3642               quit(-1)
    3643    
    3644     return realdates
    3645 
    36461795class variable_inf(object):
    36471796    """ Class to provide the information from a given variable
    … …  
    37271876                self.dimy=self.dims[2]
    37281877                self.dimx=self.dims[3]
     1878
     1879    return
    37291880
    37301881def subyrs(values, ncfile, varn):
    … …  
    39152066  print '      subyrs: File "' + ofile + '" with a subset of ' + values + ' has been created'
    39162067
     2068  return
     2069
    39172070def submns(values, ncfile, varn):
    39182071  """ Function to retrieve a series of months from a file
    … …  
    41032256  print '      submns: File "' + ofile + '" with a subset of ' + values + ' has been created'
    41042257
    4105 class statsValWeigthed(object):
    4106   """Weigthed Statistics class providing:
    4107   vals = values (can be a matrix)
    4108   wgs = weights (can be a matrix)
    4109   self.meanv: mean weigthed value
    4110   self.mean2v: mean quadratic weigthed value
    4111   self.stdv: weigthed standard deviation
    4112   self.Nokvalue non None values of a list of values
    4113   self.meanwgt: mean of the weigths
    4114   self.mean2wgt: cuadratic mean of the weigths
    4115   self.stdwgt: standard deviation of the weigths
    4116   """
    4117 
    4118   def __init__(self, vals, wgs):
    4119     if vals is None:
    4120       self.Nv = None
    4121       self.meanv = None
    4122       self.mean2v = None
    4123       self.stdv = None
    4124       self.Nokvalues = None
    4125       self.meanwgt = None
    4126       self.mean2wgt = None
    4127       self.stdwgt = None
    4128     else:
    4129       values = vals.flatten()
    4130       weights = wgs.flatten()
    4131       self.Nv=len(values)
    4132       self.meanv=0.
    4133       self.mean2v=0.
    4134       self.stdv=0.
    4135       self.meanwgt = 0.
    4136       self.mean2wgt = 0.
    4137       self.stdwgt = 0.
    4138       self.Nokvalues = 0
    4139 
    4140       for inum in range(self.Nv):
    4141         if not values[inum] is None:
    4142           self.Nokvalues = self.Nokvalues + 1
    4143           self.meanv = self.meanv+values[inum]*weights[inum]
    4144           self.mean2v = self.mean2v+values[inum]*weights[inum]*values[inum]
    4145           self.meanwgt = self.meanwgt+weights[inum]
    4146           self.mean2wgt = self.mean2wgt+weights[inum]*weights[inum]
    4147 
    4148       self.meanv = self.meanv/float(self.meanwgt)
    4149       self.mean2v = self.mean2v/float(self.meanwgt)
    4150       self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
    4151       self.meanwgt = self.meanwgt/float(self.Nokvalues)
    4152       self.mean2wgt = self.mean2wgt/float(self.Nokvalues)
    4153       self.stdwgt = np.sqrt(self.mean2wgt-self.meanwgt*self.meanwgt)
    4154 
    4155 class statsValWeighted_missVal(object):
    4156   """Weighted Statistics taking into account a missing value class providing:
    4157   vals = values (can be a matrix)
    4158   wgs = weights (can be a matrix)
    4159   missVal= missing value
    4160   self.meanv: mean weigthed value
    4161   self.mean2v: mean quadratic weigthed value
    4162   self.stdv: weigthed standard deviation
    4163   self.Nokvalue non None values of a list of values
    4164   self.meanwgt: mean of the weigths
    4165   self.mean2wgt: cuadratic mean of the weigths
    4166   self.stdwgt: standard deviation of the weigths
    4167   self.quantilesv: quantiles of the weighted values
    4168   """
    4169 
    4170   def __init__(self, vals, wgs, missVal):
    4171 
    4172     fname='statsValWeigthed_missVal'
    4173     if vals is None:
    4174       self.Nv = None
    4175       self.meanv = None
    4176       self.mean2v = None
    4177       self.stdv = None
    4178       self.Nokvalues = None
    4179       self.meanwgt = None
    4180       self.mean2wgt = None
    4181       self.stdwgt = None
    4182       self.quantilesv = None
    4183     else:   
    4184       Npt = 1
    4185       for idim in range(len(vals.shape)):
    4186         Npt = Npt * vals.shape[idim]
    4187       if np.sum(vals >= missVal) == Npt:
    4188           print errormsg
    4189           print '  ' + fname + ' all values get missing!!'
    4190           print errormsg
    4191           quit(-1)
    4192       vals0 = np.where(vals >= missVal, None, vals)
    4193       values = vals0.flatten()
    4194       weights = wgs.flatten()
    4195       self.Nv=Npt
    4196       self.meanv=0.
    4197       self.mean2v=0.
    4198       self.stdv=0.
    4199       self.meanwgt = 0.
    4200       self.mean2wgt = 0.
    4201       self.stdwgt = 0.
    4202       self.Nokvalues = 0
    4203 
    4204       valswgt = values
    4205       for inum in range(self.Nv):
    4206         if not values[inum] is None:
    4207           self.Nokvalues = self.Nokvalues + 1
    4208           valswgt[inum] = valswgt[inum]*weights[inum]
    4209           self.meanv = self.meanv+values[inum]*weights[inum]
    4210           self.mean2v = self.mean2v+values[inum]*weights[inum]*values[inum]
    4211           self.meanwgt = self.meanwgt+weights[inum]
    4212           self.mean2wgt = self.mean2wgt+weights[inum]*weights[inum]
    4213 
    4214       self.meanv = self.meanv/float(self.meanwgt)
    4215       self.mean2v = self.mean2v/float(self.meanwgt)
    4216       self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
    4217       valswgt = valswgt/np.float(self.meanwgt)
    4218       self.meanwgt = self.meanwgt/float(self.Nokvalues)
    4219       self.mean2wgt = self.mean2wgt/float(self.Nokvalues)
    4220       self.stdwgt = np.sqrt(self.mean2wgt-self.meanwgt*self.meanwgt)
    4221       valsq=Quantiles(valswgt, 20)
    4222       self.quantilesv=valsq.quantilesv
    4223 
    4224 class stats2Val(object):
    4225   """two variables Statistics class providing:
    4226   vals1 = variable 1
    4227   vals2 = variable 2
    4228   power = power of the polynomial fitting to apply between both variables
    4229   self.min[var], self.max[var], self.mean[var], self.mean2[var], self.std[var] of
    4230     [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]
    4231   self.Nokvalues1: number of correct values of variable 1
    4232   self.Nokvalues2: number of correct values of variable 2
    4233   self.Nokvalues12: number of correct coincident values of variable 1 and variable 2
    4234   self.mae=mean(abs(var1-var2))
    4235   self.rmse=sqrt((var1-var2)**2)
    4236   self.correlation (and p-value)
    4237   self.linRegress: linear regression [trend, intercept, regression coefficient, p_value, standard error]
    4238   self.polRegress: polinomial Regresion  of degree [power] [coef**[power], coef**[power-1], ...., coef**0]
    4239   """
    4240 
    4241   def __init__(self, vals1, vals2, power):
    4242     import numpy as np
    4243     from scipy import stats as sts
    4244 
    4245     if vals1 is None:
    4246       self.Nv = None
    4247       self.Nokvalues1 = None
    4248       self.Nokvalues2 = None
    4249       self.Nokvalues12 = None
    4250       self.NDvalNone = None
    4251       self.minv1Av2 = None
    4252       self.maxv1Av2 = None
    4253       self.meanv1Av2 = None
    4254       self.mean2v1Av2 = None
    4255       self.stdv1Av2 = None
    4256       self.minv1Sv2 = None
    4257       self.maxv1Sv2 = None
    4258       self.meanv1Sv2 = None
    4259       self.mean2v1Sv2 = None
    4260       self.stdv1Sv2 = None
    4261       self.minv1Dv2 = None
    4262       self.maxv1Dv2 = None
    4263       self.meanv1Dv2 = None
    4264       self.mean2v1Dv2 = None
    4265       self.stdv1Dv2 = None
    4266       self.minv1Pv2 = None
    4267       self.maxv1Pv2 = None
    4268       self.meanv1Pv2 = None
    4269       self.mean2v1Pv2 = None
    4270       self.stdv1Pv2 = None
    4271       self.mae = None
    4272       self.rmse = None
    4273       self.corr = None
    4274       self.linRegress = None
    4275       self.polRegress = None
    4276       self.polRegressResidual = None
    4277       self.polRegressRes = None
    4278       self.polRegressSingVal = None
    4279     else:
    4280       values1 = vals1.flatten()
    4281       values2 = vals2.flatten()
    4282 
    4283       if not len(values1) == len(values2):
    4284         print errormsg
    4285         print '    stats2Val: lengths of variables differ!! Lvar1: ', len(values1), ' Lvar2: ',len(values2),' statistics between them can not be computed!'
    4286         quit(-1)
    4287 
    4288       self.Nv=len(values1)
    4289       self.minv1Av2=10000000000.
    4290       self.maxv1Av2=-self.minv1Av2
    4291       self.meanv1Av2=0.
    4292       self.mean2v1Av2=0.
    4293       self.stdv1Av2=0.
    4294       self.minv1Sv2=self.minv1Av2
    4295       self.maxv1Sv2=-self.minv1Av2
    4296       self.meanv1Sv2=0.
    4297       self.mean2v1Sv2=0.
    4298       self.stdv1Sv2=0.
    4299       self.minv1Dv2=self.minv1Av2
    4300       self.maxv1Dv2=-self.minv1Av2
    4301       self.meanv1Dv2=0.
    4302       self.mean2v1Dv2=0.
    4303       self.stdv1Dv2=0.
    4304       self.minv1Pv2=self.minv1Av2
    4305       self.maxv1Pv2=-self.minv1Av2
    4306       self.meanv1Pv2=0.
    4307       self.mean2v1Pv2=0.
    4308       self.stdv1Pv2=0.
    4309       self.mae = 0.
    4310       self.rmse = 0.
    4311       self.corr = np.array([0., 0.])
    4312       self.linRegress = np.zeros(5, float)
    4313       self.polRegress = np.zeros(power+1, float)
    4314       self.polRegressResidual = 0.
    4315       self.polRegressSingVal = np.zeros(power+1, float)
    4316 
    4317 # v1 [+ / - / / / *] v2
    4318 ##
    4319       self.Nokvalues1 = 0
    4320       self.Nokvalues2 = 0
    4321       self.Nokvalues12 = 0
    4322       self.NDvalNone = 0
    4323       for inum in range(self.Nv):
    4324         if not values1[inum] is None:
    4325           self.Nokvalues1 = self.Nokvalues1 + 1
    4326         if not values2[inum] is None:
    4327           self.Nokvalues2 = self.Nokvalues2 + 1
    4328         if not values1[inum] is None and not values2[inum] is None:
    4329           self.Nokvalues12 = self.Nokvalues12 + 1
    4330           Aval = values1[inum] + values2[inum]
    4331           Sval = values1[inum] - values2[inum]
    4332           Pval = values1[inum] * values2[inum]
    4333           if np.isinf(values1[inum] / values2[inum]) or np.isnan(values1[inum] / values2[inum]):
    4334             if self.NDvalNone < 1:
    4335                print warnmsg
    4336                print '      stats2Val: val1/val2 inf or Nan!!!!'
    4337             Dval = None
    4338             self.NDvalNone = self.NDvalNone + 1
    4339           else:
    4340             Dval = values1[inum] / values2[inum]
    4341 
    4342           self.mae = self.mae + abs(Sval)
    4343           self.rmse = self.rmse + Sval**2
    4344 
    4345           if Aval < self.minv1Av2:
    4346             self.minv1Av2 = Aval
    4347           if Aval > self.maxv1Av2:
    4348             self.maxv1Av2 = Aval
    4349           if Sval < self.minv1Sv2:
    4350             self.minv1Sv2 = Sval
    4351           if Sval > self.maxv1Sv2:
    4352             self.maxv1Sv2 = Sval
    4353           if not Dval is None and Dval < self.minv1Dv2:
    4354             self.minv1Dv2 = Dval
    4355           if not Dval is None and  Dval > self.maxv1Dv2:
    4356             self.maxv1Dv2 = Dval
    4357           if Pval < self.minv1Pv2:
    4358             self.minv1Pv2 = Pval
    4359           if Pval > self.maxv1Pv2:
    4360             self.maxv1Pv2 = Pval
    4361 
    4362           self.meanv1Av2 = self.meanv1Av2+Aval
    4363           self.mean2v1Av2 = self.mean2v1Av2+Aval*Aval
    4364           self.meanv1Sv2 = self.meanv1Sv2+Sval
    4365           self.mean2v1Sv2 = self.mean2v1Sv2+Sval*Sval
    4366           if not Dval is None:
    4367             self.meanv1Dv2 = self.meanv1Dv2+Dval
    4368             self.mean2v1Dv2 = self.mean2v1Dv2+Dval*Dval
    4369           self.meanv1Pv2 = self.meanv1Pv2+Pval
    4370           self.mean2v1Pv2 = self.mean2v1Pv2+Pval*Pval
    4371 
    4372 ##      print 'Nokvalues1: ', self.Nokvalues1, 'Nokvalues2: ', self.Nokvalues2, 'Nokvalues12: ', float(self.Nokvalues12), 'NDvalNone: ',self.NDvalNone
    4373       self.meanv1Av2 = self.meanv1Av2/float(self.Nokvalues12)
    4374       self.mean2v1Av2 = self.mean2v1Av2/float(self.Nokvalues12)
    4375       self.stdv1Av2 = np.sqrt(self.mean2v1Av2-self.meanv1Av2*self.meanv1Av2)
    4376       self.meanv1Sv2 = self.meanv1Sv2/float(self.Nokvalues12)
    4377       self.mean2v1Sv2 = self.mean2v1Sv2/float(self.Nokvalues12)
    4378       self.stdv1Sv2 = np.sqrt(self.mean2v1Sv2-self.meanv1Sv2*self.meanv1Sv2)
    4379       if self.Nokvalues12 - self.NDvalNone == 0:
    4380           self.meanv1Dv2 = None
    4381           self.mean2v1Dv2 = None
    4382           self.stdv1Dv2 = None
    4383           print warnmsg
    4384           print '      stats2Val: all values of val1/val2 are None!'
    4385       else:
    4386           self.meanv1Dv2 = self.meanv1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
    4387           self.mean2v1Dv2 = self.mean2v1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
    4388           self.stdv1Dv2 = np.sqrt(self.mean2v1Dv2-self.meanv1Dv2*self.meanv1Dv2)
    4389       self.meanv1Pv2 = self.meanv1Pv2/float(self.Nokvalues12)
    4390       self.mean2v1Pv2 = self.mean2v1Pv2/float(self.Nokvalues12)
    4391       self.stdv1Pv2 = np.sqrt(self.mean2v1Pv2-self.meanv1Pv2*self.meanv1Pv2)
    4392 
    4393       self.mae = self.mae/self.Nokvalues12
    4394       self.rmse = np.sqrt(self.rmse/self.Nokvalues12)
    4395 
    4396       self.corr = sts.pearsonr(values1, values2)
    4397 
    4398       self.linRegress[0], self.linRegress[1], self.linRegress[2], self.linRegress[3], self.linRegress[4] = sts.linregress(values1, values2)
    4399 
    4400       polyfitvals=np.polyfit(values1, values2, power, full = True)
    4401 
    4402       self.polRegress = polyfitvals[0]
    4403       self.polRegressRes = polyfitvals[1]
    4404       self.polRegressSingVal = polyfitvals[3]
    4405 
    4406 class stats2Val_missVal(object):
    4407   """two variables Statistics taking into account a missing value class providing:
    4408   vals1 = variable 1
    4409   vals2 = variable 2
    4410   missVal = missing value
    4411   power = power of the polynomial fitting to apply between both variables
    4412   self.min[var], self.max[var], self.mean[var], self.mean2[var], self.std[var] of
    4413     [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]
    4414   self.Nokvalues1: number of correct values of variable 1
    4415   self.Nokvalues2: number of correct values of variable 2
    4416   self.Nokvalues12: number of correct coincident values of variable 1 and variable 2
    4417   self.mae=mean(abs(var1-var2))
    4418   self.rmse=sqrt((var1-var2)**2)
    4419   self.correlation (and p-value)
    4420   self.linRegress: linear regression [trend, intercept, regression coefficient, p_value, standard error]
    4421   self.polRegress: polinomial Regresion  of degree [power] [coef**[power], coef**[power-1], ...., coef**0]
    4422   """
    4423 
    4424   def __init__(self, vals1, vals2, power, missVal):
    4425     import numpy as np
    4426     from scipy import stats as sts
    4427 
    4428     fname='stats2Val_missVal'
    4429     if vals1 is None:
    4430       self.Nv = None
    4431       self.Nokvalues1 = None
    4432       self.Nokvalues2 = None
    4433       self.Nokvalues12 = None
    4434       self.NDvalNone = None
    4435       self.minv1Av2 = None
    4436       self.maxv1Av2 = None
    4437       self.meanv1Av2 = None
    4438       self.mean2v1Av2 = None
    4439       self.stdv1Av2 = None
    4440       self.minv1Sv2 = None
    4441       self.maxv1Sv2 = None
    4442       self.meanv1Sv2 = None
    4443       self.mean2v1Sv2 = None
    4444       self.stdv1Sv2 = None
    4445       self.minv1Dv2 = None
    4446       self.maxv1Dv2 = None
    4447       self.meanv1Dv2 = None
    4448       self.mean2v1Dv2 = None
    4449       self.stdv1Dv2 = None
    4450       self.minv1Pv2 = None
    4451       self.maxv1Pv2 = None
    4452       self.meanv1Pv2 = None
    4453       self.mean2v1Pv2 = None
    4454       self.stdv1Pv2 = None
    4455       self.mae = None
    4456       self.rmse = None
    4457       self.corr = None
    4458       self.linRegress = None
    4459       self.polRegress = None
    4460       self.polRegressResidual = None
    4461       self.polRegressRes = None
    4462       self.polRegressSingVal = None
    4463     else:
    4464       Npt1 = 1
    4465       for idim in range(len(vals1.shape)):
    4466         Npt1 = Npt1 * vals1.shape[idim]
    4467       if np.sum(vals1 >= missVal) == Npt1:
    4468           print errormsg
    4469           print '  ' + fname + ' all values 1 get missing!!'
    4470           print errormsg
    4471           quit(-1)
    4472       Npt2 = 1
    4473       for idim in range(len(vals2.shape)):
    4474         Npt2 = Npt2 * vals2.shape[idim]
    4475       if np.sum(vals2 >= missVal) == Npt2:
    4476           print errormsg
    4477           print '  ' + fname + ' all values 2 get missing!!'
    4478           print errormsg
    4479           quit(-1)
    4480       vals10 = np.where(abs(vals1) >= missVal, None, vals1)
    4481       vals20 = np.where(abs(vals2) >= missVal, None, vals2)
    4482       values1 = vals10.flatten()
    4483       values2 = vals20.flatten()
    4484 
    4485       if not len(values1) == len(values2):
    4486         print errormsg
    4487         print '    stats2Val: lengths of variables differ!! Lvar1: ', len(values1), ' Lvar2: ',len(values2),' statistics between them can not be computed!'
    4488         quit(-1)
    4489 
    4490       self.Nv=Npt1
    4491       self.minv1Av2=10000000000.
    4492       self.maxv1Av2=-self.minv1Av2
    4493       self.meanv1Av2=0.
    4494       self.mean2v1Av2=0.
    4495       self.stdv1Av2=0.
    4496       self.minv1Sv2=self.minv1Av2
    4497       self.maxv1Sv2=-self.minv1Av2
    4498       self.meanv1Sv2=0.
    4499       self.mean2v1Sv2=0.
    4500       self.stdv1Sv2=0.
    4501       self.minv1Dv2=self.minv1Av2
    4502       self.maxv1Dv2=-self.minv1Av2
    4503       self.meanv1Dv2=0.
    4504       self.mean2v1Dv2=0.
    4505       self.stdv1Dv2=0.
    4506       self.minv1Pv2=self.minv1Av2
    4507       self.maxv1Pv2=-self.minv1Av2
    4508       self.meanv1Pv2=0.
    4509       self.mean2v1Pv2=0.
    4510       self.stdv1Pv2=0.
    4511       self.mae = 0.
    4512       self.rmse = 0.
    4513       self.corr = np.array([0., 0.])
    4514       self.linRegress = np.zeros(5, float)
    4515       self.polRegress = np.zeros(power+1, float)
    4516       self.polRegressResidual = 0.
    4517       self.polRegressSingVal = np.zeros(power+1, float)
    4518 
    4519 # v1 [+ / - / / / *] v2
    4520 ##
    4521       self.Nokvalues1 = 0
    4522       self.Nokvalues2 = 0
    4523       self.Nokvalues12 = 0
    4524       self.NDvalNone = 0
    4525       for inum in range(self.Nv):
    4526         if not values1[inum] is None:
    4527           self.Nokvalues1 = self.Nokvalues1 + 1
    4528         if not values2[inum] is None:
    4529           self.Nokvalues2 = self.Nokvalues2 + 1
    4530         if not values1[inum] is None and not values2[inum] is None:
    4531           self.Nokvalues12 = self.Nokvalues12 + 1
    4532           Aval = values1[inum] + values2[inum]
    4533           Sval = values1[inum] - values2[inum]
    4534           Pval = values1[inum] * values2[inum]
    4535           if np.isinf(values1[inum] / values2[inum]) or np.isnan(values1[inum] / values2[inum]):
    4536             if self.NDvalNone < 1:
    4537                print warnmsg
    4538                print '      stats2Val: val1/val2 inf or Nan!!!!'
    4539             Dval = None
    4540             self.NDvalNone = self.NDvalNone + 1
    4541           else:
    4542             Dval = values1[inum] / values2[inum]
    4543 
    4544           self.mae = self.mae + abs(Sval)
    4545           self.rmse = self.rmse + Sval**2
    4546 
    4547           if Aval < self.minv1Av2:
    4548             self.minv1Av2 = Aval
    4549           if Aval > self.maxv1Av2:
    4550             self.maxv1Av2 = Aval
    4551           if Sval < self.minv1Sv2:
    4552             self.minv1Sv2 = Sval
    4553           if Sval > self.maxv1Sv2:
    4554             self.maxv1Sv2 = Sval
    4555           if not Dval is None and Dval < self.minv1Dv2:
    4556             self.minv1Dv2 = Dval
    4557           if not Dval is None and  Dval > self.maxv1Dv2:
    4558             self.maxv1Dv2 = Dval
    4559           if Pval < self.minv1Pv2:
    4560             self.minv1Pv2 = Pval
    4561           if Pval > self.maxv1Pv2:
    4562             self.maxv1Pv2 = Pval
    4563 
    4564           self.meanv1Av2 = self.meanv1Av2+Aval
    4565           self.mean2v1Av2 = self.mean2v1Av2+Aval*Aval
    4566           self.meanv1Sv2 = self.meanv1Sv2+Sval
    4567           self.mean2v1Sv2 = self.mean2v1Sv2+Sval*Sval
    4568           if not Dval is None:
    4569             self.meanv1Dv2 = self.meanv1Dv2+Dval
    4570             self.mean2v1Dv2 = self.mean2v1Dv2+Dval*Dval
    4571           self.meanv1Pv2 = self.meanv1Pv2+Pval
    4572           self.mean2v1Pv2 = self.mean2v1Pv2+Pval*Pval
    4573 
    4574 ##      print 'Nokvalues1: ', self.Nokvalues1, 'Nokvalues2: ', self.Nokvalues2, 'Nokvalues12: ', float(self.Nokvalues12), 'NDvalNone: ',self.NDvalNone
    4575       self.meanv1Av2 = self.meanv1Av2/float(self.Nokvalues12)
    4576       self.mean2v1Av2 = self.mean2v1Av2/float(self.Nokvalues12)
    4577       self.stdv1Av2 = np.sqrt(self.mean2v1Av2-self.meanv1Av2*self.meanv1Av2)
    4578       self.meanv1Sv2 = self.meanv1Sv2/float(self.Nokvalues12)
    4579       self.mean2v1Sv2 = self.mean2v1Sv2/float(self.Nokvalues12)
    4580       self.stdv1Sv2 = np.sqrt(self.mean2v1Sv2-self.meanv1Sv2*self.meanv1Sv2)
    4581       if self.Nokvalues12 - self.NDvalNone == 0:
    4582           self.meanv1Dv2 = None
    4583           self.mean2v1Dv2 = None
    4584           self.stdv1Dv2 = None
    4585           print warnmsg
    4586           print '      stats2Val: all values of val1/val2 are None!'
    4587       else:
    4588           self.meanv1Dv2 = self.meanv1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
    4589           self.mean2v1Dv2 = self.mean2v1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
    4590           self.stdv1Dv2 = np.sqrt(self.mean2v1Dv2-self.meanv1Dv2*self.meanv1Dv2)
    4591       self.meanv1Pv2 = self.meanv1Pv2/float(self.Nokvalues12)
    4592       self.mean2v1Pv2 = self.mean2v1Pv2/float(self.Nokvalues12)
    4593       self.stdv1Pv2 = np.sqrt(self.mean2v1Pv2-self.meanv1Pv2*self.meanv1Pv2)
    4594 
    4595       self.mae = self.mae/self.Nokvalues12
    4596       self.rmse = np.sqrt(self.rmse/self.Nokvalues12)
    4597 
    4598       vals1Nomiss = np.ones(len(values1), dtype=bool)
    4599       vals2Nomiss = np.ones(len(values1), dtype=bool)
    4600 
    4601       for i in range(len(values1)):
    4602           if values1[i] is None:
    4603               vals1Nomiss[i] = False
    4604 
    4605       for i in range(len(values2)):
    4606           if values2[i] is None:
    4607               vals2Nomiss[i] = False
    4608 
    4609       v1 = np.array(values1[vals1Nomiss], dtype=float)
    4610       v2 = np.array(values2[vals2Nomiss], dtype=float)
    4611 
    4612       if not v1.shape == v2.shape:
    4613           print errormsg
    4614           print '  ' + fname + ': variables without missing values v1: ',v1.shape , ' and v2: ',v2.shape ,' do not have the same shape! '
    4615           print errormsg
    4616           quit(-1)
    4617 
    4618       self.corr = sts.pearsonr(v1, v2)
    4619 
    4620       self.linRegress[0], self.linRegress[1], self.linRegress[2], self.linRegress[3], self.linRegress[4] = sts.linregress(v1, v2)
    4621 
    4622       polyfitvals=np.polyfit(v1, v2, power, full = True)
    4623 
    4624       self.polRegress = polyfitvals[0]
    4625       self.polRegressRes = polyfitvals[1]
    4626       self.polRegressSingVal = polyfitvals[3]
    4627 
    4628 def mask_2masked(vals1, vals2):
    4629     """ Function to provide the boolean matrix (in opposite way as it is in the mask) as combination of mask from to masked matrices
    4630     """
    4631     import numpy.ma as ma
    4632     fname = 'mask_2masked'
    4633 
    4634 #    if len(vals1.shape) != len(vals2.shape):
    4635 #        print errormsg
    4636 #        print '  ' + fname + ' matrix 1 :', len(vals1.shape), ' and matrix 2 ', len(vals2.shape), ' have different size!'
    4637 #        print errormsg
    4638 #        quit(-1)
    4639 
    4640 #    for idim in range(len(vals1.shape)):
    4641 #        if vals1.shape[idim] != vals2.shape[idim]:
    4642 #            print errormsg
    4643 #            print '  ' + fname + ' dimension ', idim,' from matrix 1 :', vals1.shape[idim], ' and matrix 2 ', \
    4644 #                vals2.shape[idim], ' have different size!'
    4645 #            print errormsg
    4646 #            quit(-1)
    4647 
    4648     if type(vals1) == type(ma.array(1)):
    4649         mask1array=np.where(ma.getmaskarray(vals1) == False, True, False)
    4650     else:
    4651         mask1array=np.ones(vals1.shape, dtype=bool)
    4652 
    4653     if type(vals2) == type(ma.array(1)):
    4654         mask2array=np.where(ma.getmaskarray(vals2) == False, True, False)
    4655     else:
    4656         mask2array=np.ones(vals2.shape, dtype=bool)
    4657 
    4658     mask12 = mask1array*mask2array
    4659 
    4660     return mask12
    4661 
    4662 def mask_pearsonr(xvals, yvals):
    4663     """ Function to compute a pearson correlation from mask matrices
    4664     """
    4665     from scipy import stats as sts
    4666     fillVal = 1.e20
    4667     maskxy = mask_2masked(xvals, yvals)
    4668 
    4669     if np.sum(maskxy) > 1:
    4670         pearsonr = sts.pearsonr(xvals[maskxy], yvals[maskxy])
    4671         if np.isnan(pearsonr[0]) or np.isnan(pearsonr[1]):
    4672            pearsonr = ( fillVal, fillVal)
    4673     else:
    4674         pearsonr = (fillVal, fillVal)
    4675 
    4676     return pearsonr
    4677 
    4678 def mask_quantiles(maskmat, Nquants):
    4679     """ Function to provide the quantiles of a masked array 20 for %5 bins (21 in total)
    4680     """
    4681     import numpy.ma as ma
    4682 
    4683     fillValue = 1.e20
    4684 
    4685     sortmat = maskmat.flatten().copy()
    4686     sortmat.sort()
    4687     quants = np.zeros(Nquants+1, dtype=type(maskmat[0]))
    4688     Nv = ma.size(maskmat)
    4689     NoMask=maskmat.count()
    4690 
    4691     if NoMask < Nquants:
    4692         quants[:] = fillValue
    4693     else:
    4694         for iq in range(Nquants):
    4695             quants[iq] = sortmat[int((NoMask-1)*iq/(Nquants))]
    4696 
    4697         quants[Nquants] = sortmat[NoMask-1]
    4698 
    4699     return quants
    4700 
    4701 def percendone(nvals,tot,percen,msg):
    4702     """ Function to provide the percentage of an action across the matrix
    4703     nvals=number of values
    4704     tot=total number of values
    4705     percen=percentage frequency for which the message is wanted
    4706     msg= message
    4707     """
    4708     from sys import stdout
    4709 
    4710     num = int(tot * percen/100)
    4711     if (nvals%num == 0):
    4712         print '\r        ' + msg + '{0:8.3g}'.format(nvals*100./tot) + ' %',
    4713         stdout.flush()
    4714 
    4715     return ''
    4716 
    4717 def mask_linregres(vals1, vals2):
    4718     """ Function to compute a linear regression from masked data
    4719     vals1: x-values for the regresion
    4720     vals2: y-values for the regresion
    4721     """
    4722     import numpy.ma as ma
    4723 
    4724     fname = 'mask_linregres'
    4725 
    4726     missval1 = vals1.get_fill_value()   
    4727     missval2 = vals2.get_fill_value()   
    4728 
    4729     vals10 = np.where(abs(vals1) >= abs(missval1*0.9), None, vals1)
    4730     vals20 = np.where(abs(vals2) >= abs(missval2*0.9), None, vals2)
    4731 
    4732     values1 = vals10.flatten()
    4733     values2 = vals20.flatten()
    4734 
    4735     vals1Nomiss = np.ones(len(values1), dtype=bool)
    4736     vals2Nomiss = np.ones(len(values2), dtype=bool)
    4737 
    4738     for i in range(len(values1)):
    4739         if values1[i] is None:
    4740             vals1Nomiss[i] = False
    4741     for i in range(len(values2)):
    4742         if values2[i] is None:
    4743             vals2Nomiss[i] = False
    4744 
    4745     v1 = np.array(values1[vals1Nomiss], dtype=float)
    4746     v2 = np.array(values2[vals2Nomiss], dtype=float)
    4747 
    4748     if len(v1) != len(v2):
    4749         print errormsg
    4750         print fname + ': length of masked matrices mat1:',len(v1),'and mat2:',len(v2),'does not match!'
    4751         print errormsg
    4752         quit(-1)
    4753 
    4754     linregres = np.array(sts.linregress(v1, v2), dtype= np.float64)
    4755 
    4756     return linregres
    4757 
    4758 def mask_space_stats(maskmat,statsvals,dim):
    4759     """ Function to give back the multi-dimensional statisitcs of a given masked array
    4760     maskmat=multidimensional masked array
    4761     statsvals=[statn]:[values]
    4762       [statn]: statistics to do:
    4763         'quant', quantiles
    4764       [values]: value for the statistics: Nquantiles
    4765     dim= dimension to run the statistics
    4766     """
    4767     from sys import stdout
    4768 
    4769     fname = 'mask_space_stats'
    4770 
    4771     statn=statsvals.split(':')[0]
    4772     if len(statsvals.split(':')) > 1:
    4773         values=statsvals.split(':')[1]
    4774 
    4775     maskshape = maskmat.shape
    4776     Ndims = len(maskshape)
    4777     if statn == 'quant':
    4778         if len(statsvals.split(':')) == 1:
    4779             print errormsg
    4780             print fname + ': statistics "' + statn + '" requires a value!!!'
    4781             print errormsg
    4782             quit(-1)
    4783         Nquants=int(values)
    4784         if Ndims == 2:
    4785             if dim == 0:
    4786                 statval = np.ones((21, maskshape[1]), dtype=np.float64)*fillValue
    4787                 for i in range(maskshape[1]):
    4788                    percendone(i, maskshape[1], 5, 'quantiles')
    4789                    statval[:,i] = mask_quantiles(maskmat[:,i],Nquants)
    4790             if dim == 1:
    4791                 statval = np.ones((21, maskshape[0]), dtype=np.float64)*fillValue
    4792                 for i in range(maskshape[0]):
    4793                    percendone(i, maskshape[0], 5, 'quantiles')
    4794                    statval[:,i] = mask_quantiles(statval[i,:],Nquants)
    4795         elif Ndims == 3:
    4796             if dim == 0:
    4797                 statval = np.ones((21, maskshape[1], maskshape[2]), dtype=np.float64)*fillValue
    4798                 for i in range(maskshape[1]):
    4799                     for j in range(maskshape[2]):
    4800                         percendone(i*maskshape[2] + j, maskshape[1]*maskshape[2], 5, 'quantiles')
    4801                         statval[:,i,j] = mask_quantiles(maskmat[:,i,j],Nquants)
    4802             if dim == 1:
    4803                 statval = np.ones((21, maskshape[0], maskshape[2]), dtype=np.float64)*fillValue
    4804                 for i in range(maskshape[0]):
    4805                     for j in range(maskshape[2]):
    4806                         percendone(i*maskshape[2] + j, maskshape[0]*maskshape[2], 5, 'quantiles')
    4807                         statval[:,i,j] = mask_quantiles(maskmat[i,:,j],Nquants)
    4808             if dim == 2:
    4809                 statval = np.ones((21, maskshape[0], maskshape[1]), dtype=np.float64)*fillValue
    4810                 for i in range(maskshape[0]):
    4811                     for j in range(maskshape[1]):
    4812                         percendone(i*maskshape[1] + j, maskshape[0]*maskshape[1], 5, 'quantiles')
    4813                         statval[:,i,j] = mask_quantiles(maskmat[i,j,:],Nquants)
    4814         elif Ndims == 4:
    4815             if dim == 0:
    4816                 statval = np.ones((21, maskshape[1], maskshape[2], maskshape[3]), dtype=np.float64)*fillValue
    4817                 for i in range(maskshape[1]):
    4818                     for j in range(maskshape[2]):
    4819                         for k in range(maskshape[3]):
    4820                             percendone(i*maskshape[1]*maskshape[2] + j*maskshape[2] + k, maskshape[1]*maskshape[2]*maskshape[3], 5, 'quantiles')
    4821                             statval[:,i,j,k] = mask_quantiles(maskmat[:,i,j,k],Nquants)
    4822             if dim == 1:
    4823                 statval = np.ones((21, maskshape[0], maskshape[2], maskshape[3]), dtype=np.float64)*fillValue
    4824                 for i in range(maskshape[0]):
    4825                     for j in range(maskshape[2]):
    4826                         for k in range(maskshape[3]):
    4827                             percendone(i*maskshape[0]*maskshape[2] + j*maskshape[2] + k, maskshape[0]*maskshape[2]*maskshape[3], 5, 'quantiles')
    4828                             statval[:,i,j,k] = mask_quantiles(maskmat[i,:,j,k],Nquants)
    4829             if dim == 2:
    4830                 statval = np.ones((21, maskshape[0], maskshape[1], maskshape[3]), dtype=np.float64)*fillValue
    4831                 for i in range(maskshape[0]):
    4832                     for j in range(maskshape[1]):
    4833                         for k in range(maskshape[3]):
    4834                             percendone(i*maskshape[0]*maskshape[1] + j*maskshape[1] + k, maskshape[0]*maskshape[1]*maskshape[3], 5, 'quantiles')
    4835                             statval[:,i,j,k] = mask_quantiles(maskmat[i,j,:,k],Nquants)
    4836             if dim == 3:
    4837                 statval = np.ones((21, maskshape[0], maskshape[1], maskshape[2]), dtype=np.float64)*fillValue
    4838                 for i in range(maskshape[0]):
    4839                     for j in range(maskshape[1]):
    4840                         for k in range(maskshape[2]):
    4841                             percendone(i*maskshape[0]*maskshape[1] + j*maskshape[1] + k, maskshape[0]*maskshape[1]*maskshape[2], 5, 'quantiles')
    4842                             statval[:,i,j,k] = mask_quantiles(maskmat[i,j,k,:],Nquants)
    4843         else:
    4844             print errormsg
    4845             print fname + ': size of matrix ', Ndims,'not ready!!!'
    4846             print errormsg
    4847             quit(-1)
    4848 
    4849     else:
    4850         print errormsg
    4851         print fname + ':  statistics "' + statn + '" not ready!!!!'
    4852         print errormsg
    4853         quit(-1)
    4854 
    4855     print stdout.write("\n")
    4856    
    4857     return statval
    4858 
    4859 class statsVal(object):
    4860   """Statistics class providing
    4861   vals = variable
    4862   self.Nv = number of values
    4863   self.minv = minimum value
    4864   self.maxv = maximum value
    4865   self.meanv = mean value
    4866   self.mean2v = cuadratic mean value
    4867   self.stdv = standard deviation value
    4868   self.Nokvalues = number of correct values of variable
    4869   self.quantilesv = quantiles (%5 bins) of the variable
    4870   """
    4871 
    4872   def __init__(self, vals):
    4873     if vals is None:
    4874       self.Nv = None
    4875       self.minv = None
    4876       self.maxv = None
    4877       self.meanv = None
    4878       self.mean2v = None
    4879       self.stdv = None
    4880       self.Nokvalues = None
    4881       self.quantilesv = None
    4882     else:
    4883       values = vals.flatten()
    4884       self.Nv=len(values)
    4885       self.minv=10000000000.
    4886       self.maxv=-100000000.
    4887       self.meanv=0.
    4888       self.mean2v=0.
    4889       self.stdv=0.
    4890 
    4891       sortedvalues = sorted(values)
    4892 
    4893       self.Nokvalues = 0
    4894       for inum in range(self.Nv):
    4895         if not values[inum] is None:
    4896           self.Nokvalues = self.Nokvalues + 1
    4897           if values[inum] < self.minv:
    4898             self.minv = values[inum]
    4899           if values[inum] > self.maxv:
    4900             self.maxv = values[inum]
    4901 
    4902           self.meanv = self.meanv+values[inum]
    4903           self.mean2v = self.mean2v+values[inum]*values[inum]
    4904 
    4905       self.meanv = self.meanv/float(self.Nokvalues)
    4906       self.mean2v = self.mean2v/float(self.Nokvalues)
    4907       self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
    4908       self.quantilesv = []
    4909       for iq in range(20):
    4910         self.quantilesv.append(sortedvalues[int((self.Nv-1)*iq/20)])
    4911 
    4912       self.quantilesv.append(sortedvalues[self.Nv-1])
    4913       self.medianv = self.quantilesv[10]
    4914 
    4915 class Quantiles(object):
    4916     """ Class to provide quantiles from a given arrayof values
    4917     """
    4918 
    4919     def __init__(self, values, Nquants):
    4920         import numpy.ma as ma
    4921 
    4922         if values is None:
    4923             self.quantilesv = None
    4924 
    4925         else:
    4926             self.quantilesv = []
    4927 
    4928             vals0 = values.flatten()
    4929             Nvalues = len(vals0)
    4930             vals = ma.masked_equal(vals0, None)
    4931             Nvals=len(vals.compressed())
    4932 
    4933             sortedvals = sorted(vals.compressed())
    4934             for iq in range(Nquants):
    4935                 self.quantilesv.append(sortedvals[int((Nvals-1)*iq/Nquants)])
    4936 
    4937             self.quantilesv.append(sortedvals[Nvals-1])
    4938 
    4939 class statsVal_missVal(object):
    4940   """Statistics class tacking into account a missing value providing
    4941   vals = variable
    4942   missval = missing value
    4943   self.Nv = number of values
    4944   self.minv = minimum value
    4945   self.maxv = maximum value
    4946   self.meanv = mean value
    4947   self.mean2v = cuadratic mean value
    4948   self.stdv = standard deviation value
    4949   self.Nokvalues = number of correct values of variable
    4950   self.quantilesv = quantiles (%5 bins) of the variable
    4951   """
    4952 
    4953   def __init__(self, vals, missVal):
    4954     fname='statsVal_missVal'
    4955  
    4956     if vals is None:
    4957       self.Nv = None
    4958       self.minv = None
    4959       self.maxv = None
    4960       self.meanv = None
    4961       self.mean2v = None
    4962       self.stdv = None
    4963       self.Nokvalues = None
    4964       self.quantilesv = None
    4965     else:
    4966       Npt = 1
    4967       for idim in range(len(vals.shape)):
    4968         Npt = Npt * vals.shape[idim]
    4969       if np.sum(vals >= missVal) == Npt:
    4970           print errormsg
    4971           print '  ' + fname + ' all values get missing!!'
    4972           print errormsg
    4973           quit(-1)
    4974 
    4975       vals1 = np.where(abs(vals) >= missVal, None, vals)
    4976       vals0 = np.where(np.isnan(vals1), None, vals1)
    4977 
    4978       values = vals0.flatten()
    4979       self.Nv=Npt
    4980       self.minv=10000000000.
    4981       self.maxv=-100000000.
    4982       self.meanv=0.
    4983       self.mean2v=0.
    4984       self.stdv=0.
    4985 
    4986       sortedvalues = sorted(values)
    4987 
    4988       self.Nokvalues = 0
    4989       for inum in range(self.Nv):
    4990         if not values[inum] is None:
    4991           self.Nokvalues = self.Nokvalues + 1
    4992           if values[inum] < self.minv:
    4993             self.minv = values[inum]
    4994           if values[inum] > self.maxv:
    4995             self.maxv = values[inum]
    4996 
    4997           self.meanv = self.meanv+values[inum]
    4998           self.mean2v = self.mean2v+values[inum]*values[inum]
    4999 
    5000       self.meanv = self.meanv/float(self.Nokvalues)
    5001       self.mean2v = self.mean2v/float(self.Nokvalues)
    5002       self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
    5003       self.quantilesv = []
    5004       for iq in range(20):
    5005         self.quantilesv.append(sortedvalues[int((self.Nv-1)*iq/20)])
    5006 
    5007       self.quantilesv.append(sortedvalues[self.Nv-1])
    5008       self.medianv = self.quantilesv[10]
     2258  return
    50092259
    50102260def spacemean(ncfile, varn):
    … …  
    53882638
    53892639  print '    spacemean: File "' + ofile + '" as space mean of "' + varn + '" has been created'
     2640
     2641  return
    53902642
    53912643def timemean(values, ncfile, varn):
    … …  
    58303082  print '    timemean: File "' + ofile + '" as time mean of "' + varn + '" has been created'
    58313083
    5832 def printing_class(classobj):
    5833     """ Function to print all the values of a given class
    5834     """
    5835 
    5836     valscls = vars(classobj)
    5837     for attrcls in valscls:
    5838         print attrcls, ':', valscls[attrcls]
    5839 
    5840 def fmtprinting_class(classobj):
    5841     """ Function to print all the values of a given class
    5842     """
    5843 
    5844     valscls = vars(classobj)
    5845     for attrcls in valscls:
    5846         print '@' + attrcls + '@', ':', valscls[attrcls]
    5847 
    5848     return
     3084  return
    58493085
    58503086def flipdim(values, filename, varn):
    … …  
    60103246    ncf.close()
    60113247
    6012 def cycl_incr(cyclval,cycle,ninc):
    6013     """ Function to increment a cyclic value [cyclval] with a cycle [cycle] a given number [ninc] of times
    6014     >>> cycl_incr(1,4,1)
    6015     2
    6016     >>> cycl_incr(3,4,1)
    6017     0
    6018     >>> cycl_incr(1,4,10)
    6019     3
     3248    return
     3249
     3250def load_ncvariable_lastdims(ncf, varn, prevdims, Nlastdims):
     3251    """ Function to load the last [Nlastdims] dimensions of a variable [varn] from a netCDF object at the other dimensions at [prevdims]
     3252    ncf= netCDF object
     3253    varn= variable name
     3254    prevdims= values at the previous dimensions
     3255    Nlastims= number of last dimensions
    60203256    """
    6021 
    6022     if ninc >= cycle:
    6023         print 'cycl_incr: WARNING -- warning -- WARNING -- warning'
    6024         print '    given increment: ', ninc,' is larger than the cycle !!'
    6025         ninc = ninc - cycle*int(ninc/cycle)
    6026         print '    reducing it to: ', ninc
    6027 
    6028     val=cyclval + ninc
    6029     if val >= cycle:
    6030         val=cyclval + ninc - cycle
    6031 
    6032     return val
    6033 
    6034 def times_4seasons(tvals, integrity):
    6035     """ Function to split a time series in matricial date format ([:,year,mon,day,hour,minute,second]) in the four stations DJF,MAM,JJA,SON
    6036     tvals= matrix withe times as [:,year,mon,day,hour,minute,second]
    6037     integrity= only give values for entire seasons [True, 3 months for the season], or just split the values by seasons [False]
    6038     """
    6039     fillVal=1.e20
    6040 
    6041 #    print tvals
    6042 
    6043     dt=tvals.shape[0]
    6044     seasons=np.ones((dt,4),dtype=bool)
    6045     seasons=seasons*False
    6046 
    6047     monseas=[12,3,6,9]
    6048     firstseas=False
    6049 
    6050     for it in range(dt):
    6051         if not firstseas:
    6052             if integrity:
    6053                 for iseas in range(4):
    6054                     if tvals[it,1] == monseas[iseas]:
    6055                         nseas=iseas
    6056                         seasons[it,nseas]=True
    6057                         firstseas=True
    6058                         begseas=it
    6059             else:
    6060                 for iseas in range(4):
    6061                     for imon in range(3):
    6062                         if tvals[it,1] == cycl_incr(monseas[iseas],12,imon):
    6063                             nseas=iseas
    6064                             seasons[it,nseas]=True
    6065                             firstseas=True                   
    6066         else:
    6067             newseas=cycl_incr(nseas,4,1)
    6068             if tvals[it,1] == monseas[newseas]:
    6069                 seasons[it,newseas] = True
    6070                 nseas=newseas
    6071                 begseas=it
    6072             else:
    6073                 seasons[it,nseas] = True
    6074 
    6075     endseas = it
    6076 ##    print 'Last season: ',nseas,' beginnig: ',begseas,' ending: ',endseas
    6077 
    6078 # Checking integrity of last season (has to have 3 months)
    6079 ##
    6080     if integrity:
    6081         fullseas=True
    6082         Nmon=np.unique(tvals[begseas:endseas+1,1])
    6083         for it in range(begseas,endseas): fullseas=fullseas*seasons[it,nseas]
    6084         if len(Nmon) < 3 or not fullseas:
    6085             seasons[begseas:endseas+1,nseas] = False
    6086 
    6087     return seasons
    6088 
    6089 def realdatetime_CFcompilant(times, Srefdate, tunits):
    6090     """ Function to transform a matrix with real time values ([year, month, day, hour, minute, second]) to a netCDF one
    6091     times= matrix with times
    6092     Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
    6093     tunits= units of time respect to Srefdate
    6094     >>> realdatetime_CFcompilant(np.array([ [1976, 2, 17, 8, 20, 0], [1976, 2, 18, 8, 20, 0]], dtype=int), '19491201000000', 'hours')
    6095     [ 229784.33333333  229808.33333333]
    6096     """
    6097 
    6098     import datetime as dt
    6099     yrref=int(Srefdate[0:4])
    6100     monref=int(Srefdate[4:6])
    6101     dayref=int(Srefdate[6:8])
    6102     horref=int(Srefdate[8:10])
    6103     minref=int(Srefdate[10:12])
    6104     secref=int(Srefdate[12:14])
    6105  
    6106     refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref)
    6107 
    6108     dimt=times.shape[0]
    6109        
    6110     cfdates = np.zeros((dimt), dtype=np.float64)
    6111     if tunits == 'weeks':
    6112         for it in range(dimt):
    6113             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6114             cfdates[it] = (cfdate.days + cfdate.seconds/(3600.*24.))/7.
    6115     elif tunits == 'days':
    6116         for it in range(dimt):
    6117             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6118             cfdates[it] = cfdate.days + cfdate.seconds/(3600.*24.)
    6119     elif tunits == 'hours':
    6120         for it in range(dimt):
    6121             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6122             cfdates[it] = cfdate.days*24. + cfdate.seconds/3600.
    6123     elif tunits == 'minutes':
    6124         for it in range(dimt):
    6125             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6126             cfdates[it] = cfdate.days*24.*60. + cfdate.seconds/60.
    6127     elif tunits == 'seconds':
    6128         for it in range(dimt):
    6129             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6130             cfdates[it] = cfdate.days*24.*3600. + cfdate.seconds
    6131     elif tunits == 'milliseconds':
    6132         for it in range(dimt):
    6133             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6134             cfdates[it] = cfdate.days*1000.*24.*3600. + cfdate.seconds*1000.
    6135     elif tunits == 'microseconds':
    6136         for it in range(dimt):
    6137             cfdate = dt.datetime(times[it,0], times[it,1], times[it,2], times[it,3], times[it,4], times[it,5]) - refdate
    6138             cfdates[it] = cfdate.days*1000000.*24.*3600. + cfdate.seconds*1000000.
    6139     else:
    6140         print errormsg
    6141         print '  ' + fname + ': time units "' + tunits + '" is not ready!!!'
    6142         quit(-1)
    6143 
    6144     return cfdates
    6145 
    6146 def realdatetime1_CFcompilant(time, Srefdate, tunits):
    6147     """ Function to transform a matrix with a real time value ([year, month, day,
    6148       hour, minute, second]) to a netCDF one
    6149         time= matrix with time
    6150         Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
    6151         tunits= units of time respect to Srefdate
    6152     >>> realdatetime1_CFcompilant([1976, 2, 17, 8, 20, 0], '19491201000000', 'hours')
    6153     229784.33333333
    6154     """
    6155 
    6156     import datetime as dt
    6157     yrref=int(Srefdate[0:4])
    6158     monref=int(Srefdate[4:6])
    6159     dayref=int(Srefdate[6:8])
    6160     horref=int(Srefdate[8:10])
    6161     minref=int(Srefdate[10:12])
    6162     secref=int(Srefdate[12:14])
    6163  
    6164     refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref)
    6165 
    6166     if tunits == 'weeks':
    6167         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5])-refdate
    6168         cfdates = (cfdate.days + cfdate.seconds/(3600.*24.))/7.
    6169     elif tunits == 'days':
    6170         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
    6171         cfdates = cfdate.days + cfdate.seconds/(3600.*24.)
    6172     elif tunits == 'hours':
    6173         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
    6174         cfdates = cfdate.days*24. + cfdate.seconds/3600.
    6175     elif tunits == 'minutes':
    6176         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
    6177         cfdates = cfdate.days*24.*60. + cfdate.seconds/60.
    6178     elif tunits == 'seconds':
    6179         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
    6180         cfdates = cfdate.days*24.*3600. + cfdate.seconds
    6181     elif tunits == 'milliseconds':
    6182         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
    6183         cfdates = cfdate.days*1000.*24.*3600. + cfdate.seconds*1000.
    6184     elif tunits == 'microseconds':
    6185         cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],times[5]) - refdate
    6186         cfdates = cfdate.days*1000000.*24.*3600. + cfdate.seconds*1000000.
    6187     else:
    6188         print errormsg
    6189         print '  ' + fname + ': time units "' + tunits + '" is not ready!!!'
    6190         quit(-1)
    6191 
    6192     return cfdates
    6193 
    6194 def netCDFdatetime_realdatetime(units, tcalendar, times):
    6195     """ Function to transfrom from netCDF CF-compilant times to real time
    6196     [units]= CF time units [tunits] since [YYYY]-[MM]-[HH] [[HH]:[MI]:[SS]]
    6197     [tcalendar]= time calendar
    6198     [times]= CF time values
    6199     """
    6200     import datetime as dt
    6201 
    6202     txtunits = units.split(' ')
    6203     tunits = txtunits[0]
    6204     Srefdate = txtunits[len(txtunits) - 1]
    6205 
    6206 # Calendar type
    6207 ##
    6208     is360 = False
    6209     if tcalendar is not None:
    6210       print '  netCDFdatetime_realdatetime: There is a calendar attribute'
    6211       if tcalendar == '365_day' or tcalendar == 'noleap':
    6212           print '    netCDFdatetime_realdatetime: No leap years!'
    6213           isleapcal = False
    6214       elif tcalendar == 'proleptic_gregorian' or tcalendar == 'standard' or tcalendar == 'gregorian':
    6215           isleapcal = True
    6216       elif tcalendar == '360_day':
    6217           is360 = True
    6218           isleapcal = False
    6219       else:
    6220           print errormsg
    6221           print '    netCDFdatetime_realdatetime: Calendar "' + tcalendar + '" not prepared!'
    6222           quit(-1)
    6223 
    6224 # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
    6225 ##
    6226     timeval = Srefdate.find(':')
    6227 
    6228     if not timeval == -1:
    6229         print '  netCDFdatetime_realdatetime: refdate with time!'
    6230         refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    6231     else:
    6232         refdate = dateStr_date(Srefdate)
    6233 
    6234     dimt = len(times)
    6235 #    datetype = type(dt.datetime(1972,02,01))
    6236 #    realdates = np.array(dimt, datetype)
    6237 #    print realdates
    6238 
    6239 ## Not in timedelta
    6240 #  if tunits == 'years':
    6241 #    for it in range(dimt):
    6242 #      realdate = refdate + dt.timedelta(years=float(times[it]))
    6243 #      realdates[it] = int(realdate.year)
    6244 #  elif tunits == 'months':
    6245 #    for it in range(dimt):
    6246 #      realdate = refdate + dt.timedelta(months=float(times[it]))
    6247 #      realdates[it] = int(realdate.year)
    6248 #    realdates = []
    6249     realdates = np.zeros((dimt, 6), dtype=int)
    6250     if tunits == 'weeks':
    6251         for it in range(dimt):
    6252             realdate = refdate + dt.timedelta(weeks=float(times[it]))
    6253             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6254     elif tunits == 'days':
    6255         for it in range(dimt):
    6256             realdate = refdate + dt.timedelta(days=float(times[it]))
    6257             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6258     elif tunits == 'hours':
    6259         for it in range(dimt):
    6260             realdate = refdate + dt.timedelta(hours=float(times[it]))
    6261 #            if not isleapcal:
    6262 #                Nleapdays = cal.leapdays(int(refdate.year), int(realdate.year))
    6263 #                realdate = realdate - dt.timedelta(days=Nleapdays)
    6264 #            if is360:
    6265 #                Nyears360 = int(realdate.year) - int(refdate.year) + 1
    6266 #                realdate = realdate -dt.timedelta(days=Nyears360*5)
    6267 #            realdates[it] = realdate
    6268 #        realdates = refdate + dt.timedelta(hours=float(times))
    6269             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6270     elif tunits == 'minutes':
    6271         for it in range(dimt):
    6272             realdate = refdate + dt.timedelta(minutes=float(times[it]))
    6273             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6274     elif tunits == 'seconds':
    6275         for it in range(dimt):
    6276             realdate = refdate + dt.timedelta(seconds=float(times[it]))
    6277             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6278     elif tunits == 'milliseconds':
    6279         for it in range(dimt):
    6280             realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
    6281             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6282     elif tunits == 'microseconds':
    6283         for it in range(dimt):
    6284             realdate = refdate + dt.timedelta(microseconds=float(times[it]))
    6285             realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    6286     else:
    6287         print errormsg
    6288         print '  netCDFdatetime_realdatetime: time units "' + tunits + '" is not ready!!!'
    6289         quit(-1)
    6290 
    6291     return realdates
    6292 
    6293 class cls_yearly_means(object):
    6294     """ Class to compute the yearly mean of a time series of values
    6295       tvals=time values (year, month, day, hour, minute, second) format
    6296       values=time-series of values
    6297       dateref= reference date [YYYY][MM][DD][HH][MI][SS] format
    6298       tunits= time units
    6299       self.dimt = Number of values
    6300       self.yrv = yearly mean values
    6301       self.yrt = years of the mean values
    6302     >>> timesv = netCDFdatetime_realdatetime('days since 1949-12-01 00:00:00', 'standard', np.arange(1129))
    6303     >>> valuesv = np.zeros((1129), dtype=np.float)
    6304     >>> valuesv[0:31] = 0.
    6305     >>> valuesv[31:394] = 1.
    6306     >>> valuesv[395:761] = 2.
    6307     >>> valuesv[761:1127] = 3.
    6308     >>> valuesv[1127:1129] = 4.
    6309 
    6310     >>> yrlmeans = cls_yearly_means(timesv, valuesv, '19491201000000', 'days')
    6311     >>> print yrlmeans.dimt,yrlmeans.yrv, yrlmeans.yrt
    6312     5 [ 0.  1.  2.  3.  4.] [ 1949.  1950.  1951.  1952.  1953.]
    6313     """
    6314     def __init__(self, tvals, values, dateref, tunits):
    6315 
    6316         import numpy.ma as ma
    6317         fillVal=1.e20
    6318         fname = 'cls_yearly_means'
    6319 
    6320         if tvals is None:
    6321             self.dimt = None
    6322             self.yrv = None
    6323             self.yrt = None
    6324 
    6325         elif tvals == 'h':
    6326             print fname + '_____________________________________________________________'
    6327             print cls_yearly_means.__doc__
    6328             quit()
    6329 
    6330         else:
    6331             dt=len(values)
    6332             if dt > 0:
    6333                 years = np.unique(tvals[:,0])
    6334                 Nyrs = len(years)
    6335                 yrmean = np.zeros((Nyrs), dtype=np.float)
    6336                 yrt = np.zeros((Nyrs), dtype=np.float)
    6337 
    6338                 iiyr = tvals[0,0]
    6339                 yrmean[0] = values[0]
    6340                 yrt[0] = iiyr
    6341                 iit = 0
    6342                 for iyr in range(Nyrs):
    6343                     nvals = 1
    6344                     for it in range(iit+1,dt):
    6345 #                        print iyr, iiyr, it, tvals[it,0],':',values[it],'->',yrmean[iyr]
    6346                         if tvals[it,0] == iiyr:
    6347                             yrmean[iyr] = yrmean[iyr] + values[it]
    6348                             nvals = nvals + 1
    6349                         else:
    6350                             yrmean[iyr] = yrmean[iyr] / (nvals*1.)
    6351                             iiyr = tvals[it,0]
    6352                             yrmean[iyr + 1] = values[it]
    6353                             yrt[iyr+1] = iiyr
    6354                             iit = it
    6355                             break
    6356 
    6357                 yrmean[Nyrs-1] = yrmean[Nyrs-1]/nvals
    6358                 self.dimt = Nyrs
    6359                 self.yrv = yrmean
    6360                 self.yrt = yrt
    6361             else:
    6362                 print errormsg
    6363                 print '  ' + fname + ': No values passed!'
    6364                 print '    values:', values
    6365                 quit(-1)
    6366 
    6367 class cls_seasonal_means(object):
    6368     """ Class to compute the seasonal mean of a time series of values
    6369     tvals=time values (year, month, day, hour, minute, second) format
    6370     values=time-series of values
    6371     dateref= reference date [YYYY][MM][DD][HH][MI][SS] format
    6372     tunits= time units
    6373     self.dimt = Number of values
    6374     self.seasv = seasonal mean values
    6375     self.seast = seasonal time mean values
    6376     """
    6377     def __init__(self, tvals, values, dateref, tunits):
    6378 
    6379         import numpy.ma as ma
    6380         fillVal=1.e20
    6381         fname = 'cls_seasonal_means'
    6382 
    6383         if tvals is None:
    6384             self.dimt = None
    6385             self.seasv = None
    6386             self.seast = None
    6387 
    6388         else:
    6389             tseasons=times_4seasons(tvals, True)
    6390 
    6391             dt=len(values)
    6392             seasvals=np.ones((dt/12,4), dtype=np.float64)*fillVal
    6393             seastimes=np.zeros((dt/12,4), dtype=np.float64)
    6394             dates = realdatetime_CFcompilant(tvals, dateref, tunits)
    6395 
    6396             for iseas in range(4):
    6397                 for it in range(dt):
    6398                     if tseasons[it,iseas]:
    6399                         tfirstseas=int(it/11)
    6400                         firstseas=True
    6401 ##                        print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
    6402                         break
    6403                 for itt in range(it, dt-1,12):
    6404                     seasvals[int(itt/12),iseas]=(values[itt] + values[itt+1] + values[itt+2])/3.
    6405                     seastimes[int(itt/12),iseas]=dates[itt] + (dates[itt+2] - dates[itt])/2.
    6406 ##                    print itt, values[itt], values[itt+1], values[itt+2], '--->', seasvals[int(itt/12),iseas]
    6407 
    6408             self.dimt = dt/12
    6409             self.seasv = ma.masked_equal(seasvals, fillVal)
    6410             self.seast = ma.masked_equal(seastimes, fillVal)
    6411 
    6412 class cls_seasonal_accums(object):
    6413     """ Class to compute the seasonal accumulations of a time series of values
    6414     tvals=time values (year, month, day, hour, minute, second) format
    6415     values=time-series of values
    6416     dateref= reference date [YYYY][MM][DD][HH][MI][SS] format
    6417     tunits= time units
    6418     self.dimt = Number of values
    6419     self.seasv = seasonal mean values
    6420     self.seast = seasonal time mean values
    6421     """
    6422     def __init__(self, tvals, values, dateref, tunits):
    6423         import numpy.ma as ma
    6424         fillVal=1.e20
    6425 
    6426         if tvals is None:
    6427             self.dimt = None
    6428             self.seasv = None
    6429             self.seast = None
    6430 
    6431         else:
    6432             tseasons=times_4seasons(tvals, True)
    6433 
    6434             dt=len(values)
    6435             seasvals=np.ones((dt/12,4), dtype=np.float64)*fillVal
    6436             seastimes=np.zeros((dt/12,4), dtype=np.float64)
    6437             dates = realdatetime_CFcompilant(tvals, dateref, tunits)
    6438 
    6439             for iseas in range(4):
    6440                 for it in range(dt):
    6441                     if tseasons[it,iseas]:
    6442                         tfirstseas=int(it/11)
    6443                         firstseas=True
    6444 #                        print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
    6445                         break
    6446                 for itt in range(it, dt-1,12):
    6447                     seasvals[int(itt/12),iseas]=values[itt] + values[itt+1] + values[itt+2]
    6448                     seastimes[int(itt/12),iseas]=dates[itt] + (dates[itt+2] - dates[itt])/2.
    6449 
    6450             self.dimt = dt/12
    6451             self.seasv = ma.masked_equal(seasvals, fillVal)
    6452             self.seast = ma.masked_equal(seastimes, fillVal)
    6453 
    6454 def seasonal_means(tvals, values):
    6455     """ Function to compute the seasonal mean of a time series of values
    6456     tvals=time values (year, month, day, hour, minute, second) format
    6457     values=time-series of values
    6458     """
    6459     fillVal=1.e20
    6460 
    6461     tseasons=times_4seasons(tvals, True)
    6462 
    6463     dt=len(values)
    6464     seasvals=np.ones((dt/4,4), dtype=np.float64)
    6465     seasvals=seasvals*fillVal
    6466 
    6467     for iseas in range(4):
    6468         for it in range(dt):
    6469             if tseasons[it,iseas]:
    6470                 tfirstseas=int(it/11)
    6471                 firstseas=True
    6472 #                print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
    6473                 break
    6474         for itt in range(it, dt-1,12):
    6475             seasvals[int(itt/12),iseas]=(values[itt] + values[itt+1] + values[itt+2])/3.
    6476 
    6477     return seasvals
    6478 
    6479 def seasonal_accum(tvals, values):
    6480     """ Function to compute the seasonal accumulation of a time series of values
    6481     tvals=time values (year, month, day, hour, minute, second) format
    6482     values=time-series of values
    6483     """
    6484     fillVal=1.e20
    6485 
    6486     tseasons=times_4seasons(tvals, True)
    6487 
    6488     dt=len(values)
    6489     seasvals=np.ones((dt/4,4), dtype=np.float64)
    6490     seasvals=seasvals*fillVal
    6491 
    6492     for iseas in range(4):
    6493         for it in range(dt):
    6494             if tseasons[it,iseas]:
    6495                 tfirstseas=int(it/11)
    6496                 firstseas=True
    6497 #                print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
    6498                 break
    6499         for itt in range(it, dt-1,12):
    6500             seasvals[int(itt/11),iseas]=values[itt] + values[itt+1] + values[itt+2]
    6501 
    6502     return seasvals
    6503 
    6504 def seasmean(timename, filename, varn):
    6505     """ Function to compute the seasonal mean of a variable
    6506     timename= name of the variable time in the file
    6507     filename= netCDF file name
    6508     varn= name of the variable
    6509     """
    6510     fillValue=1.e20
    6511 
    6512     ncf = NetCDFFile(filename,'a')
    6513     ofile = 'seasmean_' + varn + '.nc'
    6514 
    6515     if not ncf.variables.has_key(varn):
    6516         print errormsg
    6517         print '    seasmean: File  "' + filename + '" does not have variable ' + varn + ' !!!!'
    6518         print errormsg
    6519         ncf.close()
    6520         quit(-1)
    6521 
    6522     varobj = ncf.variables[varn]
    6523     varinf = variable_inf(varobj)
    6524     timeobj = ncf.variables[timename]
    6525     timeinf = cls_time_information(ncf, timename)
    6526 
    6527     timevals=timeobj[:]
    6528     netcdftimes = netCDFdatetime_realdatetime(timeinf.units + ' since ' + timeinf.Urefdate, timeinf.calendar, timevals)
    6529    
    6530     timeseasons=times_4seasons(netcdftimes, True)
    6531 
    6532     timeseasvals=seasonal_means(netcdftimes, timevals)
    6533    
    6534     ncfnew = NetCDFFile(ofile,'w')
    6535     ncfnew.createDimension('seastime', timeinf.dimt/4)
    6536     ncfnew.createDimension('seas', 4)
    6537 
    6538     newvarobj = ncfnew.createVariable('seastime', 'f4', ('seastime', 'seas',), fill_value=fillValue)
    6539     newvarobj[:]=timeseasvals
    6540     attr = newvarobj.setncattr('calendar', timeinf.calendar)
    6541     attr = newvarobj.setncattr('units', timeinf.units)
    6542 
    6543     newvarobj = ncfnew.createVariable('seasn', str, ('seas'))
    6544     attr = newvarobj.setncattr('standard_name', 'seasn')
    6545     attr = newvarobj.setncattr('long_name', 'name of the season')
    6546     attr = newvarobj.setncattr('units', '3-months')
    6547     newvarobj[0]='DJF'
    6548     newvarobj[1]='MAM'
    6549     newvarobj[2]='JJA'
    6550     newvarobj[3]='SON'
    6551 
    6552     newdims=[]
    6553     idim = 0
    6554 
    6555     for ndim in varinf.dimns:
    6556         if ndim == timename:
    6557             newdims.append(unicode('seastime'))
    6558             if not idim == 0:
    6559                 print errormsg
    6560                 print '  seasmean: File  "' + filename + '" does not have time dimension "' + timename + '" as the first one. It has it as # ', idim, ' !!!!'
    6561                 print errormsg
    6562                 ncf.close()
    6563                 quit(-1)
    6564         else:
    6565             newdims.append(ndim)
    6566 
    6567         if not ncfnew.dimensions.has_key(ndim):
    6568             ncfnew.sync()
    6569             ncfnew.close()
    6570             ncf.close()           
    6571 
    6572             fdimadd(filename + ',' + ndim, ofile)
    6573             ncf = NetCDFFile(filename,'r')
    6574             ncfnew = NetCDFFile(ofile,'a')
    6575 
    6576     print ncfnew.dimensions
    6577     namedims=tuple(newdims)
    6578     print namedims
    6579 
    6580     varobj = ncf.variables[varn]
    6581     varinf = variable_inf(varobj)
    6582 
    6583     newvarobj=ncfnew.createVariable(varn + '_seasmean', 'f4', namedims, fill_value=fillValue)
    6584     attr = newvarobj.setncattr('standard_name', varn + '_seasmean')
    6585     attr = newvarobj.setncattr('long_name', 'seasonal mean of ' + varn)
    6586     attr = newvarobj.setncattr('units', varinf.units)
    6587 
    6588     newvar = np.ones(varinf.dims[0], dtype=varinf.dtype)
    6589     if varinf.Ndims == 1:
    6590         newvar = newvarobj[:]
    6591         newvarobj[:] = seasonal_means(netcdftimes, newvar)
    6592     if varinf.Ndims == 2:
    6593         for i in range(varinf.dims[1]):
    6594             newvar = newvarobj[:,i]
    6595             newvarobj[:,i] = seasonal_means(netcdftimes, newvar)
    6596     if varinf.Ndims == 3:
    6597         for i in range(varinf.dims[1]):
    6598             for j in range(varinf.dims[2]):
    6599                 newvar = newvarobj[:,i,j]
    6600                 newvarobj[:,i,j] = seasonal_means(netcdftimes, newvar)
    6601     if varinf.Ndims == 4:
    6602         for i in range(varinf.dims[1]):
    6603             for j in range(varinf.dims[2]):
    6604                 for k in range(varinf.dims[3]):
    6605                     newvar = newvarobj[:,i,j,k]
    6606                     newvarobj[:,i,j,k] = seasonal_means(netcdftimes, newvar)
    6607     if varinf.Ndims == 5:
    6608         for i in range(varinf.dims[1]):
    6609             for j in range(varinf.dims[2]):
    6610                 for k in range(varinf.dims[3]):
    6611                     for l in range(varinf.dims[4]):
    6612                         newvar = newvarobj[:,i,j,k,l]
    6613                         newvarobj[:,i,j,k,l] = seasonal_means(netcdftimes, newvar)
    6614     if varinf.Ndims == 6:
    6615         for i in range(varinf.dims[1]):
    6616             for j in range(varinf.dims[2]):
    6617                 for k in range(varinf.dims[3]):
    6618                     for l in range(varinf.dims[4]):
    6619                         for m in range(varinf.dims[5]):
    6620                             newvar = newvarobj[:,i,j,k,l,m]
    6621                             newvarobj[:,i,j,k,l,m] = seasonal_means(netcdftimes, newvar)
    6622     else:
    6623         print errormsg
    6624         print '  seasmean: number of dimensions ', varinf.Ndims, ' not ready !!!!'
    6625         print errormsg
    6626         ncf.close()
    6627         quit(-1)
    6628 
    6629     ncfnew.sync()
    6630     ncfnew.close()
    6631 
    6632     print 'seasmean: Seasonal mean file "' + ofile + '" written!'
    6633 
    6634 def load_variable_lastdims(var, prevdims, Nlastdims):
    6635     """ Function to load the last [Nlastdims] dimensions of a variable [var] at the other dimensions at [prevdims]
    6636     >>> load_variable_lastdims(np.array(range(5*5*5)).reshape(5,5,5), [1,2], 1)
    6637     [35 36 37 38 39]
    6638     """
    6639     fname='load_variable_lastdims'
     3257    fname='load_ncvariable_lastdims'
     3258
     3259    var=ncf.variables[varn]
    66403260
    66413261    dims=var.shape
    … …  
    66993319    return loadvar
    67003320
    6701 def load_ncvariable_lastdims(ncf, varn, prevdims, Nlastdims):
    6702     """ Function to load the last [Nlastdims] dimensions of a variable [varn] from a netCDF object at the other dimensions at [prevdims]
    6703     ncf= netCDF object
    6704     varn= variable name
    6705     prevdims= values at the previous dimensions
    6706     Nlastims= number of last dimensions
    6707     """
    6708     fname='load_ncvariable_lastdims'
    6709 
    6710     var=ncf.variables[varn]
    6711 
    6712     dims=var.shape
    6713     Ndims=len(dims)
    6714 
    6715     Nprevdims = len(prevdims)
    6716     if not Nprevdims + Nlastdims == Ndims:
    6717         print erromsg
    6718         print '  ' + fname + ': number of dimensions previous (',Nprevdim,') and last (',Nlastdims, ') does not match with variable size: ',Ndims,' !!!!'
    6719         print errormsg
    6720         quit(-1)
    6721 
    6722     if Nlastdims > Ndims-1:
    6723         print errormsg
    6724         print '  ' + fname + ': number of last dimensions ', Nlastdims,' >= than the number of dimensions of the variable ', Ndims,' !!!!'
    6725         print errormsg
    6726         quit(-1)
    6727 
    6728     if Ndims == 1:
    6729         loadvar = var[:]
    6730     elif Ndims == 2:
    6731         loadvar = var[prevdims[0], :]
    6732     elif Ndims == 3:
    6733         if Nlastdims == 1:
    6734             loadvar = var[prevdims[0], prevdims[1], :]
    6735         else:
    6736             loadvar = var[prevdims[0], :, :]
    6737     elif Ndims == 4:
    6738         if Nlastdims == 1:
    6739             loadvar = var[prevdims[0], prevdims[1], prevdims[2], :]
    6740         elif Nlastdims == 2:
    6741             loadvar = var[prevdims[0], prevdims[1], :, :]
    6742         else:
    6743             loadvar = var[prevdims[0], :, :, :]
    6744     elif Ndims == 5:
    6745         if Nlastdims == 1:
    6746             loadvar = var[prevdims[0], prevdims[1], prevdims[2], prevdims[3], :]
    6747         elif Nlastdims == 2:
    6748             loadvar = var[prevdims[0], prevdims[1], prevdims[2], :, :]
    6749         elif Nlastdims == 3:
    6750             loadvar = var[prevdims[0], prevdims[1], :, :, :]
    6751         else:
    6752             loadvar = var[prevdims[0], :, :, :, :]
    6753     elif Ndims == 6:
    6754         if Nlastdims == 1:
    6755             loadvar = var[prevdims[0], prevdims[1], prevdims[2], prevdims[3], prevdims[4], :]
    6756         elif Nlastdims == 2:
    6757             loadvar = var[prevdims[0], prevdims[1], prevdims[2], prevdims[3], :, :]
    6758         elif Nlastdims == 3:
    6759             loadvar = var[prevdims[0], prevdims[1], prevdims[2], :, :, :]
    6760         elif Nlastdims == 4:
    6761             loadvar = var[prevdims[0], prevdims[1], :, :, :, :]
    6762         else:
    6763             loadvar = var[prevdims[0], :, :, :, :]
    6764     else:
    6765         print errormsg
    6766         print '  ' + fname + ' variable size ', Ndims, ' not ready!!!'
    6767         print errormsg
    6768         quit(-1)
    6769 
    6770     return loadvar
     3321AQUI AQUI AQUI AQUI AQUI AQUI AQUI
     3322
    67713323
    67723324def fill_ncvariable_lastdims(ncf, varn, prevdims, Nlastdims, fillvals):
Note: See TracChangeset for help on using the changeset viewer.