source: lmdz_wrf/trunk/tools/model_graphics.py @ 1165

Last change on this file since 1165 was 1160, checked in by lfita, 9 years ago

Seem to be fixed all the issues

File size: 179.1 KB
Line 
1# Python script to generate the files and the plots for a sensitivity study with multiple modle configurations and models
2#   L. Fita
3#   LMD, Jussieu/Palaisseau, France
4#  Script configuration get from ASCII file 'model_graphics.dat'
5
6import numpy as np
7import nc_var_tools as ncvar
8import generic_tools as gen
9import drawing_tools as drw
10import time as tim
11#  To avoid errors like: '/bin/sh: 1: Syntax error: Bad fd number'
12#    Make sure that /bin/sh directs to bash and not to dash:
13#  http://stackoverflow.com/questions/15809060/sh-syntax-error-bad-fd-number
14import subprocess as sub
15import os
16
17main = 'model_graphics.py'
18
19errmsg = ncvar.errormsg
20warnmsg = ncvar.warnmsg
21
22def verify_configuration(ExpConf, debug):
23    """ Function to verify mandatory keys from experiment configuration
24      ExpConf= dictionary with the configuration
25    """
26    fname = 'verify_configuration'
27 
28    # List with the mandatory keys
29    mandatorykeys = ['pyHOME', 'cdoHOME', 'scratch', 'filescratch', 'figscratch',    \
30      'diffscratch', 'figdiffscratch', 'figallmodexpscratch', 'addfiles',            \
31      'addfigures', 'adddiffs', 'adddifffigures', 'addallmodexpfigures', 'debug',    \
32      'models', 'modgraphchar', 'expgraphchar', 'ifold', 'ofold', 'warnmsg',         \
33      'errmsg', 'titleoperations', 'kindfig', 'CFreftime', 'CFunitstime', 'mapval',  \
34      'timekind', 'timefmt', 'timelabel' ]
35
36    # List with the optional keys
37    optionalkeys = ['specificvarplot', 'specificdiffopplot', 'specificdiffvarplot',  \
38      'RefProj']
39
40    # Dictionary with the mandatory keys which are required as function of a given
41    #   key from the dictionary
42    #     if [key] must exist [keyB] = [values]
43    ifkeys = {'WRFexps': ['models','WRF'], 'LMDZexps': ['models','LMDZ'],            \
44      'WRF_LMDZexps': ['models','WRF_LMDZ'], 'WRFheaders': ['models','WRF'],         \
45      'WRF_LMDZheaders': ['models','WRF_LMDZ'], 'LMDZheaders': ['models','LMDZ']}
46
47    # Dictionary with the optional keys (at least one) which are required as function
48    #   of a given key from the dictionary
49    #     if [key] must exist at least one of the [keyA, keyB, keyC, ...] with values
50    ifopkeys = {'RefProj': ['CDOremapvar_dis bil', 'CDOremapvar_bic',                \
51      'CDOremapvar_dis', 'CDOremapvar_nn', 'CDOremapvar_con', 'CDOremapvar_con2',    \
52      'CDOremapvar_laf']}
53
54    for kdict in mandatorykeys:
55        if not ExpConf.has_key(kdict):
56            print errmsg
57            print '  ' + fname + "' configuration without required '" + kdict + "' !!"
58            quit(-1)
59        if debug:
60            print ExpConf[kdict]
61
62    for kdict in ifkeys.keys():
63        vals = ifkeys[kdict]
64        keyn = vals[0]
65        keyv = vals[1]
66        if type(ExpConf[keyn]) == type('A'):
67            if not ExpConf.has_key(kdict) and ExpConf[keyn] == keyv:
68                print errmsg
69                print '  ' + fname + "' configuration without '" + kdict +           \
70                  "' when it is required with configuration '" + keyn + '=' + keyv + \
71                  "' !!"
72                quit(-1)
73        elif type(ExpConf[keyn]) == type(['A', 'B']):
74            keyvals = ExpConf[keyn]
75            if not ExpConf.has_key(kdict) and gen.searchInlist(keyvals, keyv):
76                print errmsg
77                print '  ' + fname + "' configuration without '" + kdict +           \
78                  "' when it is required with configuration '" + keyn + '= [..., ' + \
79                 keyv + ", ...]' !!"
80                quit(-1)
81        else:
82            print errmsg
83            print '  ' +fname+ 'dictionary type ', type(ExpConf[keyn]), ' not ready!!'
84            quit(-1)
85        if debug:
86            print ExpConf[kdict]
87
88    for kdict in ifopkeys.keys():
89        if ExpConf.has_key(kdict):
90            opkeys = ifopkeys[kdict]
91            hasopkey = False
92            opdictvals = {}
93            for opkv in opkeys:
94                if ExpConf.has_key(opkv):
95                    hasopkey = True
96                    opdictvals[opkv] = ExpConf[opkv]
97
98            if not hasopkey:
99                print errmsg
100                print '  ' + fname + "' configuration wit '" + kdict + "' but " +    \
101                  "without any key value:", opkeys, "' when it is required !!"
102                quit(-1)
103
104            if debug:
105                print "  Optional key '" + kdict + "' with values ________"
106                gen.printing_dictionary(opdictvals)
107
108    print fname + ': configuration seems to be fine'
109
110    return
111
112def scratches(config):
113    """ Function to set-up if it is needed to start from the scratch
114      config: dictionary with the configuration
115    """
116    fname = 'scratches'
117
118    if config['scratch'] == 'true':
119        scr = True
120        print warnmsg
121        print "  " + fname + ": starting from the SCRATCH !!"
122        print "    10 seconds left!!"
123        filescr = True
124        figscr = True
125        difscr = True
126        figdifscr = True
127        moddifscr = True
128        figmoddifscr = True
129        figallmodexpscr = True
130        tim.sleep(10)
131    else:
132        scr = False
133        if config['filescratch'] == 'true':
134            filescr = True
135            print warnmsg
136            print "  " + main + ": files starting from the SCRATCH !!"
137            print "    5 seconds left!!"
138            tim.sleep(5)
139        else:
140            filescr = False
141
142        if config['figscratch'] == 'true':
143            figscr = True
144            print warnmsg
145            print "  " + main + ": figures starting from the SCRATCH !!"
146            print "    5 seconds left!!"
147            tim.sleep(5)
148        else:
149            figscr = False
150
151        if config['diffscratch'] == 'true':
152            difscr = True
153            print warnmsg
154            print "  " + main + ": experiment differences starting from the SCRATCH !!"
155            print "    5 seconds left!!"
156            tim.sleep(5)
157        else:
158            difscr = False
159
160        if config['figdiffscratch'] == 'true':
161            figdifscr = True
162            print warnmsg
163            print "  " + main + ": figures experiment differences starting from the SCRATCH !!"
164            print "    5 seconds left!!"
165            tim.sleep(5)
166        else:
167            figdifscr = False
168
169        if config['moddiffscratch'] == 'true':
170            moddifscr = True
171            print warnmsg
172            print "  " + main + ": model differences starting from the SCRATCH !!"
173            print "    5 seconds left!!"
174            tim.sleep(5)
175        else:
176            moddifscr = False
177
178        if config['figmoddiffscratch'] == 'true':
179            figmoddifscr = True
180            print warnmsg
181            print "  " + main + ": figures model differences starting from the SCRATCH !!"
182            print "    5 seconds left!!"
183            tim.sleep(5)
184        else:
185            figmoddifscr = False
186
187        if config['figallmodexpscratch'] == 'true':
188            figallmodexpscr = True
189            print warnmsg
190            print "  " + main + ": figures all model-experiment starting from the SCRATCH !!"
191            print "    5 seconds left!!"
192            tim.sleep(5)
193        else:
194            figallmodexpscr = False
195
196    if config['addall'] == 'true':
197        print warnmsg
198        print "  " + fname + ": adding new values everywhere !!"
199        addfils = True
200        addfigs = True
201        adddiffs = True
202        adddifffigs = True
203        addmoddiffs = True
204        addmoddifffigs = True
205        addallmodexpfigs = True
206    else:
207        if config['addfiles'] == 'true':
208            addfils = True
209        else:
210            addfils = False
211
212        if config['addfigures'] == 'true':
213            addfigs = True
214        else:
215            addfigs = False
216
217        if config['adddiffs'] == 'true':
218            adddiffs = True
219        else:
220            adddiffs = False
221
222        if config['adddifffigures'] == 'true':
223            adddifffigs = True
224        else:
225            adddifffigs = False
226
227        if config['addmoddiffs'] == 'true':
228            addmoddiffs = True
229        else:
230            addmoddiffs = False
231
232        if config['addmoddifffigures'] == 'true':
233            addmoddifffigs = True
234        else:
235            addmoddifffigs = False
236
237        if config['addallmodexpfigures'] == 'true':
238            addallmodexpfigs = True
239        else:
240            addallmodexpfigs = False
241
242    if config['debug'] == 'true':
243        debug = True
244    else:
245        debug = False
246
247    return scr, filescr, figscr, difscr, figdifscr, moddifscr, figmoddifscr,         \
248      figallmodexpscr, addfils, addfigs, adddiffs, adddifffigs, addmoddiffs,         \
249      addmoddifffigs, addallmodexpfigs, debug
250
251
252class gc(object):
253    """ Class which holds all the graphical characteristics
254      [label]: label as it should appear in the graphics
255      [color]: specific color of the model
256      [linetype]: specific type of line of the model
257      [marker]: specific marker of the model
258      [sizes]: line width and point size for the model
259      [tmod]: time modification to apply to the model files (None for nothing)
260        'setorigin',[YYYYMMDDHHMISS]: re-set origin of times at [YYYYMMDDHHMISS]
261    """
262    def __init__( self, label, color, linetype, marker, sizes, Tmod):
263        self.label = None
264        self.color = None
265        self.linetype = None
266        self.marker = None
267        self.sizes = None
268        self.tmod = None
269        if label is not None:
270            self.label = label
271            self.color = color
272            self.linetype = linetype
273            self.marker = marker
274            self.sizes = sizes
275            if Tmod != 'None':
276                self.tmod = Tmod
277
278def graphical_characteristics(config, debug):
279    """ Function to provide the graphical characteristics of models and experiments
280      config= configuration of the experiment
281    """
282    fname = 'graphical_characteristics'
283
284    modgraphchar = {}
285    modvals = config['modgraphchar'].split(':')
286    for modval in modvals:
287        modv = modval.split('|')
288        modgraphchar[modv[0]] = gc(modv[1], modv[2], modv[3], modv[4], modv[5], modv[6])
289
290    expgraphchar = {}
291    expvals = config['expgraphchar'].split(':')
292    for expval in expvals:
293        expvs = expval.split('|')
294        modexp = expvs[0] + '/' + expvs[1]
295        if not modgraphchar.has_key(expvs[0]):
296            print errmsg
297            print '  ' + fname + ": a model called '" +expvs[0]+ "' does not exist !!"
298            print '    existing models with graphic characteristics:',               \
299              modgraphchar.keys()
300            quit(-1)
301
302        modvs = modgraphchar[expvs[0]]
303        expv = expvs[2:7]
304        if expv[1] == 'asmodel': expv[1] = modvs.color
305        if expv[2] == 'asmodel': expv[2] = modvs.linetype
306        if expv[3] == 'asmodel': expv[3] = modvs.marker
307        if expv[4] == 'asmodel': expv[4] = modvs.sizes
308
309        expgraphchar[modexp] = gc(expv[0], expv[1], expv[2], expv[3], expv[4], 'None')
310
311    # Running some tests
312    mods = config['models'].split(':')
313    for mod in mods:
314        if not modgraphchar.has_key(mod):
315            print errmsg
316            print '  ' + fname + ": model called '" + mod + "' does not have " +     \
317              "graphical characteristics !!"
318            print '    existing models with graphic characteristics:',               \
319              modgraphchar.keys()
320            quit(-1)
321
322        exps = config[mod + 'exps'].split(':')
323        for exp in exps:
324            modexp = mod + '/' + exp
325            if not expgraphchar.has_key(modexp):
326                print errmsg
327                print '  ' + fname + ": model/experiment called '" + modexp +        \
328                  "' does not have graphical characteristics !!"
329                print '    existing model/experiments with graphic ' +               \
330                  'characteristics:', expgraphchar.keys()
331                quit(-1)
332
333    if debug:
334        print '  ' + fname + ': model graphical characteristics _______'
335        for modn in modgraphchar.keys():
336            with gen.Capturing() as output:
337                gen.printing_class(modgraphchar[modn])
338            print "'" + modn + "'; " + ', '.join(output)
339        print '  ' + fname + ': experiment graphical characteristics _______'
340        for modexpn in expgraphchar.keys():
341            with gen.Capturing() as output:
342                gen.printing_class(expgraphchar[modexpn])
343            print "'" + modexpn + "'; " + ', '.join(output)
344
345    return modgraphchar, expgraphchar
346
347def exp_headers(mod,config):
348    """ Function to provide the headers and the experiments of a given model
349      model= model
350      config= configuration of the experiment
351    """
352    fname = 'exp_headers'
353
354    # No CASE in python!
355    #  case ${mod} in ...
356    if mod == 'WRF':
357        expers = config['WRFexps'].split(':')
358        fhs = config['WRFheaders'].split(':')
359    elif mod == 'LMDZ':
360        expers = config['LMDZexps'].split(':')
361        fhs = config['LMDZheaders'].split(':')
362    elif mod == 'WRF_LMDZ':
363        expers = config['WRF_LMDZexps'].split(':')
364        fhs = config['WRF_LMDZheaders'].split(':')
365    else:
366        print errmsg
367        print "  " + fname + ": model '" + mod + "' not ready!!"
368        quit(-1)
369
370    return expers, fhs
371
372class VariableInf(object):
373    """ Class which holds all information of a given variable
374      name= CF name of the variable
375      header= header of the file from which it can be computed
376      model= file-header variable from which it can be directly computed
377      diag= file-header diagnostics from which it can be computed
378    """
379    def __init__( self, name, fheader, model, diag):
380       self.name= None
381       self.fheader= None
382       self.model= None
383       self.diag= None
384       if name is not None:
385           self.name = name
386           self.fheader= fheader
387           self.model= model
388           self.diag= diag
389
390class ModelInf(object):
391    """ Class which holds all information from a given model
392      name= Acronym of the model
393      model= Long description of the model
394      dim[x/y/z/t/s]n= name of the dimensions of the model
395      vdim[x/y/z/t/s]n= name of the variable-dimensions of the model
396      dimensions= list of the dimensions
397      vardimensions= list of the variable-dimensions
398      timemodif= time modification to apply to the model files ('None' for nothing)
399        'setorigin',[YYYYMMDDHHMISS]: re-set origin of times at [YYYYMMDDHHMISS]
400    """
401    def __init__( self, name, model, dx, dy, dz, dt, ds, vdx, vdy, vdz, vdt, vds,    \
402      tmod):
403        self.name= None
404        self.model= None
405        self.dimxn= None
406        self.dimyn= None
407        self.dimzn= None
408        self.dimtn= None
409        self.dimsn= None
410        self.vardxn= None
411        self.vardyn= None
412        self.vardzn= None
413        self.vardtn= None
414        self.vardsn= None
415        self.timemodif = None
416        if name is not None:
417            self.name = name
418            self.model= model
419            self.dimxn= dx
420            self.dimyn= dy
421            self.dimzn= dz
422            self.dimtn= dt
423            self.dimsn= ds
424            self.vardxn= vdx
425            self.vardyn= vdy
426            self.vardzn= vdz
427            self.vardtn= vdt
428            self.vardsn= vds
429            if tmod is not None:
430                self.timemodif = tmod
431            self.dimensions = [dx, dy, dz, dt, ds]
432            self.vardimensions = [vdx, vdy, vdz, vdt, vds]
433
434def variable_compute(idir,var,ftests,db):
435    """ Function to retrieve the computation way of a given variable using a series of test files
436      iwdir= directory with the test files
437      var= name of the variable to test
438      filtests= dictionary with the test files for each header
439    """
440    fname='variable_compute'
441
442    cancompute = True
443    for headerf in ftests.keys():
444        filen=ftests[headerf]
445        try:
446            with gen.Capturing() as output:
447                vmod, vdiag = ncvar.computevar_model(var, idir + '/' + filen)
448        except:
449            print errmsg
450            print 'ncvar.computevar_model(' + var + ', ' + idir + '/' + filen + ')'
451            for s1out in output: print s1out
452            quit(-1)
453
454        if vmod is None and vdiag is None:
455            cancompute = False
456        else:
457            cancompute = True
458            # Should be considered that variable can also be computed by both ways?
459            break
460
461    if not cancompute:
462        print warnmsg
463        print '  ' + fname + ": there is no way to compute '" + var +                \
464          "' for model '" + mod
465# Too extrict!
466#        quit(-1)
467
468    # Getting only the first possibility of model and diagnostic
469    if vmod is not None:
470        modV = vmod[0]
471    else:
472        modV = None
473    if vdiag is not None:
474        diagV = vdiag[0]
475    else:
476        diagV = None
477
478    varcomp = VariableInf(var, headerf, modV, diagV)
479
480    # a ';' list 'varcompute' it is created for each variable giving:
481    #   [var]|[vark]|[headerf][varmod]|[vardiag]
482    # This list will be used to compute a new file for each variable
483
484    if db:
485        print 'Variable information _______'
486        gen.printing_class(varcomp)
487
488    return varcomp
489
490def get_operations_var(vc,db):
491    """ Function to provide the operations to make for each variable from 'VAR_' dictionary
492      vc= 'VAR_' dictionary, dictionary with all the parameters started with 'VAR_'
493    """
494    fname = 'get_operations_var'
495
496    LH = len('VAR_')
497    iop = 1
498    # list of operations
499    doopers = {}
500    # operations by variable
501    operationsvar = {}
502    # individual operations by variable (not repeated)
503    indivoperationsvar = {}
504
505    # Variables with a calculation whic requires vertical interpolation to all the
506    #   file (operations stating by 'VAR_pinterp')
507    LVp = len('VAR_pinterp')
508    varglobalp = []
509
510    for oper in vc.keys():
511        Loper = len(oper)
512        opn = oper[LH:Loper+1]
513        vns = vc[oper].split(':')
514        ops1 = opn.split('+')
515        doopers[opn] = ops1
516        for vn in vns:
517            if not operationsvar.has_key(vn):
518                operationsvar[vn] = [opn]
519                indivoperationsvar[vn] = ops1
520            else:
521                opers = operationsvar[vn]
522                opers.append(opn)
523                operationsvar[vn] = opers
524                opers1 = indivoperationsvar[vn]
525                for op1 in ops1:
526                    if not gen.searchInlist(opers1, op1):
527                        opers1.append(op1)
528                indivoperationsvar[vn] = opers1
529
530        if oper[0:LVp] == 'VAR_pinterp':
531            for vn in vns:
532                if not gen.searchInlist(varglobalp,vn): varglobalp.append(vn)
533
534        iop = iop + 1
535
536    if db:
537        print '  operations to make _______'
538        gen.printing_dictionary(doopers)
539        print '  operations by variable _______'
540        gen.printing_dictionary(operationsvar)
541        print '  individual operations by variable _______'
542        gen.printing_dictionary(indivoperationsvar)
543        print '  variables with the vertical interpolation for all the file _______'
544        print '    #', varglobalp
545
546    return doopers, operationsvar, indivoperationsvar, varglobalp
547
548def pinterp_var(oper):
549    """ Function to retrieve characteristics of the vertical interpolation for the operation
550      oper= operation
551        Wether vertical interpolation is:
552          # 'global': for all file
553          # 'local': at a given step of the process
554          # 'none': no vertical interpolation
555    >>> pinterp_var('last+pinterp+xmean')
556    local
557    >>> pinterp_var('pinterp+tmean+xmean')
558    global
559    >>> pinterp_var('xmean')
560    none
561    """
562    fname = 'pinterp_var'
563
564    pinn = 'pinterp'
565
566    Lpinterp = len(pinn)
567    if oper[0:Lpinterp] == pinn: 
568        pkind = 'global'
569        return pkind
570
571    if oper.find(pinn) != -1:
572        pkind = 'local'
573    else:
574        pkind = 'none'
575       
576    return pkind
577
578def compvars_listconstruct(config, minf, Files, TestFiles, idir, odir, debug):
579    """ Function to construct the list of variables to compute
580      config= dictionary with the configuration of the execution
581         variables to compute are get with all values started by `VAR_' in 'module_graphics.dat', and then
582           providing a consecutive number of calculations separated by '+'
583             VAR_[calc1]+[calc2] = tas:wss
584           will compute first [calc1] and then [calc2] for 'tas' and 'wss'
585      minf= class with information about the model
586      Files= dictionary of files for heach header
587      TestFiles= dictionary of files to test how to compute the variable
588      odir= output directory
589    """
590    fname='compvars_listconstruct'
591
592    varcomp = gen.get_specdictionary_HMT(config,H='VAR_')
593
594    if debug:
595        print '  variables to compute ________'
596        gen.printing_dictionary(varcomp)
597
598    # Getting operations by variable
599    opers, varoper, indivaroper, varglobp = get_operations_var(varcomp,debug)   
600
601    strfiles = gen.dictKeysVals_stringList(Files)
602    strtestfiles = gen.dictKeysVals_stringList(TestFiles)
603
604    Svarcompute = ''
605
606    # Main dictionary with all the variables and their calculation
607    allvarcomp = {}
608
609    ivop = 0
610    for vn in varoper:
611        vcomp = variable_compute(idir,vn,TestFiles,False)
612        for op in varoper[vn]:
613            # Creation of a String as: [CFname]|[operation]|[fheader]|[vmodel]|[vdiag]|[globalP]
614            #   [CFname]: CF-name of the variable (must appear in 'variables_values.dat')
615            #   [operation]: operation to compute
616            #   [fheader]: header of file with the required variables
617            #   [vmodel]: Variable from model output
618            #   [vdiag]: Varible as diagnostic from different variables of the model output
619            #   [globalP]: Wether vertical interpolation is:
620            #     'global': for all file
621            #     'local': at a given step of the process
622            #     'none': no vertical interpolation
623            globalP = pinterp_var(op)
624            if vcomp.model is not None:
625                if type(vcomp.model) == type(list([1, 2])):
626                    Smodel = ':'.join(vcomp.model)
627                elif type(vcomp.model) == type('A'):
628                    Smodel = vcomp.model
629            else:
630                Smodel = 'None'
631            if vcomp.diag is not None:
632                if type(vcomp.diag) == type(list([1, 2])):
633                    Sdiag = ':'.join(vcomp.diag)
634                elif type(vcomp.diag) == type('A'):
635                    Sdiag = vcomp.diag
636            else:
637                Sdiag = 'None'
638
639            # To get a combination of operations, all the precedent steps are kept
640            #   thus, they are also available !
641            seqops = op.split('+')
642            for seqop in seqops:
643                if seqop == seqops[0]: seqopS=seqop
644                else: seqopS = seqopS + '+' + seqop
645                allvarcomp[vn+'_'+seqopS] = [vcomp.fheader, vcomp.model, vcomp.diag, \
646                  globalP]
647
648                Svc = vcomp.name + '|' + seqopS + '|' + vcomp.fheader + '|' +        \
649                  Smodel + '|' + Sdiag + '|' + globalP
650                if ivop == 0:
651                    Svarcompute = Svc
652                else:
653                    Svarcompute = Svarcompute + ',' + Svc
654
655                ivop = ivop + 1
656
657    if debug:
658        print '  Variables to compute _______'
659        gen.printing_dictionary(allvarcomp)
660
661    # Outwritting the varcompute to avoid next time (if it is not filescratch!)
662    objf = open(odir + '/varcompute.inf', 'w')
663    objf.write('files: ' + strfiles + '\n')
664    objf.write('testfiles: ' + strtestfiles + '\n')
665    objf.write('varcompute: ' + Svarcompute + '\n')
666    objf.write('itotv: ' + str(ivop) + '\n')
667    objf.close() 
668
669    return allvarcomp, ivop
670
671def read_varcomp_file(filen):
672    """ Function to read 'varcompute.inf' and reconstruct the dictionaries
673    """
674    fname = 'read_varcomp_file'
675
676    allvarcomp = {}
677   
678    objf = open(filen, 'r')
679
680    for line in objf:
681        vals = line.replace('\n','').split(' ')
682        if vals[0] == 'files:':
683            Files = gen.stringList_dictKeysVals(vals[1])
684        elif vals[0] == 'testfiles:':
685            TestFiles = gen.stringList_dictKeysVals(vals[1])
686        elif vals[0] == 'varcompute:':
687            for VvalsS in vals[1].split(','):
688               
689                Vvals = VvalsS.split('|')
690                if Vvals[3] == 'None':
691                    mod = None
692                else:
693                    mod = Vvals[3].split(':')
694                if Vvals[4] == 'None':
695                    diag = None
696                else:
697                    diag = Vvals[4].split(':')
698
699                allvarcomp[Vvals[0]+'_'+Vvals[1]] = [Vvals[2], mod, diag, Vvals[5]]
700
701#                # To get a combination of operations, all the precedent steps are kept
702#                #   thus, they are also available !
703#                seqops = Vvals[1].split('+')
704#                for seqop in seqops:
705#                    if seqop == seqops[0]: seqopS=seqop
706#                    else: seqopS = seqopS + '+' + seqop
707#                    allvarcomp[Vvals[0]+'_'+seqopS] = [Vvals[2], mod, diag, Vvals[5]]
708        elif vals[0] == 'itotv:':
709            ivop = int(vals[1])
710
711    objf.close()
712
713    return Files, TestFiles, allvarcomp, ivop
714
715def pinterpS(gP):
716    """ For that variables which require vertical interpolation 'p' suffix to the file header is added
717      gP= kind of pinterp: 'global', 'local', 'none'
718    """
719    fname = 'pinterpS'
720
721    if gP != 'none':
722        gPS = 'p'
723    else:
724        gPS = ''
725
726    return gPS
727
728def compute_variable(minf, idir, usefiles, odir, cvar, gP, scr, pyH, Tref, Tunits, db):
729    """ Function to compute a variable
730      minf= class with the information of the model
731      idir= directory with the input files
732      usefiles= dictionary of files as dict([headerf]) = [file1], ..., [fileN]
733      odir= directory to write the output files
734      cvar= class with the information of the variable: 'name', 'fheader', 'varmod', 'vardiag'
735      gP= kind of vertical interpolation ('global', 'local', 'none')
736      scr= should it be done from the scratch?
737      pyH= location of the python HOME
738      Tref= CF time reference
739      Tunits= CF time units
740    """
741    fname='compute_variable'
742
743    CFvarn=cvar.name
744    headerf = cvar.fheader
745    modvar = cvar.model
746    diagvar = cvar.diag
747
748    cfiles = usefiles[headerf]
749
750    # dimensions
751    dnx = minf.dimxn
752    dny = minf.dimyn
753    # var-dimensions
754    vdnx = minf.vardxn
755    vdny = minf.vardyn
756
757    # Computing separately and then joinging for all files
758    Ntotfiles = len(cfiles)
759    Nzeros = len(str(Ntotfiles))
760    NStot = str(Ntotfiles).zfill(Nzeros)
761
762    # For that variables which require vertical interpolation 'p' suffix to the
763    #   file header is added
764    SgP = pinterpS(gP)
765
766    # Getting in working dir
767    os.chdir(odir)
768
769    # File to keep track of all operations
770    otrackf = open( odir + '/all_computevars.inf', 'a')
771
772    # Computing variable
773    ifile=1
774    for cf in cfiles:
775        ifS = str(ifile).zfill(Nzeros)
776        ifilen = odir + '/' + CFvarn + '_' + headerf + SgP + '_' + ifS + '-' +       \
777          NStot + '.nc'
778        fileon = odir + '/' + CFvarn + '_' + headerf + SgP + '.nc'
779
780        if scr:
781            sout = sub.call('rm ' + ifilen + ' >& /dev/null', shell=True)
782            sout = sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
783
784        if not os.path.isfile(ifilen) and not os.path.isfile(fileon):
785# Since model direct values are retrieved from `variables_valules.dat' which was
786#   initially coincived as a way to only give variable attributes, range and color
787#   bars, if a variable has a diagnostic way to be computed, the later one will be
788#   preferred
789            if db:
790                print '  ' + fname + ": creation of variable file '" + CFvarn +      \
791                  "' in file '" + ifilen + "' ..."
792                print '    model variable:', modvar
793                print '    diagostic variable:', diagvar
794
795            if diagvar is None:
796                # model variable
797                if db: print '  '+fname+ ": variable direct from model '" +modvar+ "'"
798                values = modvar + ',0,-1,-1'
799                vs = modvar + ',' + vdnx + ',' + vdny + ',' + vdnz + ',' + vdnt
800
801                try:
802                    with gen.Capturing() as output:
803                        ncvar.DataSetSection_multivars(values, idir + '/' + cf, vs)
804                except:
805                    print errmsg
806                    print 'ncvar.DataSetSection_multivars('+values+', '+idir+'/'+cf+ \
807                      ', '+vs+')'
808                    for s1out in output: print s1out
809                    quit(-1)
810
811                newfile, loc = gen.search_sec_list(output,'succesfull')
812                ofile = newfile[0].split(' ')[7]
813                sub.call('mv ' + ofile + ' ' + ifilen, shell=True)
814
815                # Keeping track of the operations
816                pyins = pyH + "/nc_var.py -f " + idir + '/' + cf +                   \
817                  " -o DataSetSection_multivars -v " + vs + " -S '" + values + "'"
818                otrackf.write('\n')
819                otrackf.write('# ' + CFvarn + " " + modvar + '\n')
820                otrackf.write('python ' + pyins + '\n')
821
822                # CF renaming of variable
823                ncvar.chvarname(CFvarn,ifilen,modvar)
824            else:
825                # diagnostic variable
826                if db: print '  '+fname+ ": variable as diagnostic '" +diagvar[0]+ "'"
827                dims = dnt+'@'+vdnt+','+dnz+'@'+vdnz+','+dny+'@'+vdny+','+dnx+'@'+vdnx
828                diagn = diagvar[0]
829                Ndiagvars = len(diagvar)
830                diagc = '@'.join(diagvar[1:])
831           
832                values = '-f ' + idir + '/' + cf + " -d '" + dims + "' -v '" +       \
833                  diagn + '|' + diagc + "'"
834                try:
835                    with gen.Capturing() as output:
836                        sout = sub.call('python ' +pyH+ '/diagnostics.py ' + values, \
837                          shell=True)
838                except:
839                    print errmsg
840                    print 'python ' + pyH + '/diagnostics.py ' + values
841                    for s1out in output: print s1out
842                    quit(-1)
843                if db:
844                    for s1out in output: print s1out
845
846                sout = sub.call(' mv diagnostics.nc ' + ifilen, shell=True)
847
848                # Keeping track of the operations
849                pyins = 'python ' + pyH + '/diagnostics.py ' + values
850                otrackf.write('\n')
851                otrackf.write('# ' + CFvarn + " " + diagn + '\n')
852                otrackf.write(pyins + '\n')
853
854            # Attaching necessary variables for the pressure interpolation
855            if gP != 'none':
856                if minf.name == 'WRF' or minf.name == 'WRF_LMDZ':
857                    requiredinterpvars = ['P', 'PB', 'PSFC', 'PH', 'PHB', 'HGT', 'T',\
858                      'QVAPOR', 'XLONG', 'XLAT', 'Times']
859                elif minf.name == 'LMDZ':
860                    requiredinterpvars = ['pres', 'psol', 'geop', 'phis', 'temp',    \
861                      'ovap', 'lon', 'lat', 'time_counter']
862                else:
863                    print errmsg
864                    print fname + ": for model '" + minf.name + "' required " +      \
865                      " variables for vertical interpolation are not known !!"
866                    quit(-1)
867
868                print "   " + fname + ": adding variables:", requiredinterpvars,     \
869                  ' to allow pressure interpolation'
870                for rqv in requiredinterpvars:
871                    try:
872                        with gen.Capturing() as output:
873                            ncvar.fvaradd(idir+'/'+cf+','+rqv,ifilen)
874                    except:
875                        print errmsg
876                        print 'fvaradd('+idir+'/'+cf+','+rqv+', '+ifilen+')'
877                        for s1out in output: print s1out
878                        quit(-1)
879                    if db:
880                        for s1out in output: print s1out
881
882            # CFification of files
883            if minf.name == 'WRF' or minf.name == 'WRF_LMDZ':
884                values = vdnx + ':' + vdny + ':'+ Tref + ':' + Tunits
885                try:
886                    with gen.Capturing() as output:
887                        ncvar.WRF_toCF(values, ifilen) 
888                except:
889                    print errmsg
890                    print 'WRF_toCF('+ values +', ' + ifilen + ')'
891                    for s1out in output: print s1out
892                    quit(-1)
893                if db:
894                    for s1out in output: print s1out
895            elif minf.name == 'LMDZ':
896                try:
897                    with gen.Capturing() as output:
898                        ncvar.LMDZ_toCF(Tunits + '!since!' + Tref, ifilen)
899                except:
900                    print errmsg
901                    print 'LMDZ_toCF('+ Tunits +'!since!'+ Tref + ', ' + ifilen + ')'
902                    for s1out in output: print s1out
903                    quit(-1)
904                if db:
905                    for s1out in output: print s1out
906            else:
907                print errmsg
908                print '  ' + fname + ": no CFification for model '" +minf.name+ "' !!"
909                print "    available ones: 'WRF', 'WRF_LMDZ', 'LMDZ'"
910                quit(-1)
911
912        ifile = ifile + 1
913
914    otrackf.close()
915
916    # Joining variable files
917    if not os.path.isfile(fileon):
918        if db:
919            print '  ' + fname + ": concatenating all variable files to create '" +  \
920              fileon + "..."
921        try:
922            with gen.Capturing() as output:
923                ncvar.netcdf_fold_concatenation_HMT('./,time', CFvarn+'_'+headerf+   \
924                  SgP+'_,-,.nc', 'all')
925        except:
926            print errmsg
927            print 'netcdf_fold_concatenation_HMT(./time, '+CFvarn+'_'+headerf+SgP+   \
928              '_,-,.nc, all)'
929            for s1out in output: print s1out
930            quit(-1)
931        if db:
932            for s1out in output: print s1out
933
934        sout = sub.call('mv netcdf_fold_concatenated_HMT.nc ' + fileon, shell=True)
935        if os.path.isfile(fileon):
936            sout = sub.call('rm '+CFvarn +'_'+ headerf+SgP + '_*-*.nc >& /dev/null',\
937              shell=True)
938
939    # Re-setting time of final concatenated file
940    if minf.timemodif is not None:
941        if db:
942            print '  ' + fname + ": Modifying times of the file by '" +              \
943              minf.timemodif + "' !!"
944        try:
945            with gen.Capturing() as output:
946                ncvar.time_reset(minf.timemodif, fileon,'time')
947        except:
948            print errmsg
949            print 'time_reset(' + minf.timemodif + ', ' + fileon + ', time)'
950            for s1out in output: print s1out
951            quit(-1)
952        if db:
953            for s1out in output: print s1out
954
955    return
956
957def compute_statistics(minf, config, idir, usefiles, odir, cvar, gP, Opers, scr, db):
958    """ Function to compute different statistics it will take previous steps if they
959        are availale
960      minf= class with the information of the model
961      config= dictionary with the configuration of the experiment
962      idir= directory with the input files
963      usefiles= ',' list of files to use [file1],...,[fileN]
964      odir= directory to write the output files
965      cvar= class with the information of the variable: 'name', 'fheader', 'varmod', 'vardiag'
966      gP= kind of vertical interpolation ('global', 'local', 'none')
967      Opers= kind of operation: (as possible multiple consecutive combination of operations separated by '+'
968        [calc1]+[calc2] will compute first [calc1] and then [calc2]
969        acc: temporal accumulated values
970        diff: differences between models
971        direct: no statistics
972        last: last temporal value
973        Lmean: latitudinal mean values
974        Lsec: latitudinal section (latitudinal value must be given, [var]@[lat])
975        lmean: longitudinal mean values
976        lsec: longitudinal section (longitudinal value must be given, [var]@[lat])
977        pinterp: pressure interpolation (to the given $plevels)
978        tmean: temporal mean values
979        tstd: temporal standard deviation values
980        tturb: Taylor's turbulence decomposition value (x - <x>) for time
981        tvar: temporal variance values
982        xmean: x-axis mean values
983        xvar: x-axis variance values
984        ymean: y-axis mean values
985        zsum: vertical aggregated values
986      scr= should it be done from the scratch?
987    """
988    fname='compute_statistics'
989
990    # Getting a previous file name to continue with(for combinations)
991    CFvarn=cvar.name
992    headerf = cvar.fheader
993    modvar = cvar.model
994    diagvar = cvar.diag
995
996    cfiles = usefiles[headerf]
997
998    # dimensions
999    dnx = minf.dimxn
1000    dny = minf.dimyn
1001    # var-dimensions
1002    vdnx = minf.vardxn
1003    vdny = minf.vardyn
1004    # Model dimension-variable dictionary
1005    Modeldimvardict = {dnx: vdnx, dny: vdny, dnz: vdnz, dnt: vdnt}
1006
1007    # Some experiment configuration values
1008    #  plevels= ':' separated list of pressures (in Pa) to use for the vertical interpolation
1009    plevels = config['plevels']
1010    #  pyH= location of the python HOME
1011    pyH = config['pyHOME']
1012    #  Tref= CF time reference
1013    Tref = config['CFreftime']
1014    #  Tunits= CF time units
1015    Tunits = config['CFunitstime']
1016    #  opsur = operation surnames
1017    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
1018    Opsurs = {}
1019    for opk in opsur:
1020        opn = opk.split('_')[1]
1021        vals = opsur[opk]
1022        Opsurs[opn] = vals.split(':')
1023#    if db:
1024#        print '  ' + fname + ' operation surnames _______'
1025#        gen.printing_dictionary(Opsurs)
1026
1027    # Getting in working dir
1028    os.chdir(odir)
1029
1030    # File to keep track of all operations
1031    otrackf = open( odir + '/all_statsvars.inf', 'a')
1032
1033    # For that variables which require vertical interpolation 'p' suffix to the file
1034    #   header is added.
1035    # After some operations some CF dimension-variables might not be kept
1036
1037    if gP != 'none':
1038        SgP = 'p'
1039        varnCFs = ['lon', 'lat', 'pres', 'time']
1040    else:
1041        SgP = ''
1042        varnCFs = ['lon', 'lat', 'time']
1043
1044    # CF dimension-variable dictionary
1045    CFdimvardict = {'lon': 'lon', 'lat': 'lat', 'pres': 'pres', 'time': 'time'}
1046
1047    # Input file
1048    ifilen = odir + '/' + CFvarn + '_' + headerf + SgP + '.nc'
1049
1050    # Computing variable statisitcs
1051    istats=0
1052
1053    # List of not computed statistics
1054    wrongstats = []
1055
1056    opers = Opers.split('+')
1057    Fopers = ''
1058    if db: print "    computing statistics of variable:'", CFvarn, "' operation '" + \
1059      Opers + "'...'"
1060    for op in opers:
1061
1062        # File name from previous operations, and name of the variable within the
1063        #   file (some operations change the name of the variable)
1064        if op == opers[0]:
1065            Fopers = op
1066            prevfile = ifilen
1067            vninF = CFvarn
1068        else:
1069            Fopers = Fopers + '_' + op
1070            prevfile = fileon
1071        fileon = odir + '/' + CFvarn + '_' + headerf + SgP + '_' + Fopers + '.nc'
1072
1073        # Adding required variables for the vertical interpolation in all operations
1074        SvarnCFs = ',' + ','.join(varnCFs)
1075
1076        if gP != 'none':
1077            if gP == 'local':
1078                CFvarnp = vninF + ',P,PB,PSFC,PH,PHB,HGT,T,QVAPOR,XLONG,XLAT,Times'+ \
1079                  SvarnCFs
1080            else:
1081                CFvarnp = vninF + SvarnCFs
1082        else:
1083            CFvarnp = vninF + SvarnCFs
1084
1085        if scr:
1086            sout = sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
1087
1088        hprevfile = prevfile[0:len(prevfile)-3]
1089
1090        if db: print '      ', op, "using '" + prevfile + "' "
1091        # Just produce file in case output file does not exist
1092        if not os.path.isfile(fileon):
1093            if db: 
1094                print '  ' + fname + ": creation of file '" + fileon + "' ..."
1095            dofile = True
1096        else:
1097            dofile = False
1098       
1099        # Variables to be kept in the final file
1100        varkeep = []
1101
1102        if op == 'acc':
1103            # temporal accumulated values
1104            print "  " + fname + ": kind '" + op + "' not ready !!"
1105            wrongstats.append(CFvarn + '_' + opers)
1106            break
1107        elif op == 'direct':
1108            # no statistics
1109            if dofile:
1110                sout = sub.call('mv ' + prevfile + ' ' + fileon, shell=True)
1111        elif op == 'last':
1112            # last temporal value
1113            vals='time,-9,0,0'
1114            if dofile:
1115                try:
1116                    with gen.Capturing() as output:
1117                        ncvar.DataSetSection(vals,prevfile)
1118                except:
1119                    print errmsg
1120                    print 'DataSetSection('+vals+',', prevfile, ')'
1121                    for s1out in output: print s1out
1122                    quit(-1)
1123
1124                if db:
1125                    for s1out in output: print s1out
1126               
1127                sout = sub.call('mv ' + hprevfile + '_time_B-9-E0-I0.nc ' + fileon, \
1128                  shell=True)
1129
1130                # Keeping the operations
1131                pyins=pyH + "/nc_var.py -o DataSetSection -S '" + vals + "' -f " +  \
1132                  prevfile
1133                otrackf.write("\n")
1134                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1135                otrackf.write(pyins + "\n")
1136
1137        elif op =='Lmean':
1138            # latitudinal mean values
1139            print "  " + fname + ": kind '" + op + "' not ready !!"
1140            wrongstats.append(CFvarn + '_' + opers)
1141            break
1142        elif op =='Lsec':
1143            # latitudinal section (latitudinal value must be given, [var]@[lat])
1144            print "  " + fname + ": kind '" + op + "' not ready !!"
1145            wrongstats.append(CFvarn + '_' + opers)
1146            break
1147        elif op =='lmean':
1148            # longitudinal mean values
1149            print "  " + fname + ": kind '" + op + "' not ready !!"
1150            wrongstats.append(CFvarn + '_' + opers)
1151            break
1152        elif op =='lsec':
1153            # longitudinal section (longitudinal value must be given, [var]@[lon])
1154            print "  " + fname + ": kind '" + op + "' not ready !!"
1155            wrongstats.append(CFvarn + '_' + opers)
1156            break
1157        elif op == 'pinterp':
1158            # pinterp: pressure interpolation (to the given $plevels)
1159            vals=plevels + ',1,1'
1160            if dofile:
1161                try:
1162                    with gen.Capturing() as output:
1163                        ncvar.pinterp(vals,prevfile,vninF)
1164                except:
1165                    print errmsg
1166                    print 'pinterp('+vals+',', prevfile, ','+vninF+')'
1167                    for s1out in output: print s1out
1168                    quit(-1)
1169
1170                if db:
1171                    for s1out in output: print s1out
1172
1173                sout = sub.call('mv pinterp.nc ' + fileon, shell=True)
1174
1175                # Keeping the operations
1176                pyins=pyH + "/nc_var.py -o pinterp -S '" + vals + "' -f " +          \
1177                  prevfile + "-v " + vninF
1178                otrackf.write("\n")
1179                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1180                otrackf.write(pyins + "\n")
1181
1182                # adding CF lon,lat,time in WRF files
1183                if minf.name == 'WRF':
1184                    values = vdnx + ':' + vdny + ':'+ Tref + ':' + Tunits
1185                    ncvar.WRF_toCF(values, fileon)
1186
1187            # vertical interpolation variables are no more needed
1188            CFvarnp = vninF
1189            gP = 'none'
1190
1191        elif op == 'tmean':
1192            # temporal mean values
1193            vals='time|-1,time,mean,' + ':'.join(varnCFs) + ':' + vdnz
1194            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1195            if dofile:
1196                try:
1197                    with gen.Capturing() as output:
1198                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1199                except:
1200                    print errmsg
1201                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1202                    for s1out in output: print s1out
1203                    quit(-1)
1204
1205                if db:
1206                    for s1out in output: print s1out
1207
1208                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)
1209
1210                # Keeping the operations
1211                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1212                  "' -f " + prevfile + " -v " + CFvarnp
1213                otrackf.write("\n")
1214                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1215                otrackf.write(pyins + "\n")
1216
1217            # removing dimension variable-dimension 'time'
1218            varnCFs.remove('time')
1219
1220            varkeep.append('timestats')
1221
1222        elif op == 'tstd':
1223            # temporal standard deviation values
1224            vals='time|-1,time,std,' + ':'.join(varnCFs) + ':' + vdnz
1225            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1226            if dofile:
1227                try:
1228                    with gen.Capturing() as output:
1229                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1230                except:
1231                    print errmsg
1232                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1233                    for s1out in output: print s1out
1234                    quit(-1)
1235
1236                if db:
1237                    for s1out in output: print s1out
1238
1239                sout = sub.call('mv file_oper_alongdims_std.nc '+fileon, shell=True)
1240
1241                # Keeping the operations
1242                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1243                  "' -f " + prevfile + " -v " + CFvarnp
1244                otrackf.write("\n")
1245                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1246                otrackf.write(pyins + "\n")
1247
1248            # removing dimension variable-dimension 'time'
1249            varnCFs.remove('time')
1250
1251            varkeep.append('timestats')
1252
1253        elif op == 'tturb':
1254            # temporal turbulent values
1255            vals='time|-1,time,turb,' + ':'.join(varnCFs) + ':' + vdnz
1256            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1257            if dofile:
1258                try:
1259                    with gen.Capturing() as output:
1260                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1261                except:
1262                    print errmsg
1263                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1264                    for s1out in output: print s1out
1265                    quit(-1)
1266
1267                if db:
1268                    for s1out in output: print s1out
1269
1270                sout = sub.call('mv file_oper_alongdims_turb.nc '+fileon, shell=True)
1271
1272                # Keeping the operations
1273                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1274                  "' -f " + prevfile + " -v " + CFvarnp
1275                otrackf.write("\n")
1276                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1277                otrackf.write(pyins + "\n")
1278
1279            # removing dimension variable-dimension 'time'
1280            varkeep.append('timestats')
1281
1282        elif op == 'tvar':
1283            # temporal variance values (Taylor's turbulence)
1284            vals='time|-1,time,var,' + ':'.join(varnCFs) + ':' + vdnz
1285            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1286            if dofile:
1287                try:
1288                    with gen.Capturing() as output:
1289                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1290                except:
1291                    print errmsg
1292                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1293                    for s1out in output: print s1out
1294                    quit(-1)
1295
1296                if db:
1297                    for s1out in output: print s1out
1298
1299                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)
1300
1301                # Keeping the operations
1302                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1303                  "' -f " + prevfile + " -v " + CFvarnp
1304                otrackf.write("\n")
1305                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1306                otrackf.write(pyins + "\n")
1307
1308            # removing dimension variable-dimension 'time'
1309            varnCFs.remove('time')
1310
1311            varkeep.append('timestats')
1312
1313        elif op == 'xmean':
1314            # x-axis mean values
1315            vals='lon|-1,lon,mean,' + ':'.join(varnCFs) + ':' + vdnz
1316            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1317            if dofile:
1318                try:
1319                    with gen.Capturing() as output:
1320                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1321                except:
1322                    print errmsg
1323                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1324                    for s1out in output: print s1out
1325                    quit(-1)
1326
1327                if db:
1328                    for s1out in output: print s1out
1329
1330                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)
1331
1332                # Keeping the operations
1333                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1334                  "' -f " + prevfile + " -v " + CFvarnp
1335                otrackf.write("\n")
1336                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1337                otrackf.write(pyins + "\n")
1338
1339            # removing dimension variable-dimension 'lon'
1340            varnCFs.remove('lon')
1341
1342            varkeep.append('lonstats')
1343
1344        elif op == 'xvar':
1345            # x-axis variance values
1346            vals='lon|-1,lon,var,' + ':'.join(varnCFs) + ':' + vdnz
1347            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1348            if dofile:
1349                try:
1350                    with gen.Capturing() as output:
1351                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1352                except:
1353                    print errmsg
1354                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1355                    for s1out in output: print s1out
1356                    quit(-1)
1357
1358                if db:
1359                    for s1out in output: print s1out
1360
1361                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)
1362
1363                # Keeping the operations
1364                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1365                  "' -f " + prevfile + " -v " + CFvarnp
1366                otrackf.write("\n")
1367                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1368                otrackf.write(pyins + "\n")
1369
1370            # removing dimension variable-dimension 'lon'
1371            varnCFs.remove('lon')
1372
1373            varkeep.append('lonstats')
1374
1375        elif op == 'ymean':
1376            # y-axis mean values
1377            vals='lat|-1,lat,mean,' + ':'.join(varnCFs) + ':' + vdnz
1378            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
1379            if dofile:
1380                try:
1381                    with gen.Capturing() as output:
1382                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
1383                except:
1384                    print errmsg
1385                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
1386                    for s1out in output: print s1out
1387                    quit(-1)
1388
1389                if db:
1390                    for s1out in output: print s1out
1391
1392                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)
1393
1394                # Keeping the operations
1395                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
1396                  "' -f " + prevfile + " -v " + CFvarnp
1397                otrackf.write("\n")
1398                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
1399                otrackf.write(pyins + "\n")
1400
1401            # removing dimension variable-dimension 'lat'
1402            varnCFs.remove('lat')
1403
1404            varkeep.append('latstats')
1405
1406        elif op == 'zsum':
1407            # vertical aggregated values
1408            print "  " + fname + ": kind '" + op + "' not ready !!"
1409            wrongstats.append(CFvarn + '_' + opers)
1410            break
1411        else:
1412            print errmsg
1413            print '  ' + fname + ": operation '" + op + "' not ready !!"
1414            quit(-1)
1415
1416        # End of kind of operation
1417
1418        # Variable name in file (vninF) changed due to operation
1419        #   but only if previous operation does not the same 'statistic'
1420        chvn = gen.dictionary_key_list(Opsurs, op)
1421        if chvn is not None:
1422            oldvninF = vninF
1423            vninF = vninF + chvn
1424            CFvarnp = CFvarnp.replace(oldvninF,vninF)
1425               
1426        if dofile:
1427            if len(varkeep) > 0:
1428                varkeepS = ',' + ','.join(varkeep)
1429            else:
1430                varkeepS = ''
1431
1432            # Adding such CF dimension variables which might not be in the
1433            #   post-operation file
1434            try:
1435                with gen.Capturing() as output: 
1436                    varinfile = ncvar.ivars(fileon)
1437            except:
1438                print errmsg
1439                print 'ivars('+fileon+')'
1440                for s1out in output: print s1out
1441                quit(-1)
1442
1443            for CFvn in varnCFs:
1444                if not gen.searchInlist(varinfile, CFvn):
1445                    if db:
1446                        print '  ' + fname + ": recupering CF variable '" + CFvn + "'"
1447                    try:
1448                        with gen.Capturing() as output:
1449                            ncvar.fvaradd(prevfile+','+CFvn, fileon)
1450                    except:
1451                        print errmsg
1452                        print 'fvaradd(' + prevfile + ',' + CFvn +', '+fileon+')'
1453                        for s1out in output: print s1out
1454                        quit(-1)
1455
1456                    if db:
1457                        for s1out in output: print s1out
1458
1459            totalvarkeeps = CFvarnp + ',' + ','.join(varnCFs) + varkeepS
1460            try:
1461                with gen.Capturing() as output:
1462                    oclean = ncvar.cleaning_varsfile(totalvarkeeps,fileon)
1463            except:
1464                print errmsg
1465                print 'cleaning_varsfile('+totalvarkeeps+','+fileon+')'
1466                for s1out in output: print s1out
1467                quit(-1)
1468
1469            if db:
1470                for s1out in output: print s1out
1471
1472    # End of operations
1473    if db:
1474        print fname + ": successful calculation of '" + vninF + "' from '" +        \
1475          CFvarn + "' '" + Opers + "' !!"
1476
1477    if len(wrongstats) > 1:
1478        print warnmsg
1479        print '  ' + fname + ": statisitcs not possible to compute:", wrongstats
1480
1481    return 
1482
1483def compute_vars(config, modinf, idir, odir, Files, allvarcomp, fscr, debug):
1484    """ Function to compute the variables
1485      config= Configuration of the experiment
1486      modinf= class with information about the model
1487      idir= input experiment folder
1488      odir= output experiment folder
1489      Files= dictionary with the header and the files' correspondence
1490      allvarcomp= dictionary with all the variables to compute and their information
1491      fscr= whether files should be done from the scratch or not
1492    """
1493    fname = 'compute_vars'
1494
1495    for vopn in allvarcomp.keys():
1496        # variable & operation
1497        vn = vopn.split('_')[0]
1498        op = vopn.split('_')[1]
1499        vopnV = allvarcomp[vopn]
1500        fheader = vopnV[0]
1501        if type(vopnV[1]) == type([1, 2]) and len(vopnV[1]) == 1:
1502            vals = vopnV[1]
1503            model = vals[0]
1504        else:
1505            model = vopnV[1]
1506        diag = vopnV[2]
1507        globalP = vopnV[3]
1508
1509        if debug:
1510            print '      ' + vn + ': ' + op
1511
1512        # Computing CF variable
1513        if model is not None or diag is not None:
1514            vinf = VariableInf(vn,fheader,model,diag)
1515            if op.find('pinterp') != -1: 
1516                SP = 'p'
1517            else:
1518                SP = ''
1519
1520            # Comppute variable
1521            finalfilen = odir + '/' + vn + '_' + fheader + SP + '.nc'
1522            if fscr:
1523                sout = sub.call('rm ' + finalfilen + ' >& /dev/null', shell=True)
1524
1525            if not os.path.isfile(finalfilen):
1526                compute_variable(modinf, idir, Files, odir, vinf, globalP, fscr,     \
1527                  config['pyHOME'], config['CFreftime'], config['CFunitstime'], debug)
1528
1529            # Compute variable statistics
1530            finalfilen = odir + '/' + vn + '_' + fheader + SP + '_' +                \
1531              op.replace('+','_') + '.nc'
1532            if fscr:
1533                sout = sub.call('rm ' + finalfilen + ' >& /dev/null', shell=True)
1534
1535            if not os.path.isfile(finalfilen):
1536                compute_statistics(modinf, cnf, iwdir, Files, owdir, vinf, globalP,  \
1537                  op, fscr, debug)
1538
1539        else:
1540            print errmsg
1541            print '  ' + fname + ": neither 'model' or 'diag' variables for '" + vn  \
1542              + "' !!"
1543            quit(-1)
1544
1545    return
1546
1547def get_plots_var(vc,pltk,db):
1548    """ Function to provide the plots to make for each variable from pltk ('DIRPLT_', or 'DIFFPLT_') dictionary
1549      vc= pltk dictionary, dictionary with all the parameters started with pltk
1550      pltk= kind of plots. Values started by 'DIRPLT_', or 'DIFFPLT_' in model configuration
1551    """
1552    fname = 'get_plots_var'
1553
1554    LH = len(pltk)
1555    # list of plots
1556    doplots = {}
1557    # plots by variable
1558    plotsvar = {}
1559    # individual plots by variable (not repeated)
1560    indivplotsvar = {}
1561
1562    # Variables with a plot which requires vertical interpolation
1563    varplotp = []
1564
1565    ipl = 0
1566    for plot in vc.keys():
1567        Lplot = len(plot)
1568        pln = plot[LH:Lplot+1]
1569        varops = vc[plot].split(':')
1570        for vn in varops:
1571            VOn = vn.split('#')
1572            # Variables and operations in plot
1573            vplt = []
1574            oplt = []
1575            voplt = []
1576            for vop in VOn:
1577                vv = vop.split('|')[0]
1578                oo = vop.split('|')[1]
1579                vplt.append(vv)
1580                oplt.append(oo)
1581                if not gen.searchInlist(voplt,vop): voplt.append(vop)
1582            vpltS = '#'.join(vplt)
1583            opltS = '#'.join(oplt)
1584           
1585            if not doplots.has_key(pln):
1586                doplots[pln] = [vn]
1587            else:
1588                plots = doplots[pln]
1589                plots.append(vn)
1590                doplots[pln] = plots
1591
1592            if not plotsvar.has_key(vn):
1593                plotsvar[vn] = [pln]
1594            else:
1595                plots = plotsvar[vn]
1596                plots.append(pln)
1597                plotsvar[vn] = plots
1598
1599            for VOv in voplt:
1600                if not indivplotsvar.has_key(VOv):
1601                    indivplotsvar[VOv] = [pln]
1602                else:
1603                    plots = indivplotsvar[VOv]
1604                    plots.append(pln)
1605                    indivplotsvar[VOv] = plots
1606
1607        if vc[plot].find('pinterp') != -1:
1608            if not gen.searchInlist(varplotp,vn): varplotp.append(pln)
1609
1610        ipl = ipl + 1
1611
1612    if db:
1613        print '  plots to make _______'
1614        gen.printing_dictionary(doplots)
1615        print '  plots by variable|operation _______'
1616        gen.printing_dictionary(plotsvar)
1617        print '  individual plots by variable|operation _______'
1618        gen.printing_dictionary(indivplotsvar)
1619        print '  plots with the vertical interpolation _______'
1620        print '    #', varplotp
1621
1622    return doplots, plotsvar, indivplotsvar, varplotp
1623
1624def plots_listconstruct(config, pkind, finf, odir, debug):
1625    """ Function to create the list of plots to draw
1626      config= Configuration of the experiment
1627      pkind= kind of plots ('DIRPLT', 'DIFFPLT')
1628      finf= file with the instructions of the plots
1629      odir= output experiment folder
1630    """
1631    fname='plots_listconstruct'
1632 
1633    dirplot = gen.get_specdictionary_HMT(config, H=pkind+'_')
1634
1635    # No plots because there are not plots of this kind
1636    if dirplot is None: 
1637        print warnmsg
1638        print '  ' + fname + ": no plots for type '" + pkind + "' !!"
1639        return {}, 0
1640
1641    if debug:
1642        print "  '" + pkind + "' plots to draw ________"
1643        gen.printing_dictionary(dirplot)
1644
1645    # Getting plots by variable
1646    plots, varplot, indivarplot, plotp = get_plots_var(dirplot, pkind+'_', debug)
1647
1648    Svarplot = gen.dictKeysVals_stringList(plots,cV=':')
1649
1650    Nplots = 0
1651    for pl in Svarplot.split(','):
1652        Nplots = Nplots + len(pl.split(':'))
1653
1654    # Outwritting the varcompute to avoid next time (if it is not filescratch!)
1655    objf = open(finf, 'w')
1656    objf.write('plots: ' + Svarplot + '\n')
1657    objf.write('itotp: ' + str(Nplots) + '\n')
1658    objf.close() 
1659
1660    return plots, Nplots
1661
1662def read_plot_file(figfile):
1663    """ Function to read the file with the information about the  plots to draw
1664      figfile= file with the information
1665    """
1666    fname = 'read_plot_file'
1667
1668    if not os.path.isfile(figfile):
1669        print errormsg
1670        print '  ' + fname + ": figures file '" + figfile + "' does not exist !!"
1671        quit(-1)
1672
1673    objf = open(figfile, 'r')
1674    for line in objf:
1675        if line[0:1] != '#' and len(line) > 1:
1676            values = line.replace('\\n','').split(' ')
1677            if values[0] == 'plots:':
1678                if values[1] != 'none':
1679                    plots = gen.stringList_dictKeysVals(values[1],cV=':')
1680                else:
1681                    pots = {}
1682            elif values[0] == 'itotp:':
1683                Nplots = int(values[1])
1684
1685    objf.close()
1686
1687    return plots, Nplots
1688
1689def varnoper(vn, oper, opsurdict):
1690    """ Function to provide name of the variable after operation as result of addition of the 'surnames'
1691      vn= variable name
1692      oper= '+' separated list of operations
1693      opsurdict= dictionary with the surname of operations
1694    >>> OpSurDict = {'mean': ['tmean', 'xmean', 'ymean']}
1695    >>> varnoper('tas', 'pinterp+tmean+xmean', OpSurDict)
1696    tasmeanmean
1697    """
1698    fname = 'varnoper'
1699
1700    newvn = vn
1701
1702    for op in oper.split('+'):
1703        surname = gen.dictionary_key_list(opsurdict,op)
1704        if surname is not None:
1705            newvn = newvn + surname
1706
1707    return newvn
1708
1709def draw_plot(kplot, CFvplot, fplot, vplot, dplot, pplot, finame, tfig, kfig, mapval,\
1710  tvals, expgraphc, od, pyH, fscr, db):
1711    """ Function to draw a plot
1712      kplot= kind of plot
1713      CFvplot= list of the CFnames of the variables
1714      fplot= list with the files in the plot
1715      vplot= list with the name of the variables in each file
1716      dplot= list of the dimensions in each file
1717      pplot= list of pictoric characteristics for each variable
1718      finame= name of the figure
1719      tfig= title of the figure
1720      kfig= kind of the figure
1721      mapval= value of the map
1722      tvals= list with time values
1723        tunits: units of time ('!' for spaces)
1724        timekind: kind of time ticks in the plot
1725        timefmt: format of time ticks in the plot
1726        timelabel: label of time-axis in the plot
1727      expgraphc= dictionary with expriment graphical characteristics
1728      od= output directory
1729      pyH= python HOME
1730      fscr= whether figure should be done from scratch
1731    """
1732    fname = 'draw_plot'
1733
1734    if fscr:
1735        sout = sub.call('rm ' + finame + ' >& /dev/null', shell=True)
1736
1737    os.chdir(od)
1738
1739    # CF variable-dimensions
1740    CFvardims = {'lon':'lon', 'lat': 'lat', 'pres': 'pres', 'time': 'time'}
1741
1742    # Tracking file with all the figures
1743    trkobjf = open('all_figures.inf', 'a')
1744
1745    # Time-characteristics
1746    tunits = tvals[0]
1747    timekind = tvals[1]
1748    timefmt = tvals[2]
1749    timelabel = tvals[3]
1750
1751    # Graphical labels and configuration
1752   
1753
1754    if not os.path.isfile(finame):
1755        if db:
1756            print "    drawing '" + finame + "' ..."
1757            print '      plot kind:', kplot
1758            print '      CFvars:', CFvplot
1759            print '      files:', fplot
1760            print '      in file vars:', vplot
1761            print '      dims:', dplot
1762            print '      pictoric values:', pplot
1763            print '      figure name:', finame
1764            print '      title:', tfig.replace('!',' ')
1765            print '      kind:', kfig
1766            print '      map:', mapval
1767            print '      time units:', tunits
1768
1769        if kplot == 'diffmap2Dsfc':
1770            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
1771            quit(-1)
1772        elif kplot == 'diffmap2Dz':
1773            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
1774            quit(-1)
1775        elif kplot == 'map2Dsfc':
1776            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
1777            quit(-1)
1778        elif kplot == 'map3D':
1779            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
1780            quit(-1)
1781        elif kplot == '2lines':
1782            figtit = tfig.replace('!','|')
1783
1784            lineAstdn = CFvplot[0]
1785            lineBstdn = CFvplot[1]
1786            figfs = ','.join(fplot)
1787
1788            lineA = pplot[0]
1789            Arange = str(lineA[0]) + ',' + str(lineA[1])
1790            # Values are changed from specific values in `model_graphics.dat'
1791            if lineA[3].find('%') == -1:
1792                Aline = lineA[2]
1793                Akind = lineA[3]
1794                Amark = lineA[4]
1795                Asize = lineA[5]
1796            else:
1797                Aline = 'red'
1798                Akind = '-'
1799                Amark = ','
1800                Asize = 2.
1801
1802            lineB = pplot[0]
1803            Brange = str(lineB[0]) + ',' + str(lineB[1])
1804            # Values are changed from specific values in `model_graphics.dat'
1805            if lineB[3].find('%') == -1:
1806                Bline = lineB[2]
1807                Bkind = lineB[3]
1808                Bmark = lineB[4]
1809                Bsize = lineB[5]
1810            else:
1811                Bline = 'blue'
1812                Bkind = '-'
1813                Bmark = ','
1814                Bsize = '2.'
1815
1816            # It is assumed that if the space variable is 'lon': is desired a
1817            #   (lon, vals) plot if it is 'lat': then (vals, lat) plot
1818            if gen.searchInlist(dplot[0],'lon'):
1819                spacedim ='lon'
1820                axisvals = 'y'
1821                figmid = 'longitudinal'
1822            elif gen.searchInlist(dplot[0],'lat'):
1823                spacedim = 'lat'
1824                axisvals = 'x'
1825                figmid = 'meridional'
1826            elif gen.searchInlist(dplot[0],'pres'):
1827                spacedim = 'pres'
1828                axisvals = 'x'
1829                figmid = 'vertical'
1830            else:
1831                print errmsg
1832                print '  ' + fname + ": in '2lines' only ready for: 'lon', 'lat'," + \
1833                 " 'pres' as common dimension !!"
1834                print '  dimensions in plot:', dplot[0]
1835                quit(-1)
1836
1837            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
1838              '!'.join(tfig.split('!')[2:])
1839
1840            graphvals = spacedim + ':' + Arange + ':' + Brange + ':Extrs:' +         \
1841              axisvals + ':' + ','.join(CFvplot) + ':' + Aline + ',' + Bline +       \
1842              ':2.,2.:' + Akind + ',' + Bkind + ':2.,2.:,:' + figtit + ':' +         \
1843              spacedim + ':0|auto:' +  finame.replace('.'+kfig, '') + ':' + kfig +   \
1844              ':True'
1845
1846            fvarS = ','.join(vplot)
1847
1848            drwins = 'draw_2lines'
1849            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
1850              " -S '" + graphvals + "' -v " + fvarS
1851
1852            try:
1853                with gen.Capturing() as output:
1854                    sout = sub.call(plotins, shell=True)
1855            except:
1856                print errmsg
1857                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
1858                print sout
1859                for s1out in output: print s1out
1860                quit(-1)
1861
1862            # keeping all figures
1863            trkobjf.write('\n')
1864            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
1865            trkobjf.write(plotins + '\n')
1866
1867        elif kplot == '2linesTime':
1868            figtit = tfig.replace('!','|')
1869
1870            lineAstdn = CFvplot[0]
1871            lineBstdn = CFvplot[1]
1872            figfs = ','.join(fplot)
1873
1874            lineA = pplot[0]
1875            Arange = str(lineA[0]) + ',' + str(lineA[1])
1876            # Values are changed from specific values in `model_graphics.dat'
1877            if lineA[3].find('%') == -1:
1878                Aline = lineA[2]
1879                Akind = lineA[3]
1880                Amark = lineA[4]
1881                Asize = lineA[5]
1882            else:
1883                Aline = 'red'
1884                Akind = '-'
1885                Amark = ','
1886                Asize = '2.'
1887
1888            lineB = pplot[1]
1889            Brange = str(lineB[0]) + ',' + str(lineB[1])
1890            # Values are changed from specific values in `model_graphics.dat'
1891            if lineB[3].find('%') == -1:
1892                Bline = lineB[2]
1893                Bkind = lineB[3]
1894                Bmark = lineB[4]
1895                Bsize = lineB[5]
1896            else:
1897                Bline = 'blue'
1898                Bkind = '-'
1899                Bmark = ','
1900                Bsize = '2.'
1901
1902            timeaxis = 'x'
1903            figmid = 'evolution'
1904
1905            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
1906              '!'.join(tfig.split('!')[2:])
1907
1908            graphvals = 'time:' + Arange + ':' + Brange + ':' + timekind + ';' +     \
1909              timefmt + ':' + timeaxis + ':' + ','.join(CFvplot) + ':' + Aline + ',' \
1910              + Bline + ':2.,2.:' + Akind + ',' + Bkind + ':2.,2.:,:' + figtit + ':' \
1911              + timelabel + ':0|auto:' +  finame.replace('.'+kfig, '') + ':' +       \
1912              kfig + ':True' 
1913
1914            fvarS = ','.join(vplot)
1915
1916            drwins = 'draw_2lines_time'
1917            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
1918              " -S '" + graphvals + "' -v " + fvarS
1919
1920            try:
1921                with gen.Capturing() as output:
1922                    sout = sub.call(plotins, shell=True)
1923            except:
1924                print errmsg
1925                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
1926                print sout
1927                for s1out in output: print s1out
1928                quit(-1)
1929
1930            # keeping all figures
1931            trkobjf.write('\n')
1932            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
1933            trkobjf.write(plotins + '\n')
1934
1935        elif kplot == 'Nlines':
1936            figtit = tfig.replace('!','|')
1937
1938            linestdn = CFvplot[0]
1939            figfs = ','.join(fplot)
1940
1941            Nlines = len(fplot)
1942
1943            lablines = []
1944            for il in range(Nlines):
1945                linev = pplot[il]
1946                # Values are changed from specific values in `model_graphics.dat'
1947                if linev[3].find('%') == -1:
1948                    line = linev[2]
1949                    kind = linev[3]
1950                    mark = linev[4]
1951                    size = linev[5]
1952                    if il == 0:
1953                        vrange = str(linev[0]) + ',' + str(linev[1])
1954                        kinds = kind
1955                        lines = line
1956                        marks = mark
1957                        sizes = str(size)
1958                    else:
1959                        kinds = kinds + ',' + kind
1960                        lines = lines + ',' + line
1961                        marks = marks + '@' + mark
1962                        sizes = sizes + ',' + str(size)
1963                else:
1964                    lines = 'None'
1965                    kinds = '-'
1966                    marks = ','
1967                    sizes = '2.'
1968
1969                secsf = fplot[il].split('/')
1970                Nsecsf = len(secsf)
1971                expvs = expgraphc[secsf[Nsecsf-3] + '/' + secsf[Nsecsf-2]]
1972                lablines.append(expvs.label)
1973
1974            # It is assumed that if the space variable is 'lon': is desired a
1975            #   (lon, vals) plot if it is 'lat': then (vals, lat) plot
1976            if gen.searchInlist(dplot[0],'lon'):
1977                spacedim ='lon'
1978                sdlab ='lon ($degrees\ East$)'
1979                axisvals = 'y'
1980                figmid = 'longitudinal'
1981            elif gen.searchInlist(dplot[0],'lat'):
1982                spacedim = 'lat'
1983                sdlab = 'lat ($degrees\ North$)'
1984                axisvals = 'x'
1985                figmid = 'meridional'
1986            elif gen.searchInlist(dplot[0],'pres'):
1987                spacedim = 'pres'
1988                sdlab = 'pres ($Pa$)'
1989                axisvals = 'x'
1990                figmid = 'vertical'
1991            else:
1992                print errmsg
1993                print '  ' + fname + ": in 'Nlines' only ready for: 'lon', 'lat'," + \
1994                 " 'pres' as common dimension !!"
1995                print '  dimensions in plot:', dplot[0]
1996                quit(-1)
1997
1998            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
1999              '!'.join(tfig.split('!')[2:])
2000            figtit = figtit.replace('!',' ')
2001            leglabels = ','.join(lablines)
2002
2003            graphvals = spacedim + ':' + axisvals + ':' + sdlab + ':' + leglabels +  \
2004              ':' + CFvplot[0] +':'+ figtit + ':0|auto:' + lines +':'+ kinds + ':' + \
2005              marks +':'+ sizes +':'+ sizes +':all:'+ finame.replace('.'+kfig, '') + \
2006              ':' + kfig + ':True'
2007
2008            fvarS = vplot[0]
2009
2010            drwins = 'draw_lines'
2011            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
2012              " -S '" + graphvals + "' -v " + fvarS
2013
2014            try:
2015                with gen.Capturing() as output:
2016                    sout = sub.call(plotins, shell=True)
2017            except:
2018                print errmsg
2019                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
2020                print sout
2021                for s1out in output: print s1out
2022                quit(-1)
2023
2024            # keeping all figures
2025            trkobjf.write('\n')
2026            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
2027            trkobjf.write(plotins + '\n')
2028
2029        elif kplot == 'Nlines_time':
2030            figtit = tfig.replace('!','|')
2031
2032            linestdn = CFvplot[0]
2033            figfs = ','.join(fplot)
2034
2035            Nlines = len(fplot)
2036
2037            lablines = []
2038            for il in range(Nlines):
2039                linev = pplot[il]
2040                # Values are changed from specific values in `model_graphics.dat'
2041                if linev[3].find('%') == -1:
2042                    line = linev[2]
2043                    kind = linev[3]
2044                    mark = linev[4]
2045                    size = linev[5]
2046                    if il == 0:
2047                        vrange = str(linev[0]) + ',' + str(linev[1])
2048                        kinds = kind
2049                        lines = line
2050                        marks = mark
2051                        sizes = str(size)
2052                    else:
2053                        kinds = kinds + ',' + kind
2054                        lines = lines + ',' + line
2055                        marks = marks + '@' + mark
2056                        sizes = sizes + ',' + str(size)
2057                else:
2058                    lines = 'None'
2059                    kinds = '-'
2060                    marks = ','
2061                    sizes = '2.'
2062
2063                secsf = fplot[il].split('/')
2064                Nsecsf = len(secsf)
2065                expvs = expgraphc[secsf[Nsecsf-3] + '/' + secsf[Nsecsf-2]]
2066                lablines.append(expvs.label)
2067
2068            valsaxis = 'y'
2069            figmid = 'evolution'
2070
2071            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
2072              '!'.join(tfig.split('!')[2:])
2073            figtit = figtit.replace('!',' ')
2074            leglabels = ','.join(lablines)
2075
2076            graphvals = 'time;' + valsaxis + ';time;' + leglabels + ';' +            \
2077              CFvplot[0] + ';' + figtit + ';None;time|' + '|'.join(tvals) +          \
2078              ';0|auto;' + kfig + ';' + kinds + ';' + lines + ';' + marks + ';' +    \
2079              sizes + ';' + sizes + ';all;-1;True' 
2080
2081            fvarS = vplot[0]
2082
2083            drwins = 'draw_lines_time'
2084            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
2085              " -S '" + graphvals + "' -v " + fvarS
2086
2087            try:
2088                with gen.Capturing() as output:
2089                    sout = sub.call(plotins, shell=True)
2090            except:
2091                print errmsg
2092                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
2093                print sout
2094                for s1out in output: print s1out
2095                quit(-1)
2096
2097            sout = sub.call('mv lines_time.' + kfig + ' ' + finame, shell=True)
2098
2099            # keeping all figures
2100            trkobjf.write('\n')
2101            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
2102            trkobjf.write(plotins + '\n')
2103
2104        elif kplot == 'shadcont2Dsfc':
2105            figtit = tfig.replace('!','|')
2106            shdstdn = CFvplot[0]
2107            cntstdn = CFvplot[1]
2108            figfs = ','.join(fplot)
2109
2110            shad = pplot[0]
2111            srange = str(shad[0]) + ',' + str(shad[1])
2112            cbar = shad[2]
2113            cnt = pplot[1]
2114            crange = str(cnt[0]) + ',' + str(cnt[1])
2115            cfmt = cnt[3]
2116            ckind = cnt[4]
2117            cline = cnt[5]
2118
2119            graphvals = ','.join(CFvplot)
2120            graphvals = graphvals + ':lon|-1,lat|-1:lon|-1,lat|-1:lon:lat:auto:' +   \
2121              cbar + ',auto,auto:' + ckind + ',' + cline + ':' + cfmt + ':' +        \
2122              srange + ':' + crange + ',9:' + figtit + ':' + kfig + ':None:' +       \
2123              mapval + ':True'
2124
2125            fvarS = ','.join(vplot)
2126
2127            drwins = 'draw_2D_shad_cont'
2128            plotins = "python " + pyH + "/drawing.py -f " + figfs + " -o " + drwins +\
2129              " -S '" + graphvals + "' -v " + fvarS
2130            try:
2131                with gen.Capturing() as output:
2132                    sout = sub.call(plotins, shell=True)
2133            except:
2134                print errmsg
2135                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
2136                print sout
2137                for s1out in output: print s1out
2138                quit(-1)
2139
2140            sout = sub.call('mv 2Dfields_shadow-contour.'+kfig+' '+finame, shell=True)
2141
2142            # keeping all figures
2143            trkobjf.write('\n')
2144            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
2145            trkobjf.write(plotins + '\n')
2146
2147        elif kplot == 'shadconthovmsfc':
2148            figtit = tfig.replace('!','|')
2149            shdstdn = CFvplot[0]
2150            cntstdn = CFvplot[1]
2151            figfs = ','.join(fplot)
2152
2153            shad = pplot[0]
2154            srange = str(shad[0]) + ',' + str(shad[1])
2155            cbar = shad[2]
2156            cnt = pplot[1]
2157            crange = str(cnt[0]) + ',' + str(cnt[1])
2158            cfmt = cnt[3]
2159            ckind = cnt[4]
2160            cline = cnt[5]
2161            # It is assumed that if the space variable is 'lon': is desired a
2162            #   (lon, time) plot if it is 'lat': then (time, lat) plot
2163            if gen.searchInlist(dplot,'lon'):
2164                spacedim ='lon'
2165                figmid = 'longitudinal|evolution|of'
2166                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
2167                  '|'.join(tfig.split('!')[2:])
2168                reverse= 'None'
2169                dims= spacedim+'|-1,time|-1;'+spacedim+'|-1,time|-1;'+spacedim+';time'
2170            else:
2171                spacedim='lat'
2172                figmid = 'meridional|evolution|of'
2173                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
2174                  '|'.join(tfig.split('!')[2:])
2175                reverse='None'
2176                dims= spacedim+'|-1,time|-1;'+spacedim+'|-1,time|-1;time;'+spacedim
2177
2178            graphvals = ','.join(CFvplot)
2179            graphvals = graphvals + ';' + dims +';auto;' + cbar + ',auto,auto' + ';' \
2180              + ckind + ',' + cline + ';' + cfmt + ';' + srange + ';' + crange +     \
2181              ',9;' + figtit + ';' + kfig + ';' + reverse + ';' + 'time|' + tunits + \
2182              '|'+ timekind + '|' + timefmt + '|' + timelabel + ';True'
2183
2184            fvarS = ','.join(vplot)
2185
2186            drwins = 'draw_2D_shad_cont_time'
2187            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
2188              " -S '" + graphvals + "' -v " + fvarS
2189
2190            try:
2191                with gen.Capturing() as output:
2192                    sout = sub.call(plotins, shell=True)
2193            except:
2194                print errmsg
2195                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
2196                print sout
2197                for s1out in output: print s1out
2198                quit(-1)
2199
2200            sout = sub.call('mv 2Dfields_shadow-contour.'+kfig+' '+finame, shell=True)
2201
2202            # keeping all figures
2203            trkobjf.write('\n')
2204            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
2205            trkobjf.write(plotins + '\n')
2206
2207        elif kplot == 'shadcont2Dzsec':
2208            figtit = tfig.replace('!','|')
2209            shdstdn = CFvplot[0]
2210            cntstdn = CFvplot[1]
2211            figfs = ','.join(fplot)
2212
2213            shad = pplot[0]
2214            srange = str(shad[0]) + ',' + str(shad[1])
2215            cbar = shad[2]
2216            cnt = pplot[1]
2217            crange = str(cnt[0]) + ',' + str(cnt[1])
2218            cfmt = cnt[3]
2219            ckind = cnt[4]
2220            cline = cnt[5]
2221
2222            # It is assumed that if the space variable is 'lon': is desired a
2223            #   (lon, pres) plot it it is 'lat': then (pres, lat) plot
2224            if gen.searchInlist(dplot, 'lon'):
2225                spacedim='lon'
2226                figmid = 'long.|vert.|cross|sec.|of'
2227                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
2228                  '|'.join(tfig.split('!')[2:])         
2229                dims= spacedim+'|-1,pres|-1:'+spacedim+'|-1,pres|-1:'+spacedim+':pres'
2230            else:
2231                spacedim='lat'
2232                figmid = 'mer.|vert.|cross|sec.|of'
2233                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
2234                  '|'.join(tfig.split('!')[2:])
2235# NOT WORKING?   dims= spacedim+'|-1,pres|-1:'+spacedim+'|-1,pres|-1:'+'pres:'+spacedim
2236                dims= spacedim+'|-1,pres|-1:'+spacedim+'|-1,pres|-1:'+spacedim+':pres'
2237
2238            reverse = 'flip@y'
2239
2240            graphvals = ','.join(CFvplot)
2241            graphvals = graphvals + ':' + dims + ':auto:' + cbar + ',auto,auto:' +   \
2242              ckind + ',' + cline + ':' + cfmt + ':' + srange + ':' + crange + ',9:' \
2243              + figtit + ':' + kfig + ':' + reverse + ':None:True'
2244 
2245            fvarS = ','.join(vplot)
2246
2247            drwins = 'draw_2D_shad_cont'
2248            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs + ' -o ' + drwins + \
2249              " -S '" + graphvals + "' -v " + fvarS
2250            try:
2251                with gen.Capturing() as output:
2252                    sout0 = sub.call(plotins, shell=True)
2253            except:
2254                print errmsg
2255                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
2256                for s1out in output: print s1out
2257                quit(-1)
2258
2259            sout = sub.call('mv 2Dfields_shadow-contour.'+kfig+' '+finame, shell=True)
2260
2261            # keeping all figures
2262            trkobjf.write('\n')
2263            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
2264            trkobjf.write(plotins + '\n')
2265
2266        else:
2267            print errmsg
2268            print "  " + fname + ": plot kind '" + kplot + "' not ready !!"
2269            quit(-1)
2270
2271    trkobjf.close()
2272
2273    return
2274
2275def optitle(op,opexp):
2276    """ Function to creat the operation section in the title of a figure ('!' for spaces)
2277      op= '+' separated list of operations
2278      opexp= dictionary with the `explanation'(text in title) to appear for each operation
2279    >>> optitle('pinterp+xmean+tmean',{'pinterp':'pinterp', 'xmean':'xmean', 'tmean':'tmean'})
2280    $_{[pinterp\!xmean\!tmean]}$
2281    """
2282    fname = 'optitle'
2283
2284    opvals = op.split('+')
2285    for op1 in opvals:
2286        if not opexp.has_key(op1):
2287            print errmsg
2288            print '  '+fname+": no explanation for operation '"+op1+"' !!"
2289            print '    provided:', opexp.keys()
2290            quit(-1)
2291        if op1 == opvals[0]:
2292            titop = '$_{[' + opexp[op1].replace(' ','!')
2293            if len(opvals) == 1: titop = titop + ']}$'
2294        elif op1 == opvals[len(opvals)-1]:
2295            titop = titop + '\!' + opexp[op1].replace(' ','!') + ']}$'
2296        else:
2297            titop = titop + '\!' + opexp[op1].replace(' ','!')
2298
2299    return titop
2300
2301def create_figure_title(mod,exp,varops,explop):
2302    """ Function to create the title of a figure ('!' for spaces)
2303      mod: model name
2304      exp: experiment name
2305      varops: list with [var]|[op] values from wich create the title
2306      explop: dictionary with the `explanation'(text in title) to appear for each operation
2307    >>> expops = {'pinterp':'pinterp', 'xmean':'xmean', 'tmean':'tmean', 'last':'last'}
2308    >>> create_figure_title('WRF','current',['ua|pinterp+xmean+tmean', 'va|pinterp+xmean+tmean'], expops)
2309    WRF!current!ua!&!va$_{[pinterp\!xmean\!tmean]}$
2310    >>> create_figure_title('WRF','current',['ua|pinterp+xmean+tmean', 'ua|pinterp+xmean+last'], expops)
2311    WRF!current!ua$_{[pinterp\!xmean\!tmean]}$!&!!$_{[pinterp\!xmean\!last]}$
2312    >>> create_figure_title('WRF','current',['ua|pinterp+xmean+tmean', 'va|pinterp+xmean+last'], expops)
2313    WRF!current!ua$_{[pinterp\!xmean\!tmean]}$!&!va$_{[pinterp\!xmean\!last]}$
2314    """
2315    fname = 'create_figure_title'
2316
2317    titfig = mod + '!' + exp + '!'
2318
2319    varns = []
2320    opns = []
2321    for varop in varops:
2322        vn = varop.split('|')[0]
2323        op = varop.split('|')[1]
2324        varns.append(vn)
2325        opns.append(op)
2326        if varop == varops[0]:
2327            allvarequal = True
2328            allopequal = True
2329            varinfig = vn
2330            opinfig = op
2331        else:
2332            if vn != varinfig: allvarequal = False
2333            if op != opinfig: allopequal = False
2334
2335    Nvars = len(varns)
2336    Nops = len(opns)
2337
2338    if allvarequal:
2339        titfig = titfig + '!' + varns[0]
2340        if allopequal:
2341            opS = optitle(op,explop)
2342            titfig = titfig + opS
2343        else:
2344            for opn in opns:
2345                opS = optitle(opn,explop)
2346                if opn == opns[0]:
2347                    titfig = titfig + opS
2348                elif opn == opns[Nops-1]:
2349                    titfig = titfig + '!&!' + opS
2350                else:
2351                    titfig = titfig + ',!' + opS
2352    else: 
2353        if allopequal:
2354            opS = optitle(op,explop)
2355            for vn in varns:
2356                if vn == varns[0]:
2357                    titfig = titfig + '!' + vn
2358                elif vn == varns[Nvars-1]:
2359                    titfig = titfig + '!&!' + vn
2360                else:
2361                    titfig = titfig + ',!' + vn
2362            titfig = titfig + opS
2363        else:
2364            for ivop in range(Nvars):
2365                vn = varns[ivop]
2366                op = opns[ivop]
2367                vnop = vn + optitle(op,explop)
2368               
2369                if vn == varns[0]:
2370                    titfig = titfig + '!' + vnop
2371                elif vn == varns[Nvars-1]:
2372                    titfig = titfig + '!&!' + vnop
2373                else:
2374                    titfig = titfig + ',!' + vnop
2375
2376    return titfig
2377
2378def draw_plots(config, plots, mod, exp, odir, allvarcomp, figscr, debug):
2379    """ Function to draw all plots
2380      config= Configuration of the experiment
2381      plots= dictionary with the plots
2382      mod= model of the plots
2383      exp= experiment of the plots
2384      odir= output experiment folder
2385      allvarcomp= dictionary with all the variables to compute and their information
2386      figscr= whether figures should be done from the scratch or not
2387
2388    * Plot as
2389      {[kplot]} = [varn1]|[op1]#[varn2]|[op2]#[...[varnN]|[opN]], ...
2390        [kplot] ___
2391          diffmap2Dsfc: 2D map of surface differences values of 1 variable
2392          diffmap2Dz: 2D map of 3D differences values of 1 variable
2393          map2Dsfc: 2D map of surface values of 1 variable
2394          map3D: 2D map of 3D values of 1 variable
2395          shadconthovmsfc: Hovmoeller diagrams of 2 variable at the surface in shadow and the other in contourn
2396          shadcont2Dsfc: 2D map of shadow (1st variable) and countour (2nd variable) [stvar1]#[stvar2]
2397          shadcont2Dzsec: 2D map of vertical section of 2 variables one in shadow and the other in contourn
2398        [varn]
2399          variable
2400        [op]
2401          '+' separated list of operations
2402        in figures with more than 1 variable, use '#' to separate the [varn]|[op]
2403    """ 
2404    fname = 'draw_plots'
2405
2406    os.chdir(odir)
2407
2408    # Dictionary with the operations with surnames for the operated variable
2409    opersurnames = {}
2410    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
2411    for opsr in opsur.keys():
2412        opn = opsr.split('_')[1]
2413        vls = opsur[opsr].split(':')
2414        opersurnames[opn] = vls
2415
2416    # time values
2417    # Units time for the plots
2418    rd = config['CFreftime']
2419    tunits = config['CFunitstime'] + '!since!' + rd[0:4] + '-' + rd[4:6] + '-' +  \
2420      rd[6:8] + '!' + rd[8:10] + ':' + rd[10:12] + ':' + rd[12:14]
2421    # time ticks kind
2422    tkind = config['timekind']
2423    # time ticks format
2424    tfmt = config['timefmt']
2425    # time axis label
2426    tlab = config['timelabel']
2427    timevals = [tunits, tkind, tfmt, tlab]
2428
2429    # Dictionary of plot specificities
2430    specplotkeyn = 'specificvarplot'
2431    plotspecifics = {}
2432    if config.has_key(specplotkeyn):
2433        #   [minval]: minimum value
2434        #   [maxval]: minimum value
2435        #   [colorbar]: name of the colorbar (from matplotlib) to use
2436        #   [cntformat]: format of the contour labels
2437        #   [colorcnt]: color for the countor lines
2438        plotspecs = config[specplotkeyn].split(':')
2439        for pltspc in plotspecs:
2440            pltvls = pltspc.split('|')
2441            vn = pltvls[0]
2442            op = pltvls[1]
2443            fn = pltvls[2]
2444            plotspecifics[fn + '_' + vn + '_' + op] = pltvls[3:]
2445        if debug:
2446            print 'Specific values for plots _______'
2447            gen.printing_dictionary(plotspecifics)
2448
2449    # Kind of figures
2450    kindfigure = config['kindfig']
2451
2452    # Graphical labels and configuration
2453    modgc, expgc = graphical_characteristics(config, debug)
2454    modv = modgc[mod]
2455    modlab = modv.label
2456    expv = expgc[mod+'/'+exp]
2457    explab = expv.label
2458
2459    # Map value
2460    mapvalue = config['mapval']
2461
2462    # pythone scripts HOME
2463    pyHOME = config['pyHOME']
2464
2465    # Title-text of operations
2466    opexplained = {}
2467    optits = config['titleoperations'].split(':')
2468    for optit in optits:
2469        opn = optit.split('|')[0]
2470        opt = optit.split('|')[1]
2471        opexplained[opn] = opt
2472    if debug:
2473        print 'Titles for operations  _______'
2474        gen.printing_dictionary(opexplained)
2475
2476    for kplot in plots.keys():
2477        varsplt = plots[kplot]
2478        for varplt in varsplt:
2479            if debug:
2480                print "  printing '" + kplot + "' var ':" + varplt + "'..."
2481            varops = varplt.split('#')
2482
2483            # CF variables in plot
2484            CFvarsplot = []
2485            # Files in plot
2486            filesplot = []
2487            # Variables in plot within the files
2488            varsplot = []
2489            # Dims in figure
2490            dimsplot = []
2491            # pictoric values in figure
2492            pictplot = []
2493            # Name of the figure
2494            figname = ''
2495            # Title of the figure
2496            titfigure = ''
2497
2498            ivp = 0
2499            for varop in varops:
2500                vn = varop.split('|')[0]
2501                op = varop.split('|')[1]
2502
2503                # CF variables in plot
2504                CFvarsplot.append(vn)
2505 
2506                vnopS = vn + '_' + op
2507                if not allvarcomp.has_key(vnopS):
2508                    print errmsg
2509                    print '  ' + fname + ": no file for variable-operation '" +     \
2510                      vnopS + "' !!"
2511                vopvals = allvarcomp[vnopS]
2512                headf = vopvals[0]
2513                globalP = vopvals[3]
2514                gP = pinterpS(globalP)
2515
2516                filen = vn + '_' + headf + gP + '_' + op.replace('+','_') + '.nc'
2517                filesplot.append(filen)
2518                # Do we have processed the given variable?
2519                if not os.path.isfile(filen):
2520                    print warnmsg
2521                    print "  " + fname + ": there is no file for variable '" + varop \
2522                      + "' skiping it !!"
2523                    break
2524
2525                # Name of the variable inside the file
2526                vnsur = varnoper(vn, op, opersurnames)
2527                varsplot.append(vnsur)
2528
2529                # Dimensions in file
2530                try:
2531                    with gen.Capturing() as output:
2532                        dims = ncvar.idims(filen)
2533                except:
2534                    print errmsg
2535                    print 'ncvar.idims('+filen+')'
2536                    for s1out in output: print s1out
2537                    quit(-1)
2538
2539                dimsplot.append(dims)
2540
2541                # pictoric values for the figure
2542                Sfivaop = kplot + '_' + vn + '_' + op
2543                if plotspecifics.has_key(Sfivaop):
2544                    pictvals = plotspecifics[Sfivaop]
2545                else:
2546                    Vvals = gen.variables_values(vn)
2547                    pictvals = [Vvals[2], Vvals[3], Vvals[6], '%g', 'fixc', 'black']
2548
2549                pictplot.append(pictvals)
2550
2551                # Header of the name of the figure
2552                if ivp == 0:
2553                    figname = kplot +'_'+ vn + '_' + headf + '_' + op.replace('+','-')
2554                else:
2555                    figname = figname + '-'+vn+'_'+headf+'_'+op.replace('+','-')
2556
2557                ivp = ivp + 1
2558            # End of variable-operation
2559            figname = figname + '.' + kindfigure
2560
2561            # Title of figure
2562            titfigure = create_figure_title(modlab, explab, varops, opexplained)
2563
2564            draw_plot(kplot, CFvarsplot, filesplot, varsplot, dimsplot, pictplot,    \
2565              figname, titfigure, kindfigure, mapvalue, timevals, expgc, odir,       \
2566              pyHOME, figscr, debug)
2567
2568        # End of variables-operations
2569
2570    # End of kind of plots
2571
2572    return
2573
2574def get_differences_var(vc,kdiff,db):
2575    """ Function to provide the operations to make for each variable from kdiff ('DIFFOP_' or 'DIFFVAR_') dictionary
2576      vc= kdiff dictionary, dictionary with all the parameters started with kdiff
2577      kdiff= kind of differences 'op', 'var'
2578    """
2579    fname = 'get_differences_var'
2580
2581    if kdiff == 'op':
2582        diffn = 'DIFFOP'
2583    elif kdiff == 'var':
2584        diffn = 'DIFFVAR'
2585    else:
2586         print errmsg
2587         print '  ' + fname + ": differences '" + kdiff + "' not ready !!"
2588         quit(-1)
2589    LD = len(diffn + '_')
2590
2591    iop = 1
2592    # list of operations by difference
2593    doopers = {}
2594    # operations by variable
2595    operationsvar = {}
2596    # individual operations by variable (not repeated)
2597    indivoperationsvar = {}
2598
2599    # Variables with a calculation whic requires vertical interpolation to all the
2600    #   file (operations stating by [kdiff]'_pinterp')
2601    LVp = len(diffn + '_pinterp')
2602    varglobalp = []
2603
2604    for oper in vc.keys():
2605        Loper = len(oper)
2606        opn = oper[LD:Loper+1]
2607        vns = vc[oper].split(':')
2608        ops1 = opn.split('+')
2609        doopers[opn] = ops1
2610        for vn in vns:
2611            if not operationsvar.has_key(vn):
2612                operationsvar[vn] = [opn]
2613                indivoperationsvar[vn] = ops1
2614            else:
2615                opers = operationsvar[vn]
2616                opers.append(opn)
2617                operationsvar[vn] = opers
2618                opers1 = indivoperationsvar[vn]
2619                for op1 in ops1:
2620                    if not gen.searchInlist(opers1, op1):
2621                        opers1.append(op1)
2622                indivoperationsvar[vn] = opers1
2623
2624        if oper[0:LVp] == diffn + '_pinterp':
2625            for vn in vns:
2626                if not gen.searchInlist(varglobalp,vn): varglobalp.append(vn)
2627
2628        iop = iop + 1
2629
2630    if db:
2631        print '  operations to make _______'
2632        gen.printing_dictionary(doopers)
2633        print '  operations by variable _______'
2634        gen.printing_dictionary(operationsvar)
2635        print '  individual operations by variable _______'
2636        gen.printing_dictionary(indivoperationsvar)
2637        print '  variables with the vertical interpolation for all the file _______'
2638        print '    #', varglobalp
2639
2640    return doopers, operationsvar, indivoperationsvar, varglobalp
2641
2642def diffvarop_listconstruct(config, mods, exps, reprj, odir, debug):
2643    """ Function to create the list of differences to compute and draw
2644      config= Configuration of the experiment
2645      mods= list with the models to use
2646      exps= list with the experiments to use
2647      reprj= list with whether re-rpojection needs to be done
2648      odir= output differences folder
2649    """
2650    fname='diffvarop_listconstruct'
2651 
2652    diffop = gen.get_specdictionary_HMT(config, H='DIFFOP_')
2653    diffvar = gen.get_specdictionary_HMT(config, H='DIFFVAR_')
2654
2655    if debug:
2656        if diffop is not None:
2657            print "  'op' differences ________"
2658            gen.printing_dictionary(diffop)
2659        if diffvar is not None:
2660            print "  'var' differences ________"
2661            gen.printing_dictionary(diffvar)
2662
2663    # Getting variable characteristics from each model/experiment
2664    idir1 = config['ofold'] + '/' + mods[0] + '/' + exps[0]
2665    idir2 = config['ofold'] + '/' + mods[1] + '/' + exps[1]
2666
2667    fvarcomp = 'varcompute.inf'
2668    if reprj[0]: fvc1 = 'reproj_' + fvarcomp
2669    else: fvc1 = fvarcomp
2670
2671    if reprj[1]: fvc2 = 'reproj_' + fvarcomp
2672    else: fvc2 = fvarcomp
2673
2674    files1, testfiles1, allcompvar1, Nvar1= read_varcomp_file(idir1 + '/' + fvc1)
2675    files2, testfiles2, allcompvar2, Nvar2= read_varcomp_file(idir2 + '/' + fvc2)
2676
2677    diffops = {}
2678    if diffop is not None:
2679        # Getting 'op' differences characteristics
2680        doops, opvar, indivopvar, varglobalp = get_differences_var(diffop,'op',debug)
2681
2682        Ndiffop = 0
2683        for vn in opvar.keys():
2684            vnops = opvar[vn]
2685            for op in vnops:
2686                vnop = vn + '_' + op
2687                if not allcompvar1.has_key(vnop):
2688                    print warnmsg
2689                    print '  ' + fname + ": DIFFOP variable+operation '" + vnop +    \
2690                      "' in '" + mods[0] + '/' + exps[0] + "' was not computed !!"
2691                    print '    skipping it!'
2692                    print '    variables computed:', allcompvar1
2693                    break
2694                if not allcompvar2.has_key(vnop):
2695                    print warnmsg
2696                    print '  ' + fname + ": DIFFOP variable+operation '" + vnop +    \
2697                      "' in '" + mods[1] + '/' + exps[1] + "' was not computed !!"
2698                    print '    skipping it!'
2699                    print '    variables computed:', allcompvar2
2700                    break
2701                vals1 = allcompvar1[vnop]
2702                vals2 = allcompvar2[vnop]
2703
2704                headerf1 = vals1[0]
2705                headerf2 = vals2[0]
2706
2707                diffops[vnop] = [mods[0],mods[1],exps[0],exps[1],headerf1,headerf2]
2708
2709                Ndiffop = Ndiffop + 1       
2710    else:
2711        Ndiffop = 0
2712
2713    diffvars = {}
2714    if diffvar is not None:
2715        # Getting 'var' differences characteristics
2716        doops,opvar,indivopvar,varglobalp = get_differences_var(diffvar,'var',debug)
2717
2718        Ndiffvar = 0
2719        for vn in opvar.keys():
2720            vnops = opvar[vn]
2721            for op in vnops:
2722                vnop = vn + '_' + op
2723                if not allcompvar1.has_key(vnop):
2724                    print warnmsg
2725                    print '  ' + fname + ": DIFFVAR variable+operation '" + vnop +   \
2726                      "' in '" + mods[0] + '/' + exps[0] + "' was not computed !!"
2727                    print '    skipping it!'
2728                    print '    variables computed:', allcompvar1
2729                    break
2730                if not allcompvar2.has_key(vnop):
2731                    print warnmsg
2732                    print '  ' + fname + ": DIFFVAR variable+operation '" + vnop +   \
2733                      "' in '" + mods[1] + '/' + exps[1] + "' was not computed !!"
2734                    print '    skipping it!'
2735                    print '    variables computed:', allcompvar2
2736                    break
2737                vals1 = allcompvar1[vnop]
2738                vals2 = allcompvar2[vnop]
2739   
2740                headerf1 = vals1[0]
2741                headerf2 = vals2[0]
2742
2743                diffvars[vnop] = [mods[0],mods[1],exps[0],exps[1],headerf1,headerf2]
2744
2745                Ndiffvar = Ndiffvar + 1
2746    else:
2747        Ndiffvar = 0
2748
2749    if diffop is not None:
2750        Sopdiffs = gen.dictKeysVals_stringList(diffops)
2751    else:
2752        Sopdiffs = 'none'
2753    # Outwritting the op diffferences file to avoid next time
2754    objf = open(odir + '/diffop.inf', 'w')
2755    objf.write('differences: ' + Sopdiffs + '\n')
2756    objf.write('Ndiff: ' + str(Ndiffop) + '\n')
2757    objf.close() 
2758
2759    if diffvar is not None:
2760        Svardiffs = gen.dictKeysVals_stringList(diffvars)
2761    else:
2762        Svardiffs = 'none'
2763    # Outwritting the var diffferences file to avoid next time
2764    objf = open(odir + '/diffvar.inf', 'w')
2765    objf.write('differences: ' + Svardiffs + '\n')
2766    objf.write('Ndiff: ' + str(Ndiffvar) + '\n')
2767    objf.close() 
2768
2769    return diffops, diffvars, Ndiffop, Ndiffvar
2770
2771def read_diff_file(difffile):
2772    """ Function to read the file with the information about the differences
2773      difffile= file with the information
2774    """
2775    fname = 'read_diff_file'
2776
2777    if not os.path.isfile(difffile):
2778        print errormsg
2779        print '  ' + fname + ": differences file '" + difffile + "' does not exist !!"
2780        quit(-1)
2781
2782    objf = open(difffile, 'r')
2783    for line in objf:
2784        if line[0:1] != '#' and len(line) > 1:
2785            values = line.replace('\n','').split(' ')
2786            if values[0] == 'differences:':
2787                if values[1] != 'none':
2788                    diffs = gen.stringList_dictKeysVals(values[1])
2789                else:
2790                    diffs = {} 
2791            elif values[0] == 'Ndiff:':
2792                Ndiffs = int(values[1])
2793
2794    objf.close()
2795
2796    return diffs, Ndiffs
2797
2798def compute_op_diffs(config, diffallop, odir, diffscr, debug):
2799    """ Function to compute operation differences
2800      config= experiment configuration
2801      diffallop= dictonary with the differences to compute
2802      odir= output folder
2803      diffscr= wether it should be done from the scratch
2804    """
2805    fname = 'compute_op_diffs'
2806
2807    # output folder
2808    ofold = config['ofold']
2809
2810    # Home of the python scripts
2811    pyH = config['pyHOME']
2812
2813    #  opsur = operation surnames
2814    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
2815    Opsurs = {}
2816    for opk in opsur:
2817        opn = opk.split('_')[1]
2818        vals = opsur[opk]
2819        Opsurs[opn] = vals.split(':')
2820
2821    # Getting in working dir
2822    os.chdir(odir)
2823
2824    # CF dimensions
2825    CFdims = ['lon', 'lat', 'pres', 'time']
2826
2827    # Trakcing file
2828    trkobjf = open(odir + '/all_diffop.inf', 'a')
2829
2830    for vnop in diffallop.keys():
2831        vn = vnop.split('_')[0]
2832        op = vnop.split('_')[1]
2833        opS = op.replace('+','_')
2834
2835        # `p' header add related to 'pinterp'
2836        if opS.find('pinterp') != -1: SgP = 'p'
2837        else: SgP = ''
2838
2839        # Variable name within the file
2840        vninF = varnoper(vn, op, Opsurs)
2841
2842        vals = diffallop[vnop]
2843        mod1 = vals[0]
2844        mod2 = vals[1]
2845        exp1 = vals[2]
2846        exp2 = vals[3]
2847        headerf1 = vals[4]
2848        headerf2 = vals[5]
2849        modexpdiff = mod2 + '/' + exp2 + ' - ' + mod1 + '/' + exp1
2850        modexpdiffS = mod2 + '-' + exp2 + '_' + mod1 + '-' + exp1
2851
2852        ifile1 = ofold + '/' + mod1 + '/' + exp1 + '/' + vn + '_' + headerf1 + SgP + \
2853          '_' + opS + '.nc'
2854        ifile2 = ofold + '/' + mod2 + '/' + exp2 + '/' + vn + '_' + headerf2 + SgP + \
2855          '_' + opS + '.nc'
2856        difffilen = odir + '/' + vn + '_diffop_' + modexpdiffS + '_' + opS + '.nc'
2857
2858        if diffscr:
2859            sout = sub.call('rm ' + difffilen + ' >& /dev/null', shell=True)
2860
2861        if not os.path.isfile(difffilen):
2862            if debug:
2863                print "    Computing operation difference '" + op + "' of '" + vn +  \
2864                  "' for: '" + modexpdiff
2865                print "    var in file: '" + vninF + "'"
2866
2867            dfiles = 'add|' + ifile2 + '|' + vninF + ',sub|' + ifile1 + '|' + vninF
2868            values = '|'.join(CFdims) + '@' + dfiles
2869            pyins = 'python ' + pyH + "/nc_var.py -o compute_opersvarsfiles -S '" +  \
2870              values + "' -v " + vninF
2871
2872            try:
2873                with gen.Capturing() as output:
2874                    ncvar.compute_opersvarsfiles(values,vninF)
2875            except:
2876                print errmsg
2877                print 'ncvar.compute_opersvarsfiles(' + values + ',' + vninF + ')'
2878                for s1out in output: print s1out
2879                quit(-1)
2880
2881            sout = sub.call('mv opersvarsfiles_'+vninF+'.nc ' + difffilen, shell=True)
2882
2883            # keeping all differences
2884            trkobjf.write('\n')
2885            trkobjf.write("#" + vn + ' ' + op + '\n')
2886            trkobjf.write( pyins + '\n')
2887
2888    trkobjf.close()
2889
2890    return
2891
2892def compute_diff_statistics(Ecnf, filen, ofold, varn, Opers, scr, db):
2893    """ Function to compute statistics from var diff files
2894      Ecnf= dictionary with the configuration of the experiment
2895      filen= variable difference file to use
2896      ofold= directory to write the output files
2897      varn= name of the variable
2898      Opers= kind of operation: (as possible multiple consecutive combination of operations separated by '+'
2899        [calc1]+[calc2] will compute first [calc1] and then [calc2]
2900        acc: temporal accumulated values
2901        diff: differences between models
2902        direct: no statistics
2903        last: last temporal value
2904        Lmean: latitudinal mean values
2905        Lsec: latitudinal section (latitudinal value must be given, [var]@[lat])
2906        lmean: longitudinal mean values
2907        lsec: longitudinal section (longitudinal value must be given, [var]@[lat])
2908        pinterp: pressure interpolation (to the given $plevels)
2909        tmean: temporal mean values
2910        tstd: temporal standard deviation values
2911        tturb: Taylor's turbulence decomposition value (x - <x>) for time
2912        tvar: temporal variance values
2913        xmean: x-axis mean values
2914        xvar: x-axis variance values
2915        ymean: y-axis mean values
2916        zsum: vertical aggregated values
2917      scr= should it be done from the scratch?
2918    """
2919    fname='compute_diff_statistics'
2920
2921    # Full header of file
2922    Lfilen = len(filen)
2923    Hfile = filen[0:Lfilen-3]
2924
2925    # CF dimension-variables name
2926    OpersS = '+'.join(Opers)
2927    if OpersS.find('pinterp') != -1 or filen.find('pinterp') != -1:
2928        varnCFs = ['lon', 'lat', 'pres', 'time']
2929    else:
2930        varnCFs = ['lon', 'lat', 'time']
2931
2932    # Variable-dimension CF dictionary
2933    CFdimvardict = {'lon': 'lon', 'lat': 'lat', 'pres': 'pres', 'time': 'time'}
2934
2935    # Python scripts HOME
2936    pyH = Ecnf['pyHOME']
2937
2938    #  opsur = operation surnames
2939    opsur = gen.get_specdictionary_HMT(Ecnf, H='opsur_',M='',T='')
2940    Opsurs = {}
2941    for opk in opsur:
2942        opn = opk.split('_')[1]
2943        vals = opsur[opk]
2944        Opsurs[opn] = vals.split(':')
2945
2946    # Getting in working dir
2947    os.chdir(ofold)
2948
2949    # File to keep track of all the statistics of var differences
2950    otrackf = open(ofold + '/all_vardiffstatistics.inf', 'a')
2951
2952    if db: 
2953        print "    computing statistics of diff variable: '", varn ,"' operation '"+ \
2954      OpersS + "' using file '" + filen + " ...'"
2955
2956    for op in Opers:
2957        # File name from previous operations, and name of the variable within the
2958        #   file (some operations change the name of the variable)
2959        if op == Opers[0]:
2960            Fopers = op
2961            prevfile = filen
2962            vninF = varn
2963        else:
2964            Fopers = Fopers + '_' + op
2965            prevfile = fileon
2966        fileon = Hfile + '_' + Fopers + '.nc'
2967
2968        SvarnCFs = ',' + ','.join(varnCFs) 
2969        CFvarnp = vninF + SvarnCFs
2970
2971        if scr:
2972            sout = sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
2973
2974        # Just produce file in case output file does not exist
2975        if not os.path.isfile(fileon):
2976            if db: 
2977                print '  ' + fname + ": creation of file '" + fileon + "' ..."
2978            dofile = True
2979        else:
2980            dofile = False
2981       
2982        # Variables to be kept in the final file
2983        varkeep = []
2984
2985        if op == 'acc':
2986            # temporal accumulated values
2987            print "  " + fname + ": kind '" + op + "' not ready !!"
2988            wrongstats.append(varn + '_' + opers)
2989            break
2990        elif op == 'direct':
2991            # no statistics
2992            if dofile:
2993                sout = sub.call('mv ' + prevfile + ' ' + fileon, shell=True)
2994        elif op == 'last':
2995            # last temporal value
2996            vals='time,-9,0,0'
2997            if dofile:
2998                try:
2999                    with gen.Capturing() as output:
3000                        ncvar.DataSetSection(vals,prevfile)
3001                except:
3002                    print errmsg
3003                    print 'DataSetSection('+vals+',', prevfile, ')'
3004                    for s1out in output: print s1out
3005                    quit(-1)
3006
3007                if db:
3008                    for s1out in output: print s1out
3009
3010                hprevfile = prevfile[0:len(prevfile)-3]
3011
3012                sout = sub.call('mv ' + hprevfile + '_time_B-9-E0-I0.nc ' + fileon, \
3013                  shell=True)
3014
3015                # Keeping the operations
3016                pyins=pyH + "/nc_var.py -o DataSetSection -S '" + vals + "' -f " +  \
3017                  prevfile
3018                otrackf.write("\n")
3019                otrackf.write("# " + varn + " " + Fopers + "\n")
3020                otrackf.write(pyins + "\n")
3021
3022        elif op =='Lmean':
3023            # latitudinal mean values
3024            print "  " + fname + ": kind '" + op + "' not ready !!"
3025            wrongstats.append(varn + '_' + opers)
3026            break
3027        elif op =='Lsec':
3028            # latitudinal section (latitudinal value must be given, [var]@[lat])
3029            print "  " + fname + ": kind '" + op + "' not ready !!"
3030            wrongstats.append(varn + '_' + opers)
3031            break
3032        elif op =='lmean':
3033            # longitudinal mean values
3034            print "  " + fname + ": kind '" + op + "' not ready !!"
3035            wrongstats.append(varn + '_' + opers)
3036            break
3037        elif op =='lsec':
3038            # longitudinal section (longitudinal value must be given, [var]@[lon])
3039            print "  " + fname + ": kind '" + op + "' not ready !!"
3040            wrongstats.append(varn + '_' + opers)
3041            break
3042        elif op == 'tmean':
3043            # temporal mean values
3044            vals='time|-1,time,mean,' + ':'.join(varnCFs)
3045            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
3046            if dofile:
3047                try:
3048                    with gen.Capturing() as output:
3049                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3050                except:
3051                    print errmsg
3052                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3053                    for s1out in output: print s1out
3054                    quit(-1)
3055
3056                if db:
3057                    for s1out in output: print s1out
3058
3059                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)
3060
3061                # Keeping the operations
3062                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3063                  "' -f " + prevfile + " -v " + CFvarnp
3064                otrackf.write("\n")
3065                otrackf.write("# " + varn + " " + Fopers + "\n")
3066                otrackf.write(pyins + "\n")
3067
3068            # removing dimension variable-dimension 'time'
3069            varnCFs.remove('time')
3070
3071            varkeep.append('timestats')
3072
3073        elif op == 'tstd':
3074            # temporal standard deviation values
3075            vals='time|-1,time,std,' + ':'.join(varnCFs) + ':' + vdnz
3076            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
3077            if dofile:
3078                try:
3079                    with gen.Capturing() as output:
3080                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3081                except:
3082                    print errmsg
3083                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3084                    for s1out in output: print s1out
3085                    quit(-1)
3086
3087                if db:
3088                    for s1out in output: print s1out
3089
3090                sout = sub.call('mv file_oper_alongdims_std.nc '+fileon, shell=True)
3091
3092                # Keeping the operations
3093                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3094                  "' -f " + prevfile + " -v " + CFvarnp
3095                otrackf.write("\n")
3096                otrackf.write("# " + varn + " " + Fopers + "\n")
3097                otrackf.write(pyins + "\n")
3098
3099            # removing dimension variable-dimension 'time'
3100            varnCFs.remove('time')
3101
3102            varkeep.append('timestats')
3103
3104        elif op == 'tturb':
3105            # temporal turbulent values
3106            vals='time|-1,time,turb,' + ':'.join(varnCFs)
3107            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
3108            if dofile:
3109                try:
3110                    with gen.Capturing() as output:
3111                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3112                except:
3113                    print errmsg
3114                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3115                    for s1out in output: print s1out
3116                    quit(-1)
3117
3118                if db:
3119                    for s1out in output: print s1out
3120
3121                sout = sub.call('mv file_oper_alongdims_turb.nc '+fileon, shell=True)
3122
3123                # Keeping the operations
3124                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3125                  "' -f " + prevfile + " -v " + CFvarnp
3126                otrackf.write("\n")
3127                otrackf.write("# " + varn + " " + Fopers + "\n")
3128                otrackf.write(pyins + "\n")
3129
3130            varkeep.append('timestats')
3131
3132        elif op == 'tvar':
3133            # temporal variance values (Taylor's turbulence)
3134            vals='time|-1,time,var,' + ':'.join(varnCFs) + ':' + vdnz
3135            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
3136            if dofile:
3137                try:
3138                    with gen.Capturing() as output:
3139                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3140                except:
3141                    print errmsg
3142                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3143                    for s1out in output: print s1out
3144                    quit(-1)
3145
3146                if db:
3147                    for s1out in output: print s1out
3148
3149                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)
3150
3151                # Keeping the operations
3152                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3153                  "' -f " + prevfile + " -v " + CFvarnp
3154                otrackf.write("\n")
3155                otrackf.write("# " + varn + " " + Fopers + "\n")
3156                otrackf.write(pyins + "\n")
3157
3158            # removing dimension variable-dimension 'time'
3159            varnCFs.remove('time')
3160
3161            varkeep.append('timestats')
3162
3163        elif op == 'xmean':
3164            # x-axis mean values
3165            vals='lon|-1,lon,mean,' + ':'.join(varnCFs)
3166            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
3167            if dofile:
3168                try:
3169                    with gen.Capturing() as output:
3170                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3171                except:
3172                    print errmsg
3173                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3174                    for s1out in output: print s1out
3175                    quit(-1)
3176
3177                if db:
3178                    for s1out in output: print s1out
3179
3180                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)
3181
3182                # Keeping the operations
3183                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3184                  "' -f " + prevfile + " -v " + CFvarnp
3185                otrackf.write("\n")
3186                otrackf.write("# " + varn + " " + Fopers + "\n")
3187                otrackf.write(pyins + "\n")
3188
3189            # removing dimension variable-dimension 'lon'
3190            varnCFs.remove('lon')
3191
3192            varkeep.append('lonstats')
3193
3194        elif op == 'xvar':
3195            # x-axis variance values
3196            vals='lon|-1,lon,var,' + ':'.join(varnCFs)
3197            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
3198            if dofile:
3199                try:
3200                    with gen.Capturing() as output:
3201                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3202                except:
3203                    print errmsg
3204                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3205                    for s1out in output: print s1out
3206                    quit(-1)
3207
3208                if db:
3209                    for s1out in output: print s1out
3210
3211                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)
3212
3213                # Keeping the operations
3214                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3215                  "' -f " + prevfile + " -v " + CFvarnp
3216                otrackf.write("\n")
3217                otrackf.write("# " + varn + " " + Fopers + "\n")
3218                otrackf.write(pyins + "\n")
3219
3220            # removing dimension variable-dimension 'lon'
3221            varnCFs.remove('lon')
3222
3223            varkeep.append('lonstats')
3224
3225        elif op == 'ymean':
3226            # y-axis mean values
3227            vals='lat|-1,lat,mean,' + ':'.join(varnCFs)
3228            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
3229            if dofile:
3230                try:
3231                    with gen.Capturing() as output:
3232                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
3233                except:
3234                    print errmsg
3235                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
3236                    for s1out in output: print s1out
3237                    quit(-1)
3238
3239                if db:
3240                    for s1out in output: print s1out
3241
3242                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)
3243
3244                # Keeping the operations
3245                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
3246                  "' -f " + prevfile + " -v " + CFvarnp
3247                otrackf.write("\n")
3248                otrackf.write("# " + varn + " " + Fopers + "\n")
3249                otrackf.write(pyins + "\n")
3250
3251            # removing dimension variable-dimension 'lat'
3252            varnCFs.remove('lat')
3253
3254            varkeep.append('latstats')
3255
3256        elif op == 'zsum':
3257            # vertical aggregated values
3258            print "  " + fname + ": kind '" + op + "' not ready !!"
3259            wrongstats.append(varn + '_' + opers)
3260            break
3261        else:
3262            print errmsg
3263            print '  ' + fname + ": operation '" + op + "' not ready !!"
3264            quit(-1)
3265
3266        # End of kind of operation
3267
3268        # Variable name in file (vninF) changed due to operation
3269        #   but only if previous operation does not the same 'statistic'
3270        chvn = gen.dictionary_key_list(Opsurs, op)
3271        if chvn is not None:
3272            oldvninF = vninF
3273            vninF = vninF + chvn
3274            CFvarnp = CFvarnp.replace(oldvninF,vninF)
3275               
3276        if dofile:
3277            if len(varkeep) > 0:
3278                varkeepS = ',' + ','.join(varkeep)
3279            else:
3280                varkeepS = ''
3281
3282            # Adding such CF dimension variables which might not be in the
3283            #   post-operation file
3284            try:
3285                with gen.Capturing() as output: 
3286                    varinfile = ncvar.ivars(fileon)
3287            except:
3288                print errmsg
3289                print 'ivars('+fileon+')'
3290                for s1out in output: print s1out
3291                quit(-1)
3292
3293            for CFvn in varnCFs:
3294                if not gen.searchInlist(varinfile, CFvn):
3295                    if db:
3296                        print '  ' + fname + ": recupering CF variable '" + CFvn + "'"
3297                    try:
3298                        with gen.Capturing() as output:
3299                            ncvar.fvaradd(prevfile+','+CFvn, fileon)
3300                    except:
3301                        print errmsg
3302                        print 'fvaradd(' + prevfile + ',' + CFvn +', ' + fileon + ')'
3303                        for s1out in output: print s1out
3304                        quit(-1)
3305
3306                    if db:
3307                        for s1out in output: print s1out
3308
3309            totalvarkeeps = CFvarnp + ',' + ','.join(varnCFs) + varkeepS
3310            try:
3311                with gen.Capturing() as output:
3312                    oclean = ncvar.cleaning_varsfile(totalvarkeeps,fileon)
3313            except:
3314                print errmsg
3315                print 'cleaning_varsfile('+totalvarkeeps+','+fileon+')'
3316                for s1out in output: print s1out
3317                quit(-1)
3318
3319            if db:
3320                for s1out in output: print s1out
3321
3322    # End of operations
3323    otrackf.close()
3324
3325    return
3326
3327def compute_var_diffs(config, diffallvar, odir, diffscr, debug):
3328    """ Function to compute variable differences
3329      config= experiment configuration
3330      diffallvar= dictonary with the differences to compute
3331      odir= output folder
3332      diffscr= wether it should be done from the scratch
3333    """
3334    fname = 'compute_var_diffs'
3335
3336    # output folder
3337    ofold = config['ofold']
3338
3339    # Home of the python scripts
3340    pyH = config['pyHOME']
3341
3342    #  opsur = operation surnames
3343    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
3344    Opsurs = {}
3345    for opk in opsur:
3346        opn = opk.split('_')[1]
3347        vals = opsur[opk]
3348        Opsurs[opn] = vals.split(':')
3349
3350    # Getting in working dir
3351    os.chdir(odir)
3352
3353    # CF dimensions
3354    CFdims = ['lon', 'lat', 'pres', 'time']
3355
3356    # Trakcing file
3357    trkobjf = open(odir + '/all_diffvar.inf', 'a')
3358
3359    for vnop in diffallvar.keys():
3360        vn = vnop.split('_')[0]
3361        op = vnop.split('_')[1]
3362        opSfinal = op.replace('+','_')
3363
3364        lop = op.split('+')
3365        # `p' header add related to 'pinterp'
3366        if op.find('pinterp') != -1:
3367            print warnmsg
3368            print '  ' + fname + ': there is the vertical interpolation as operation'
3369            print '    variable differences will be computed from that files where '
3370            print '    vertical interpolation was just done'
3371            ipop = lop.index('pinterp')
3372            opS = '_' + '_'.join(lop[0:ipop+1])
3373            doops = lop[ipop+1:]
3374            SgP = 'p'
3375            vninF = varnoper(vn, '+'.join(lop[0:ipop+1]), Opsurs)
3376            print '    operations to perform:', doops
3377        else: 
3378            SgP = ''
3379            opS = ''
3380            doops = lop
3381            vninF = vn
3382
3383        Sdoops = '+'.join(doops)
3384
3385        # Variable name within the file
3386
3387        vals = diffallvar[vnop]
3388        mod1 = vals[0]
3389        mod2 = vals[1]
3390        exp1 = vals[2]
3391        exp2 = vals[3]
3392        headerf1 = vals[4]
3393        headerf2 = vals[5]
3394        modexpdiff = mod2 + '/' + exp2 + ' - ' + mod1 + '/' + exp1
3395        modexpdiffS = mod2 + '-' + exp2 + '_' + mod1 + '-' + exp1
3396
3397        ifile1 = ofold + '/' + mod1 + '/' + exp1 + '/' + vn + '_' + headerf1 + SgP + \
3398          opS + '.nc'
3399        ifile2 = ofold + '/' + mod2 + '/' + exp2 + '/' + vn + '_' + headerf2 + SgP + \
3400          opS + '.nc'
3401        difffilen = odir + '/' + vn + '_diffvar_' + modexpdiffS + opS + '.nc'
3402        difffinalfilen = odir + '/' + vn + '_diffvar_' + modexpdiffS + opS + '_' +   \
3403          Sdoops + '.nc'
3404
3405        if diffscr:
3406            sout = sub.call('rm ' + difffilen + ' >& /dev/null', shell=True)
3407            sout = sub.call('rm ' + difffinalfilen + ' >& /dev/null', shell=True)
3408
3409        if not os.path.isfile(difffinalfilen):
3410            if debug:
3411                print "    Computing variable difference of '" + vn + "' for: '" +   \
3412                  modexpdiff + "' and performing operation '" + Sdoops + "'"
3413                print "    var in file: '" + vninF + "'"
3414
3415            dfiles = 'add|' + ifile2 + '|' + vninF + ',sub|' + ifile1 + '|' + vninF
3416            values = '|'.join(CFdims) + '@' + dfiles
3417            pyins = 'python ' + pyH + "/nc_var.py -o compute_opersvarsfiles -S '" +  \
3418              values + "' -v " + vninF
3419
3420            try:
3421                with gen.Capturing() as output:
3422                    ncvar.compute_opersvarsfiles(values,vninF)
3423            except:
3424                print errmsg
3425                print 'ncvar.compute_opersvarsfiles(' + values + ',' + vninF + ')'
3426                for s1out in output: print s1out
3427                quit(-1)
3428
3429            sout = sub.call('mv opersvarsfiles_'+vninF+'.nc ' + difffilen, shell=True)
3430
3431            # keeping all differences
3432            trkobjf.write('\n')
3433            trkobjf.write("#" + vn + ' ' + op + '\n')
3434            trkobjf.write( pyins + '\n')
3435
3436            # Computing statisitcs
3437            compute_diff_statistics(config, difffilen, odir, vninF, doops, diffscr,  \
3438              debug)
3439
3440    trkobjf.close()
3441
3442    return
3443
3444
3445def draw_diff_plots(config, plots, odir, allvarcomp, kdiff, figscr, debug):
3446    """ Function to draw all plots
3447      config= Configuration of the experiment
3448      plots= dictionary with the plots
3449      odir= output experiment folder
3450      allvarcomp= dictionary with all the variables to compute and their information
3451      kdiff= kind of differences:
3452        'diffop': plots from operation differences
3453        'diffvar': plots from variable differences
3454      figscr= whether figures should be done from the scratch or not
3455
3456    * Plot as
3457      {[kplot]} = [varn1]|[op1]#[varn2]|[op2]#[...[varnN]|[opN]], ...
3458        [kplot] ___
3459          diffmap2Dsfc: 2D map of surface differences values of 1 variable
3460          diffmap2Dz: 2D map of 3D differences values of 1 variable
3461          map2Dsfc: 2D map of surface values of 1 variable
3462          map3D: 2D map of 3D values of 1 variable
3463          shadconthovmsfc: Hovmoeller diagrams of 2 variable at the surface in shadow and the other in contourn
3464          shadcont2Dsfc: 2D map of shadow (1st variable) and countour (2nd variable) [stvar1]#[stvar2]
3465          shadcont2Dzsec: 2D map of vertical section of 2 variables one in shadow and the other in contourn
3466        [varn]
3467          variable
3468        [op]
3469          '+' separated list of operations
3470        in figures with more than 1 variable, use '#' to separate the [varn]|[op]
3471    """ 
3472    fname = 'draw_diff_plots'
3473
3474    os.chdir(odir)
3475
3476    # Dictionary with the operations with surnames for the operated variable
3477    opersurnames = {}
3478    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
3479    for opsr in opsur.keys():
3480        opn = opsr.split('_')[1]
3481        vls = opsur[opsr].split(':')
3482        opersurnames[opn] = vls
3483
3484    # time values
3485    # Units time for the plots
3486    rd = config['CFreftime']
3487    tunits = config['CFunitstime'] + '!since!' + rd[0:4] + '-' + rd[4:6] + '-' +  \
3488      rd[6:8] + '!' + rd[8:10] + ':' + rd[10:12] + ':' + rd[12:14]
3489    # time ticks kind
3490    tkind = config['timekind']
3491    # time ticks format
3492    tfmt = config['timefmt']
3493    # time axis label
3494    tlab = config['timelabel']
3495    timevals = [tunits, tkind, tfmt, tlab]
3496
3497    if kdiff == 'diffop':
3498         specplotkeyn = 'specificdiffopplot'
3499    elif kdiff == 'diffvar':
3500         specplotkeyn = 'specificdiffvarplot'
3501    else:
3502         print errmsg
3503         print '  ' + fname + ": differences kind '" + kdiff + "' not ready !!"
3504         quit(-1)
3505
3506    # Dictionary of plot specificities
3507    plotspecifics = {}
3508    if config.has_key(specplotkeyn):
3509        #   [minval]: minimum value
3510        #   [maxval]: minimum value
3511        #   [colorbar]: name of the colorbar (from matplotlib) to use
3512        #   [cntformat]: format of the contour labels
3513        #   [colorcnt]: color for the countor lines
3514        plotspecs = config[specplotkeyn].split(':')
3515        for pltspc in plotspecs:
3516            pltvls = pltspc.split('|')
3517            vn = pltvls[0]
3518            op = pltvls[1]
3519            fn = pltvls[2]
3520            plotspecifics[fn + '_' + vn + '_' + op] = pltvls[3:]
3521        if debug:
3522            print 'Specific values for plots _______'
3523            gen.printing_dictionary(plotspecifics)
3524
3525    # Kind of figures
3526    kindfigure = config['kindfig']
3527
3528    # Graphical labels and configuration
3529    modgc, expgc = graphical_characteristics(config, debug)
3530
3531    # Map value
3532    mapvalue = config['mapval']
3533
3534    # pythone scripts HOME
3535    pyHOME = config['pyHOME']
3536
3537    # Title-text of operations
3538    opexplained = {}
3539    optits = config['titleoperations'].split(':')
3540    for optit in optits:
3541        opn = optit.split('|')[0]
3542        opt = optit.split('|')[1]
3543        opexplained[opn] = opt
3544    if debug:
3545        print 'Titles for operations  _______'
3546        gen.printing_dictionary(opexplained)
3547
3548    for kplot in plots.keys():
3549        varsplt = plots[kplot]
3550        for varplt in varsplt:
3551            if debug:
3552                print "  printing '" + kplot + "' var ':" + varplt + "'..."
3553            varops = varplt.split('#')
3554
3555            # CF variables in plot
3556            CFvarsplot = []
3557            # Files in plot
3558            filesplot = []
3559            # Variables in plot within the files
3560            varsplot = []
3561            # Dims in figure
3562            dimsplot = []
3563            # pictoric values in figure
3564            pictplot = []
3565            # Name of the figure
3566            figname = ''
3567            # Title of the figure
3568            titfigure = ''
3569
3570            ivp = 0
3571            for varop in varops:
3572                vn = varop.split('|')[0]
3573                op = varop.split('|')[1]
3574
3575                # CF variables in plot
3576                CFvarsplot.append(vn)
3577 
3578                vnopS = vn + '_' + op
3579                if not allvarcomp.has_key(vnopS):
3580                    print errmsg
3581                    print '  ' + fname + ": no file for variable-operation '" +     \
3582                      vnopS + "' !!"
3583                vopvals = allvarcomp[vnopS]
3584                mod1 = vopvals[0]
3585                mod2 = vopvals[1]
3586                mod3 = vopvals[2]
3587                mod4 = vopvals[3]
3588                headf = vopvals[4]
3589
3590                modv = modgc[mod1]
3591                expv = expgc[mod1+'/'+exp1]
3592                modlab1 = modv.label
3593                explab1 = expv.label
3594                modv = modgc[mod2]
3595                expv = expgc[mod2+'/'+exp2]
3596                modlab2 = modv.label
3597                explab2 = expv.label
3598
3599                modexpdiff = explab2 + '-' + explab1
3600                modexpdiffS = mod2 + '-' + exp2 + '_' + mod1 + '-' + exp1
3601
3602                difffilen = odir + '/' + vn + '_' + kdiff + '_' + modexpdiffS +      \
3603                  '_' + op.replace('+','_') + '.nc'
3604
3605                filesplot.append(difffilen)
3606                # Do we have processed the given difference?
3607                if not os.path.isfile(difffilen):
3608                    print warnmsg
3609                    print "  " + fname + ": there is no file for '" + kdiff +        \
3610                      "' difference '" + varop + "' skiping it !!"
3611                    break
3612
3613                # Name of the variable inside the file
3614                vnsur = varnoper(vn, op, opersurnames)
3615                varsplot.append(vnsur)
3616
3617                # Dimensions in file
3618                try:
3619                    with gen.Capturing() as output:
3620                        dims = ncvar.idims(difffilen)
3621                except:
3622                    print errmsg
3623                    print 'ncvar.idims('+difffilen+')'
3624                    for s1out in output: print s1out
3625                    quit(-1)
3626
3627                dimsplot.append(dims)
3628
3629                # pictoric values for the figure
3630                Sfivaop = kplot + '_' + vn + '_' + op
3631                if plotspecifics.has_key(Sfivaop):
3632                    pictvals = plotspecifics[Sfivaop]
3633                else:
3634                    Vvals = gen.variables_values(vn)
3635                    pictvals = [Vvals[2], Vvals[3], Vvals[6], '%g', 'fixsigc','black']
3636
3637                pictplot.append(pictvals)
3638
3639                # Header of the name of the figure
3640                if ivp == 0:
3641                    figname = kplot + '_' + vn + '_' + kdiff + '_' + modexpdiffS +   \
3642                      '_' + op.replace('+','-')
3643                else:
3644                    figname = figname + '-' + vn + '_' + op.replace('+','-')
3645
3646                # Title of the figure:
3647                if mod1 == mod2:
3648                    modexptit = mod1 + '!' + exp2 + '-' + exp1
3649                else:
3650                    modexptit = modexpdiff
3651
3652                ivp = ivp + 1
3653            # End of variable-operation
3654            figname = figname + '.' + kindfigure
3655
3656            # Title of figure
3657            titfigure = create_figure_title(modexptit, kdiff, varops, opexplained)
3658
3659            if len(titfigure) > 80:
3660                print warnmsg
3661                print ' ' + fname + ": figure title '" + titfigure.replace('!', ' ')+\
3662                  "' larger than 80 characters (actually", len(titfigure), ') !!'
3663                print "  simplifying it to '" + modexptit.replace('!',' ') + '!' +   \
3664                  kdiff + '!' + vnS + "'"
3665                titfigure = modexptit.replace(' ','!') + '!' + kdiff + '!' + vnS
3666
3667            draw_plot(kplot, CFvarsplot, filesplot, varsplot, dimsplot, pictplot,    \
3668              figname, titfigure, kindfigure, mapvalue, timevals, expgc, odir,       \
3669              pyHOME, figscr, debug)
3670
3671        # End of variables-operations
3672
3673    # End of kind of plots
3674
3675    return
3676
3677def reproject_modexp(mod,exp,config,CDOops,fscr,fadd,debug):
3678    """ Function to re-project and re-compute files from a given model/experiment
3679      mod= model
3680      exp= experiment
3681      config= configuration
3682      CDOops= dictionary with the CDO operators to use for reproject and for which variable
3683      fscr= re-projected files to compute from the scratch
3684      fadd= re-projected files to be added
3685    """
3686    fname = 'reproject_modexp'
3687
3688    odir = config['ofold'] + '/' + mod + '/' + exp
3689
3690    # NON re-projectable operations
3691    opNOreproj = config['NOreprojops'].split(':')
3692
3693    # Need to pass to analyze all the data?
3694    inffiles = ['reproj_varcompute.inf', 'all_reproj_computevars.inf']
3695    if fscr:
3696        for inff in inffiles:
3697            if debug: print "    removing information file '" + inff + "' ..."
3698            ins = 'rm ' + odir + '/' + inff + ' >& /dev/null'
3699            sout = sub.call(ins, shell=True)
3700
3701        objf = open(odir+'/all_reproj_computevars.inf','w')
3702        objf.write("## Computation of reprojected variables \n")
3703        objf.close()
3704
3705    varcompf = odir + '/reproj_varcompute.inf'
3706    if fadd:
3707        sub.call('rm ' + varcompf +' >& /dev/null', shell=True)
3708
3709    if not os.path.isfile(varcompf):
3710        # Getting original computing information
3711        origvarcompf = odir + '/varcompute.inf'
3712        origfiles,origtestfiles,origallcompvar, Nvar = read_varcomp_file(origvarcompf)
3713             
3714        Svarcompute = ''
3715        allcompvar = dict(origallcompvar)
3716        ivop = 0
3717        for vnop in origallcompvar.keys():
3718            vn = vnop.split('_')[0]
3719            op = vnop.split('_')[1]
3720            values = allcompvar[vnop]
3721            values[0] = 'reproj-' + values[0]
3722            allcompvar[vnop] = values
3723
3724            # Building strings from list values
3725            model = values[1]
3726            diag = values[2]
3727            if model is not None:
3728                if type(model) == type(list([1, 2])):
3729                    Smodel = ':'.join(model)
3730                elif type(model) == type('A'):
3731                    Smodel = model
3732            else:
3733                Smodel = 'None'
3734            if diag is not None:
3735                if type(diag) == type(list([1, 2])):
3736                    Sdiag = ':'.join(diag)
3737                elif type(diag) == type('A'):
3738                    Sdiag = diag
3739            else:
3740                Sdiag = 'None'
3741
3742            Svc = vn + '|' + op + '|' + values[0] + '|' + Smodel + '|' + Sdiag +     \
3743              '|' + values[3]
3744            if ivop == 0:
3745                Svarcompute = Svc
3746            else:
3747                Svarcompute = Svarcompute + ',' + Svc
3748            ivop = ivop + 1
3749
3750        origstrfiles = gen.dictKeysVals_stringList(origfiles)
3751        origstrtestfiles = gen.dictKeysVals_stringList(origtestfiles)
3752
3753        # Outwritting the varcompute to avoid next time (if it is not filescratch!)
3754        objf = open(odir + '/reproj_varcompute.inf', 'w')
3755        objf.write('files: ' + origstrfiles + '\n')
3756        objf.write('testfiles: ' + origstrtestfiles + '\n')
3757        objf.write('varcompute: ' + Svarcompute + '\n')
3758        objf.write('itotv: ' + str(ivop) + '\n')
3759        objf.close() 
3760    else: 
3761        print warnmsg
3762        print '  ' + main + ": getting variables to compute already from file !!"
3763        origfiles, origtestfiles, allcompvar, Nvar = read_varcomp_file(varcompf)
3764
3765    # End of avoiding to repeat all the experiment search
3766
3767    print "    For experiment '" + exp + "' is required to re-project:", Nvar,       \
3768      "variables"
3769
3770    if dbg:
3771        print 'Variables to re-project _______'
3772        gen.printing_dictionary(allcompvar)
3773
3774    # Re-projection variables
3775    print "    Reprojecting variables ..."
3776    os.chdir(odir)
3777
3778    # Keeping track of reprojection
3779    trckobjf = open(odir + '/all_reproj_computevars.inf','a')
3780
3781    for vnop in allcompvar.keys():
3782        vn = vnop.split('_')[0]
3783        op = vnop.split('_')[1]
3784        values = allcompvar[vnop]
3785
3786        # Only are reprojected that files which do contain 'lon' and 'lat',
3787        #   otherways, it is only re-projected the original file
3788        inNO, iopNO, iop = gen.list_coincidences(opNOreproj,op.split('+'))
3789
3790        cdoop = gen.dictionary_key_list(CDOops,vn).split('_')[1]
3791        if cdoop is None:
3792            print errmsg
3793            print '  ' + fname + ": no reprojecting CDO operator found for " +       \
3794              " variable '" + vn + "' !!"
3795            print '    available values _______'
3796            gen.printing_dictionary(CDOops)
3797            quit(-1)
3798       
3799        # Re-projecting also the 'original' (at least with 'pinterp') file for varsdif
3800        if op.find('pinterp') != -1:
3801            ifileorig = vn +'_'+ values[0].replace('reproj-','') + 'p_pinterp.nc'
3802            ofileorig = vn +'_'+ values[0] + 'p_pinterp.nc'
3803            ofilen = vn +'_'+ values[0] + 'p_pinterp.nc'
3804        else:
3805            ifileorig = vn +'_'+ values[0].replace('reproj-','') + '.nc'
3806            ofileorig = vn +'_'+ values[0] + '.nc'
3807
3808        ifileop = vn+'_'+values[0].replace('reproj-','')+'_'+op.replace('+','_')+'.nc'
3809        ofileop = vn + '_' + values[0] + '_' + op.replace('+','_') + '.nc'
3810
3811        if fscr:
3812            sout = sub.call('rm ' + odir+'/'+ofileorig + ' >& /dev/null', shell=True)
3813            sout = sub.call('rm ' + odir+'/'+ofileop + ' >& /dev/null', shell=True)
3814
3815        ifiles = [odir+'/'+ifileorig, odir+'/'+ifileop]
3816        ofiles = [odir+'/'+ofileorig, odir+'/'+ofileop]
3817
3818        for iif in range(2):
3819            ifile = ifiles[iif]
3820            ofile = ofiles[iif]
3821
3822            # Only projecting original file-variable if there are not 'lon' and 'lat'
3823            doreproj = True
3824            if iif == 1 and len(inNO) > 0: 
3825                doreproj = False
3826                print warnmsg
3827                print '  ' + fname + ": operation '" + op + "' can not be " +        \
3828                  "reprojected due to the lack of 'lon' or 'lat' skipping it !!"
3829                print "    NOTE: difference can only be computed as 'variable' " +   \
3830                  "difference"
3831                print "      As operating along the differences from the original" + \
3832                  " values from each model"
3833
3834                # Building up resources to compute re-projected statistics
3835                print '  ' + fname + ': computing from the original reprojected file'
3836                ModI = ModelInf(mod, mod, 'lon', 'lat', 'pres', 'time', 'depth',     \
3837                  'lon', 'lat', 'pres', 'time', 'depth', None)
3838                VarI = VariableInf(vn, values[0], values[2], values[3])
3839                usefs = {values[0]: ofileorig}
3840
3841                compute_statistics(ModI, config, config['ifold'], usefs,             \
3842                  config['ofold']+'/'+mod+'/'+exp,VarI,pinterp_var(op), op, fscr, dbg)
3843
3844            if not os.path.isfile(ofile) and doreproj:
3845                if debug and ofile == ofiles[0]:
3846                    print "      reprojecting '" + vn + "' '" + op + "' into '" +    \
3847                      ofile + "'"
3848
3849                ins = config['cdoHOME'] + '/cdo remap' + cdoop + ',' +               \
3850                  config['RefProjfile'] + ' ' + ifile + ' ' + ofile
3851                try:
3852                    with gen.Capturing() as output:
3853                        sout = sub.call(ins, shell=True)
3854                except:
3855                    print errmsg
3856                    print '  ' + fname + ': $ ' + ins
3857                    for s1out in output: print s1out
3858                    print sout
3859                    quit(-1)
3860                if sout != 0:
3861                    print errmsg
3862                    print '  ' + fname + ': $ ' + ins
3863                    sub.call(ins, shell=True)
3864                    quit(-1)
3865
3866                if debug:
3867                    for s1out in output: print s1out
3868                    print sout
3869
3870                # Keeping track of re-prjection
3871                trckobjf.write('\n')
3872                if iif == 0:
3873                    trckobjf.write('# ' + vn + '\n')
3874                else:
3875                    trckobjf.write('# ' + vn + ' ' + op + '\n')
3876
3877                trckobjf.write(ins + '\n')
3878
3879                # CFing the reprojected file
3880                try:
3881                    with gen.Capturing() as output:
3882                        ncvar.CDO_toCF(ofile)
3883                except:
3884                    print errmsg
3885                    print 'ncvar.CDO_toCF(' + ofile + ')'
3886                    for s1out in output: print s1out
3887                    quit(-1)
3888
3889                if debug:
3890                    for s1out in output: print s1out
3891               
3892    trckobjf.close()
3893
3894    return
3895
3896def allmodexps_listconstruct(config, modexps, odir, debug):
3897    """ Function to create the list of all model-experiments to draw
3898      config= Configuration of the experiment
3899      modexps= list with the models-experiments pairs to use
3900      odir= output differences folder
3901    """
3902    fname='allmodexps_listconstruct'
3903 
3904    allplts = gen.get_specdictionary_HMT(config, H='PLOTALLMODEXP_')
3905
3906    if debug:
3907        if allplts is not None:
3908            print "  'all mod-exps' plots ________"
3909            gen.printing_dictionary(allplts)
3910
3911    # Getting Model-Reference Projection
3912    ModComProj = cnf['RefProj'].split(':')
3913
3914    allmodexp = {}
3915    Nallmodexp = 0
3916    for modexp in modexps:
3917        if debug:
3918            print "    " +  fname + ": '" + modexp + "' ..."
3919        mod = modexp.split('/')[0]
3920        exp = modexp.split('/')[1]
3921
3922        # Getting variable characteristics from each model/experiment
3923        idir = config['ofold'] + '/' + modexp
3924
3925        reproj = not gen.searchInlist(ModComProj,mod)
3926
3927        fvarcomp = 'varcompute.inf'
3928        if reproj: fvc = 'reproj_' + fvarcomp
3929        else: fvc = fvarcomp
3930
3931        files, testfiles, allcompvar, Nvar= read_varcomp_file(idir + '/' + fvc)
3932
3933        if allplts is not None:
3934            # Getting all model-experiments plots characteristics
3935
3936            for allplt in allplts.keys():
3937                pltn = allplt.split('_')[1]
3938                if debug and modexp == modexps[0]:
3939                    print '    ' + fname + ":      plot '" + pltn + "' "
3940                varops = allplts[allplt].split(':')
3941                for varop in varops:
3942                    vn = varop.split('|')[0]
3943                    op = varop.split('|')[1]
3944                    vnop = vn + '_' + op
3945                    if debug and modexp == modexps[0]:
3946                        print '    ' + fname + ":        varop '" + vnop + "' "
3947                    if not allcompvar.has_key(vnop):
3948                        print warnmsg
3949                        print '  ' + fname + ": variable+operation '" + vnop +       \
3950                          "' in '" + modexp + "' was not computed !!"
3951                        print '    skipping it!'
3952                        print '    computed values:', allcompvar.keys()
3953                        break
3954
3955                    vals = allcompvar[vnop]
3956
3957                    headerf = vals[0]
3958                    if modexp == modexps[0]:
3959                        allmodexp[vnop] = [modexp,headerf]
3960                        Nallmodexp = Nallmodexp + 1
3961                    else:
3962                        prevals = allmodexp[vnop]
3963                        prevals = prevals + [modexp,headerf]
3964                        allmodexp[vnop] = prevals
3965
3966        else:
3967            Nallmodexp = 0
3968
3969    if allplts is not None:
3970        Sallmodexps = gen.dictKeysVals_stringList(allmodexp)
3971    else:
3972        Sallmodexps = 'none'
3973
3974    # Outwritting the all model-experiment plots file to avoid next time
3975    objf = open(odir + '/allmodexp.inf', 'w')
3976    objf.write('allmodexps: ' + Sallmodexps + '\n')
3977    objf.write('Nallmodexp: ' + str(Nallmodexp) + '\n')
3978    objf.close()
3979
3980    return allmodexp, Nallmodexp
3981
3982def allmodexps_read(filen):
3983    """ Function to read the file with the information about the all model-experiment plots
3984      filen= file with the information
3985    """
3986    fname = 'allmodexps_read'
3987
3988    if not os.path.isfile(filen):
3989        print errormsg
3990        print '  ' + fname + ": all model-experiment file '" + filen +               \
3991          "' does not exist !!"
3992        quit(-1)
3993
3994    objf = open(filen, 'r')
3995    for line in objf:
3996        if line[0:1] != '#' and len(line) > 1:
3997            values = line.replace('\n','').split(' ')
3998            if values[0] == 'allmodexps:':
3999                allplts = gen.stringList_dictKeysVals(values[1],cV=':')
4000            elif values[0] == 'Nallmodexp:':
4001                Nallmodexp = int(values[1])
4002
4003    objf.close()
4004
4005    return allplts, Nallmodexp
4006
4007def draw_allmodexp_plots(config, allvarcomp, plots, modexp, odir, figscr, debug):
4008    """ Function to draw all model-experiment plots
4009      config= Configuration of the experiment
4010      allvarcomp = dictionary with the file headers for each variable
4011      plots= dictionary with the plots
4012      modexps= list with the names of the model-experiment pairs
4013      odir= output experiment folder
4014      figscr= whether figures should be done from the scratch or not
4015
4016    * Plot as
4017      {[kplot]} = [varn1]|[op1]#[varn2]|[op2]#[...[varnN]|[opN]], ...
4018        [kplot] ___
4019          Nlines: Multiple lines with one line by model-experiment pairs
4020        [varn]
4021          variable
4022        [op]
4023          '+' separated list of operations
4024        in figures with more than 1 variable, use '#' to separate the [varn]|[op]
4025    """ 
4026    fname = 'draw_allmodexp_plots'
4027
4028    os.chdir(odir)
4029
4030    # Folder with the processed files
4031    infold = config['ofold']
4032
4033    # Dictionary with the operations with surnames for the operated variable
4034    opersurnames = {}
4035    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
4036    for opsr in opsur.keys():
4037        opn = opsr.split('_')[1]
4038        vls = opsur[opsr].split(':')
4039        opersurnames[opn] = vls
4040
4041    # time values
4042    # Units time for the plots
4043    rd = config['CFreftime']
4044    tunits = config['CFunitstime'] + '!since!' + rd[0:4] + '-' + rd[4:6] + '-' +  \
4045      rd[6:8] + '!' + rd[8:10] + ':' + rd[10:12] + ':' + rd[12:14]
4046    # time ticks kind
4047    tkind = config['timekind']
4048    # time ticks format
4049    tfmt = config['timefmt']
4050    # time axis label
4051    tlab = config['timelabel']
4052    timevals = [tunits, tkind, tfmt, tlab]
4053
4054    # Dictionary of plot specificities
4055    specplotkeyn = 'specificvarplot'
4056    plotspecifics = {}
4057    if config.has_key(specplotkeyn):
4058        #   [minval]: minimum value
4059        #   [maxval]: minimum value
4060        #   [colorbar]: name of the colorbar (from matplotlib) to use
4061        #   [cntformat]: format of the contour labels
4062        #   [colorcnt]: color for the countor lines
4063        plotspecs = config[specplotkeyn].split(':')
4064        for pltspc in plotspecs:
4065            pltvls = pltspc.split('|')
4066            vn = pltvls[0]
4067            op = pltvls[1]
4068            fn = pltvls[2]
4069            plotspecifics[fn + '_' + vn + '_' + op] = pltvls[3:]
4070        if debug:
4071            print 'Specific values for plots _______'
4072            gen.printing_dictionary(plotspecifics)
4073
4074    # Graphical labels and configuration
4075    modgc, expgc = graphical_characteristics(config, debug)
4076
4077    # Kind of figures
4078    kindfigure = config['kindfig']
4079
4080    # Map value
4081    mapvalue = config['mapval']
4082
4083    # pythone scripts HOME
4084    pyHOME = config['pyHOME']
4085
4086    # Title-text of operations
4087    opexplained = {}
4088    optits = config['titleoperations'].split(':')
4089    for optit in optits:
4090        opn = optit.split('|')[0]
4091        opt = optit.split('|')[1]
4092        opexplained[opn] = opt
4093    if debug:
4094        print 'Titles for operations  _______'
4095        gen.printing_dictionary(opexplained)
4096
4097    for kplot in plots.keys():
4098        varsplt = plots[kplot]
4099        for varplt in varsplt:
4100            if debug:
4101                print "  printing '" + kplot + "' var ':" + varplt + "'..."
4102            varops = varplt.split('#')
4103
4104            # CF variables in plot
4105            CFvarsplot = []
4106            # Files in plot
4107            filesplot = []
4108            # Variables in plot within the files
4109            varsplot = []
4110            # Dims in figure
4111            dimsplot = []
4112            # pictoric values in figure
4113            pictplot = []
4114            # Name of the figure
4115            figname = ''
4116            # Title of the figure
4117            titfigure = ''
4118
4119            ivp = 0
4120            for varop in varops:
4121                vn = varop.split('|')[0]
4122                op = varop.split('|')[1]
4123
4124                # CF variables in plot
4125                CFvarsplot.append(vn)
4126 
4127                if op.find('pinterp') != -1: gP = 'p'
4128                else: gP = ''
4129
4130                vnopS = vn + '_' + op                   
4131                if not allvarcomp.has_key(vnopS):
4132                    print errmsg
4133                    print '  ' + fname + ": no file for variable-operation '" +      \
4134                      vnopS + "' !!"
4135                modexpvals = allvarcomp[vnopS]
4136#                modexpvalsS = allvarcomp[vnopS]
4137#                modexpvals = modexpvalsS[0].split('|')
4138                Nmodexps = len(modexpvals)
4139
4140                for imodexp in range(0,Nmodexps,2):
4141                    modexp = modexpvals[imodexp]
4142                    headf = modexpvals[imodexp+1]
4143                    if debug:
4144                        print '  ' + fname + "    model/experiment: '" +  modexp + "'"
4145
4146                    filen = infold + '/' + modexp + '/' + vn +'_'+ headf+gP + '_' +  \
4147                      op.replace('+','_') + '.nc'
4148
4149                    filesplot.append(filen)
4150                    # Do we have processed the given variable?
4151                    if not os.path.isfile(filen):
4152                        print warnmsg
4153                        print "  " + fname + ": there is no file for variable '" +   \
4154                          varop + "' skiping it !!"
4155                        break
4156
4157                    if imodexp == 0:
4158                        # Name of the variable inside the file
4159                        vnsur = varnoper(vn, op, opersurnames)
4160                        varsplot.append(vnsur)
4161
4162                    # Dimensions in file
4163                    try:
4164                        with gen.Capturing() as output:
4165                            dims = ncvar.idims(filen)
4166                    except:
4167                        print errmsg
4168                        print 'ncvar.idims('+filen+')'
4169                        for s1out in output: print s1out
4170                        quit(-1)
4171
4172                    if imodexp == 0:
4173                        # Dimensions in plot
4174                        dimsplot.append(dims)
4175
4176                        # Header of the name of the figure
4177                        if ivp == 0:
4178                            figname = kplot+'_'+vn+'_'+headf+'_'+ op.replace('+','-')
4179                        else:
4180                            figname = figname+'-'+vn+'_'+headf+'_'+op.replace('+','-')
4181
4182                    # pictoric values for the figure
4183                    Sfivaop = kplot + '_' + vn + '_' + op
4184                    if plotspecifics.has_key(Sfivaop):
4185                        pictvals = plotspecifics[Sfivaop]
4186                    else:
4187                        Vvals = gen.variables_values(vn)
4188                        pictvals = [Vvals[2],Vvals[3],Vvals[6], '%g', 'fixc', 'black']
4189
4190                    pictplot.append(pictvals)
4191
4192                ivp = ivp + 1
4193            # End of variable-operation
4194            figname = figname + '.' + kindfigure
4195
4196            # Title of figure
4197            titfigure = 'all!model-experiments!' + vn + optitle(op,opexplained)
4198
4199            draw_plot(kplot, CFvarsplot, filesplot, varsplot, dimsplot, pictplot,    \
4200              figname, titfigure, kindfigure, mapvalue, timevals, expgc, odir,       \
4201              pyHOME, figscr, debug)
4202
4203        # End of variables-operations
4204
4205    # End of kind of plots
4206
4207    return
4208
4209# Files with information about the configuration of the script
4210inffiles = ['varcompute.inf', 'all_computevars.inf', 'all_statsvars.inf']
4211
4212#######    #######
4213## MAIN
4214    #######
4215if not os.path.isfile('model_graphics.dat'):
4216    print errmsg
4217    print main + ": configuration file 'model_graphics.dat' does not exist!!"
4218    quit(-1)
4219
4220# Getting configuration from external ASCII file 'model_graphics.dat'
4221cnf = gen.get_configuration('model_graphics.dat', False)
4222
4223verify_configuration(cnf, gen.Str_Bool(cnf['debug']))
4224
4225# scratches
4226scratch, filescratch, figscratch, diffscratch, figdiffscratch, moddiffscratch,       \
4227  figmoddiffscratch, figallmodexpscratch, addfiles, addfigures, adddiffs,            \
4228  adddifffigures, addmoddiffs, addmoddifffigures, addallmodexpfigures, dbg =         \
4229  scratches(cnf)
4230
4231if dbg:
4232    print 'Experiment configuration scratches _______'
4233    print '  scratch:', scratch, 'filescratch:', filescratch,'figscratch:', figscratch
4234    print '  diffscratch:', diffscratch, 'figdiffscratch:', figdiffscratch
4235    print '  moddiffscratch:', moddiffscratch, 'figmoddiffscratch:', figmoddiffscratch
4236    print 'Experiment configuration adding _______'
4237    print '  addfiles:', addfiles, 'addfigures:', addfigures
4238    print '  adddiffs:', adddiffs, 'adddifffigures:', adddifffigures
4239    print '  addmoddiffs:', addmoddiffs, 'addmoddifffigures:', addmoddifffigures
4240
4241# Getting models
4242mods = cnf['models'].split(':')
4243
4244# Getting graphical characeristics
4245modGC, expGC = graphical_characteristics(cnf, dbg)
4246
4247# Models loop
4248##
4249
4250# dictionary with the experiments of each model
4251modexps = {}
4252
4253for mod in mods:
4254    print mod
4255    if scratch:
4256        if cnf['ofold'] != cnf['ifold']:
4257            sout = sub.call('rm -rf '+cnf['ofold']+'/'+mod+' > /dev/null', shell=True)
4258        else:
4259            print warnmsg
4260            print '  ' + main + ": input '" + cnf['ifold'] + "' and output '" +      \
4261              cnf['ofold'] + "' are the same folder !!"
4262            print "    Keeping output folder although it has 'scratch' start"
4263
4264    # Get experiments and headers of model
4265    exps, fheaders = exp_headers(mod,cnf)
4266
4267    modexps[mod] = exps
4268    # Characteristics of the model
4269    Modinf = ncvar.model_characteristics(mod,'None','False')
4270    dnx = Modinf.dimxn
4271    dny = Modinf.dimyn
4272    dnz = Modinf.dimzn
4273    dnt = Modinf.dimtn
4274    dns = Modinf.dimsn
4275    vdnx = Modinf.vardxn
4276    vdny = Modinf.vardyn
4277    vdnz = Modinf.vardzn
4278    vdnt = Modinf.vardtn
4279    vdns = Modinf.vardsn
4280
4281    modgraphv = modGC[mod]
4282    Modinf = ModelInf(mod, mod, dnx, dny, dnz, dnt, dns, vdnx, vdny, vdnz, vdnt,     \
4283      vdns, modgraphv.tmod)
4284
4285    if dbg:
4286        print '  model characteristics _______'
4287        gen.printing_class(Modinf)
4288
4289    moddims = dnx + ',' + dny + ',' + dnz + ',' + dnt
4290    modvdims = vdnx + ',' + vdny + ',' + vdnz + ',' + vdnt
4291
4292# Experiments loop
4293##
4294    for exp in exps:
4295        print '  ' + exp + '...'
4296
4297        # input folder
4298        iwdir = cnf['ifold'] + '/' + mod + '/' + exp
4299
4300        # Does input folder exist?
4301        if not os.path.isdir(iwdir):
4302            if filescratch*addfiles:
4303                print errmsg
4304                print "  " + main + ": folder '" + iwdir + "' does not exist !!"
4305                quit(-1)
4306            else:
4307                print warnmsg
4308                print "  " + main + ": folder '" + iwdir + "' does not exist !!"
4309                print "    but start of files from scratch 'filescratch':",          \
4310                  filescratch, "and there is not the add of any file 'addfiles':",   \
4311                  addfiles
4312                print '    so, keep moving'
4313
4314        owdir = cnf['ofold'] + '/' + mod + '/' + exp
4315        sout = sub.call('mkdir -p ' + owdir, shell=True)
4316
4317        # Need to pass to analyze all the data?
4318        if filescratch:
4319            for inff in inffiles:
4320                if dbg: print "    removing information file '" + inff + "' ..."
4321                ins = 'rm ' + owdir + '/' + inff + ' >& /dev/null'
4322                sout = sub.call(ins, shell=True)
4323
4324            objf = open(owdir+'/all_computevars.inf','w')
4325            objf.write("## Computation of variables \n")
4326            objf.close()
4327            objf = open(owdir+'/all_statsvars.inf','w')
4328            objf.write("## Computation of statistics \n")
4329            objf.close()
4330
4331        varcompf = owdir + '/varcompute.inf'
4332        if addfiles:
4333            if dbg:
4334                print '  ' + main + ": adding variables to copmute removing " +      \
4335                  "file '" + varcompf + "' ..."
4336            sub.call('rm ' + varcompf +' >& /dev/null', shell=True)
4337
4338        if not os.path.isfile(varcompf):
4339            # Does input folder has header files?
4340            ih=1
4341            # Dictionary with the list of files for each headers
4342            files = {}
4343            # Dictionary with a file fromor each headers
4344            testfiles = {}
4345
4346            for fh in fheaders:
4347                if filescratch:
4348                    ins = 'rm '+ owdir+'/*_' + fh + '*.nc >& /dev/null'
4349                    sout = sub.call(ins, shell=True)
4350                files1h = gen.files_folder(iwdir,fh)
4351                if len(files1h) < 1:
4352                    print errmsg
4353                    print '  ' + main + ": folder '" + iwdir + "' does not " +   \
4354                      "contain files '" + fh + "*' !!"
4355                    quit(-1)
4356                files[fh] = files1h
4357                testfiles[fh] = files1h[0]
4358
4359            if dbg:
4360                print '  Dictionary of files _______'
4361                gen.printing_dictionary(files)
4362
4363                print '  Dictionary of test-files _______'
4364                gen.printing_dictionary(testfiles)
4365         
4366            allcompvar, Nvar = compvars_listconstruct(cnf, Modinf, files, testfiles, \
4367              iwdir, owdir, dbg)
4368        else: 
4369            print warnmsg
4370            print '  ' + main + ": getting variables to compute already from file !!"
4371            files, testfiles, allcompvar, Nvar = read_varcomp_file(varcompf)
4372
4373        # End of avoiding to repeat all the experiment search
4374
4375        print "    For experiment '" + exp + "' is required to compute:", Nvar,      \
4376          "variables"
4377
4378        if dbg:
4379            print 'Variables to compute _______'
4380            gen.printing_dictionary(allcompvar)
4381
4382        # Computing variables
4383        print "    Computing variables ..."
4384        compute_vars(cnf, Modinf, iwdir, owdir, files, allcompvar, filescratch, dbg)
4385
4386        # Figures of variables direct from files
4387        print "    Plotting direct figures ..."
4388        dirfigf = owdir + '/directplotsdraw.inf'
4389        if figscratch:
4390            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4391
4392            objf = open(owdir+'/all_figures.inf','w')
4393            objf.write("## Drawing of all figures \n")
4394            objf.close()
4395
4396        if addfigures:
4397            if dbg:
4398                print '  ' + main + ": adding direct figures removing file '" +      \
4399                  dirfigf + "' ..."
4400            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4401
4402        if not os.path.isfile(dirfigf):
4403            listplots, Nplt = plots_listconstruct(cnf, 'DIRPLT', dirfigf, owdir, dbg)
4404        else:
4405            print warnmsg
4406            print '  ' + main + ": getting plots to draw already from file !!"
4407            listplots, Nplt = read_plot_file(dirfigf)
4408
4409        # End of avoiding to repeat all the plots search
4410
4411        print "    For experiment '" + exp + "' is required to plot:", Nplt, "plots"
4412
4413        if dbg:
4414            print 'Plots to draw _______'
4415            gen.printing_dictionary(listplots)
4416
4417        draw_plots(cnf, listplots, mod, exp, owdir, allcompvar, figscratch, dbg)
4418
4419    # end of experiments loop
4420
4421### ## #
4422# Experiment differences
4423## # ## #
4424
4425# Experiments loop
4426##
4427    print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
4428    print "  ** '" + mod + "': Inter experiments differences  "
4429    print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
4430    print '    experiments:', exps
4431    # There are two kind of differences:
4432    #   DIFFop: differences between operations of each given variable.
4433    #     Differences are computed directly from the last stage of the operation
4434    #   DIFFvar: operations of differences between variables
4435    #     First are computed the differences from the initial variable file
4436    #     and then operations are made
4437    # NOTE: remember that: meanvar2 - meanvar1 = mean(var2 - var1)
4438    difffiles = ['diffop.inf', 'diffvar.inf']
4439
4440    Nexps = len(exps)
4441    for iexp1 in range(0,Nexps-1):
4442        exp1 = exps[iexp1]
4443        for iexp2 in range(iexp1+1,Nexps):
4444            exp2 = exps[iexp2]
4445            Sexps = exp2 + '-' + exp1
4446            print '  ' + Sexps + '...'
4447            owdir = cnf['ofold'] + '/' + mod + '/' + exp2 + '-' + exp1
4448            sout = sub.call('mkdir -p ' + owdir, shell=True)
4449
4450            # Removing files with list of differences if should be started from scratch or add differences
4451            for fdiff in difffiles:
4452                difff = owdir + '/' + fdiff
4453                if diffscratch:
4454                    sub.call('rm -rf' + owdir +' >& /dev/null', shell=True)
4455                    sub.call('rm ' + difff +' >& /dev/null', shell=True)
4456                    objf = open(owdir+'/all_' + fdiff,'w')
4457                    objf.write("## Computation and drawing of differences " +        \
4458                      "between '" + exp2 + "'-'" + exp1 + "'\n")
4459                    objf.close()
4460                    difff = owdir + '/all_vardiffstatistics.inf'
4461                    sub.call('rm ' + difff +' >& /dev/null', shell=True)
4462                    objf = open(owdir+'/all_' + fdiff,'w')
4463                    objf.write("## Computation of all differences statistics " +     \
4464                      "between '" + exp2 + "'-'" + exp1 + "'\n")
4465                    objf.close()
4466
4467                if adddiffs:
4468                    if dbg:
4469                        print '  ' + main + ": adding differences removing file '" + \
4470                          difff + "' ..."
4471                    sub.call('rm ' + difff +' >& /dev/null', shell=True)
4472
4473            # op and var differences
4474            diffvarcompf = owdir + '/' + difffiles[0]
4475            if not os.path.isfile(diffvarcompf):
4476                alldiffop, alldiffvar,Nopdiffs,Nvardiffs=diffvarop_listconstruct(cnf,\
4477                  [mod, mod], [exp1, exp2], [False, False], owdir, dbg)
4478            else: 
4479                print warnmsg
4480                print '  ' + main + ": getting 'op' and 'var'differences to " +      \
4481                  "compute already from file !!"
4482                alldiffop, Nopdiffops = read_diff_file(diffvarcompf)
4483                alldiffvar, Nopdiffvars = read_diff_file(owdir+'/'+difffiles[1])
4484
4485        # End of avoiding to repeat all the experiment search
4486
4487        print "    For experiments '"+exp2+"'-'"+exp1+"' is required to compute:",   \
4488          Nvar, "differences"
4489
4490        if dbg:
4491            print 'Differences to compute _______'
4492            gen.printing_dictionary(alldiffop)
4493            gen.printing_dictionary(alldiffvar)
4494
4495        # Computing differences
4496        ##
4497        print "    Computing operation differences ..."
4498        compute_op_diffs(cnf, alldiffop, owdir, diffscratch, dbg)
4499        print "    Computing variable differences ..."
4500        compute_var_diffs(cnf, alldiffvar, owdir, diffscratch, dbg)
4501
4502        # Plotting operation differences
4503        ##
4504        print "  " + main + ": Plotting operation differences' figures ..."
4505        dirfigf = owdir + '/diffopplotsdraw.inf'
4506        if figscratch:
4507            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4508
4509            objf = open(owdir+'/all_diffopfigures.inf','w')
4510            objf.write("## Drawing of all operation difference figures \n")
4511            objf.close()
4512
4513        if adddifffigures:
4514            if dbg:
4515                print '  ' + main + ": adding differences' figures removing file '" +\
4516                  dirfigf + "' ..."
4517
4518            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4519
4520        if not os.path.isfile(dirfigf):
4521            listplots, Nplt = plots_listconstruct(cnf, 'PLOTDIFFOP', dirfigf, owdir, \
4522              dbg)
4523        else:
4524            print warnmsg
4525            print '  ' + main + ": getting plots to draw already from file !!"
4526            listplots, Nplt = read_plot_file(dirfigf)
4527
4528        # End of avoiding to repeat all the plots search
4529
4530        print "  For experiment 'operation' differences '" + Sexps + "' is " +       \
4531          "required to plot:" , Nplt, "plots"
4532
4533        if dbg:
4534            print '    Plots to draw _______'
4535            gen.printing_dictionary(listplots)
4536
4537        draw_diff_plots(cnf, listplots, owdir, alldiffop, 'diffop', figdiffscratch,  \
4538          dbg)
4539
4540        # Plotting variable differences
4541        ##
4542        print "  " + main + ": Plotting variable differences' figures ..."
4543        dirfigf = owdir + '/diffvarplotsdraw.inf'
4544        if figscratch:
4545            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4546
4547            objf = open(owdir+'/all_diffvarfigures.inf','w')
4548            objf.write("## Drawing of all variables difference figures \n")
4549            objf.close()
4550
4551        if adddifffigures:
4552            if dbg:
4553                print '  ' + main + ": adding differences' figures removing file '" +\
4554                  dirfigf + "' ..."
4555            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4556
4557        if not os.path.isfile(dirfigf):
4558            listplots, Nplt = plots_listconstruct(cnf, 'PLOTDIFFVAR', dirfigf, owdir,\
4559              dbg)
4560        else:
4561            print warnmsg
4562            print '  ' + main + ": getting plots to draw already from file !!"
4563            listplots, Nplt = read_plot_file(dirfigf)
4564
4565        # End of avoiding to repeat all the plots search
4566
4567        print "  For experiment 'variables' differences '" + Sexps + "' is " +       \
4568          "required to plot:" , Nplt, "plots"
4569
4570        if dbg:
4571            print '    Plots to draw _______'
4572            gen.printing_dictionary(listplots)
4573
4574        draw_diff_plots(cnf, listplots, owdir, alldiffvar, 'diffvar', figdiffscratch,\
4575          dbg)
4576
4577# end of mods loop
4578
4579### ## #
4580# Model differences
4581### ## #
4582
4583# Models already at the reference projection
4584ModComProj = cnf['RefProj'].split(':')
4585if dbg:
4586    print 'Models to use as reference for the projection:', ModComProj
4587
4588# CDO operators to use for the remapping of each variable
4589CDOremapvar0 = gen.get_specdictionary_HMT(cnf, H='CDOremapvar_')
4590CDOremapvar = {}
4591for Skey in CDOremapvar0.keys():
4592    CDOremapvar[Skey] = gen.str_list(CDOremapvar0[Skey],':')
4593
4594if dbg:
4595    print 'Remaping operators to use _______'
4596    gen.printing_dictionary(CDOremapvar)
4597
4598# Getting common base wights for grid point remapping with CDO
4599frefgridn = 'RefProj.nc'
4600for mod in mods:
4601    frefgriddes = cnf['ofold'] + '/' + mod + '/' + frefgridn
4602    if moddiffscratch: 
4603        sout = sub.call('rm ' + frefgriddes + ' >& /dev/null', shell=True)
4604    if gen.searchInlist(ModComProj,mod) and not os.path.isfile(frefgriddes):
4605        # Looking for a file to use (all files are CF compilant, but we need to
4606        #   be sure that at least they have 'lon', 'lat' variables)
4607        exps = modexps[mod]
4608        exp = exps[0]
4609        cdofdir = cnf['ofold'] + '/' + mod + '/' + exp
4610        cdofile = None
4611        for cdoop in CDOremapvar.keys():
4612            cdovars = CDOremapvar[cdoop]
4613            for cdovar in cdovars:
4614                cdofs= gen.files_folder_HMT(folder=cdofdir,head=cdovar+'_',tail='.nc')
4615                for cdof in cdofs:
4616                    # Getting direct files as [var]_[fhead].nc
4617                    if len(cdof.split('_')) == 2:
4618                        cdofile = cdofdir + '/' + cdof
4619                        cdovn = cdovar
4620                        break
4621        if cdofile is None:
4622            print errmsg
4623            print '  ' + main + ": no proper file to get projection information " +  \
4624              "from '" + cdofdir + "' has been found !!"
4625            quit(-1)
4626
4627        print "  Creation of reference projection information file from '" + mod+ "'"
4628        if dbg:
4629            print "    using file: '" + cdofile
4630
4631        # Getting only 1 time-step for disk space issues
4632        try:
4633             with gen.Capturing() as output:
4634                 ncvar.DataSetSection('time,0,1,1', cdofile)
4635        except:
4636            print errmsg
4637            print 'ncvar.DataSetSection(time,0,1,1, ' + cdofile + ')'
4638            for s1out in output: print s1out
4639            quit(-1)
4640
4641        of = cdofile.split('.')[0] + '_time_B0-E1-I1.nc'
4642
4643        # Cleaning variables
4644        try:
4645             with gen.Capturing() as output:
4646                 ncvar.selvar('lon@lon,lat@lat,time@time', of, cdovn)
4647        except:
4648            print errmsg
4649            print 'ncvar.selvar(lon@lon,lat@lat, ' + of + ', ' + cdovn + ')'
4650            for s1out in output: print s1out
4651            quit(-1)
4652        sout = sub.call('mv selvar_new.nc ' + frefgriddes, shell=True)
4653        sout = sub.call('rm ' + of + ' >& /dev/null', shell=True)
4654
4655        cnf['RefProjfile'] = frefgriddes
4656
4657        if dbg:
4658            for s1out in output: print s1out
4659        break
4660
4661    elif gen.searchInlist(ModComProj,mod) and os.path.isfile(frefgriddes):
4662        # Including reprojecting reference file in configuration dictionary
4663        cnf['RefProjfile'] = frefgriddes
4664
4665Nmods = len(mods)
4666for imod1 in range(0,Nmods-1):
4667    mod1 = mods[imod1]
4668    exps1 = modexps[mod1]
4669
4670    for imod2 in range(imod1+1,Nmods):
4671        mod2 = mods[imod2]
4672        exps2 = modexps[mod2]
4673
4674        Smods = mod2 + '-' + mod1
4675        print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
4676        print "  ** '" + Smods + "': Inter models differences  "
4677        print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
4678        print '    experiments mod1:', exps1
4679        print '    experiments mod2:', exps2
4680
4681# Experiments loop
4682##
4683        difffiles = ['diffop.inf', 'diffvar.inf']
4684
4685        Nexps1 = len(exps1)
4686        Nexps2 = len(exps2)
4687        for iexp1 in range(0,Nexps1):
4688            exp1 = exps1[iexp1]
4689            for iexp2 in range(0,Nexps2):
4690                exp2 = exps2[iexp2]
4691                Sexps = exp2 + '-' + exp1
4692                print '    ' + Sexps + '...'
4693                owdir = cnf['ofold'] + '/' + Smods + '/' + Sexps
4694                sout = sub.call('mkdir -p ' + owdir, shell=True)
4695
4696                # Getting the right projection in order to perform differences
4697                difmods = [mod1, mod2]
4698                difexps = [exp1, exp2]
4699                difreproj = [False, False]
4700               
4701                for imod in range(2):
4702                    mod = difmods[imod]
4703                    exp = difexps[imod]
4704                    if not gen.searchInlist(ModComProj,mod):
4705                        print  "  Projecting files of model '" + mod + "' exp: '" +  \
4706                          exp + "' with file '" + cnf['RefProjfile'] + "'"
4707                        reproject_modexp(mod, exp, cnf, CDOremapvar, moddiffscratch, \
4708                          addmoddiffs, dbg)
4709                        difreproj[imod] = True
4710
4711                # Removing files with list of differences if should be started from
4712                #   scratch or add differences
4713                for fdiff in difffiles:
4714                    difff = owdir + '/' + fdiff
4715                    if moddiffscratch:
4716                        sub.call('rm -rf' + owdir +' >& /dev/null', shell=True)
4717                        sub.call('rm ' + difff +' >& /dev/null', shell=True)
4718                        objf = open(owdir+'/all_' + fdiff,'w')
4719                        objf.write("## Computation and drawing of differences " +    \
4720                          "between '" + mod2+'@'+exp2 + "'-'" + mod1+'@'+exp1 + "'\n")
4721                        objf.close()
4722                        difff = owdir + '/all_vardiffstatistics.inf'
4723                        sub.call('rm ' + difff +' >& /dev/null', shell=True)
4724                        objf = open(owdir+'/all_' + fdiff,'w')
4725                        objf.write("## Computation of all differences statistics " + \
4726                          "between '" + mod2+'@'+exp2 + "'-'" + mod1+'@'+exp1 + "'\n")
4727                        objf.close()
4728
4729                    if addmoddiffs:
4730                        if dbg:
4731                            print '  ' + main + ": adding model differences " +      \
4732                              "removing file '" + difff + "' ..."
4733                        sub.call('rm ' + difff +' >& /dev/null', shell=True)
4734
4735                # op and var differences
4736                diffvarcompf = owdir + '/' + difffiles[0]
4737                if not os.path.isfile(owdir+'/'+difffiles[0]) or                     \
4738                  not os.path.isfile(owdir+'/'+difffiles[1]):
4739                    alldiffop, alldiffvar, Nopdiffs, Nvardiffs =                     \
4740                      diffvarop_listconstruct(cnf, difmods, difexps, difreproj,owdir,\
4741                      dbg)
4742                else: 
4743                    print warnmsg
4744                    print '  ' + main + ": getting 'op' and 'var'differences to " +  \
4745                      "compute already from file !!"
4746                    alldiffop, Nopdiffops = read_diff_file(diffvarcompf)
4747                    alldiffvar, Nopdiffvars = read_diff_file(owdir+'/'+difffiles[1])
4748
4749                # End of avoiding to repeat all the experiment search
4750
4751                print "    For experiments '" + mod2+'@'+exp2 + "'-'" + mod+'@'+exp1+\
4752                  "' is required to compute:", Nvar, "differences"
4753
4754                if dbg:
4755                    print '    op differences to compute _______'
4756                    gen.printing_dictionary(alldiffop)
4757                    print '    var differences to compute _______'
4758                    gen.printing_dictionary(alldiffvar)
4759
4760# Computing differences
4761##
4762                print "    Computing operation differences ..."
4763                compute_op_diffs(cnf, alldiffop, owdir, moddiffscratch, dbg)
4764                print "    Computing variable differences ..."
4765                compute_var_diffs(cnf, alldiffvar, owdir, moddiffscratch, dbg)
4766
4767# Plotting operation differences
4768##
4769                print "  " + main + ": Plotting operation differences' figures ..."
4770                dirfigf = owdir + '/diffopplotsdraw.inf'
4771                if figmoddiffscratch:
4772                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4773
4774                objf = open(owdir+'/all_diffopfigures.inf','w')
4775                objf.write("## Drawing of all operation difference figures \n")
4776                objf.close()
4777
4778                if addmoddifffigures:
4779                    if dbg:
4780                        print '  ' + main + ": adding model differences' figures " + \
4781                          "removing file '" + dirfigf + "' ..."
4782                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4783
4784                if not os.path.isfile(dirfigf):
4785                    listplots, Nplt = plots_listconstruct(cnf, 'PLOTDIFFOP', dirfigf,\
4786                      owdir, dbg)
4787                else:
4788                    print warnmsg
4789                    print '  ' + main + ": getting plots to draw already from file !!"
4790                    listplots, Nplt = read_plot_file(dirfigf)
4791
4792                # End of avoiding to repeat all the plots search
4793
4794                print "  For experiment 'operation' differences '" + Sexps + "' is "+\
4795                  "required to plot:" , Nplt, "plots"
4796
4797                if dbg:
4798                    print 'Plots to draw _______'
4799                    gen.printing_dictionary(listplots)
4800
4801                draw_diff_plots(cnf, listplots, owdir, alldiffop, 'diffop',          \
4802                  figmoddiffscratch, dbg)
4803
4804                # Plotting variable differences
4805                ##
4806                print "  " + main + ": Plotting variable differences' figures ..."
4807                dirfigf = owdir + '/diffvarplotsdraw.inf'
4808                if figscratch:
4809                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4810
4811                    objf = open(owdir+'/all_diffvarfigures.inf','w')
4812                    objf.write("## Drawing of all variables difference figures \n")
4813                    objf.close()
4814       
4815                if adddifffigures:
4816                    if dbg:
4817                        print '  ' + main + ": adding differences's figures " +      \
4818                          "removing file '" + dirfigf + "' ..."
4819                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4820
4821                if not os.path.isfile(dirfigf):
4822                    listplots, Nplt = plots_listconstruct(cnf,'PLOTDIFFVAR',dirfigf, \
4823                      owdir, dbg)
4824                else:
4825                    print warnmsg
4826                    print '  ' + main + ": getting plots to draw already from file !!"
4827                    listplots, Nplt = read_plot_file(dirfigf)
4828
4829                # End of avoiding to repeat all the plots search
4830
4831                print "  For experiment 'variables' differences '" + Sexps+ "' is "+ \
4832                  "required to plot:" , Nplt, "plots"
4833
4834                if dbg:
4835                    print '    Plots to draw _______'
4836                    gen.printing_dictionary(listplots)
4837
4838                draw_diff_plots(cnf, listplots, owdir, alldiffvar, 'diffvar',        \
4839                  figdiffscratch, dbg)
4840
4841            # end of exp2 loop
4842        # end of exp1 loop
4843    # end of mod2 loop
4844# end of mod1 loop
4845
4846### ## #
4847# All model/exps linear plots
4848### ## #
4849
4850allmodexps = []
4851for mod in mods:
4852    for exp in modexps[mod]:
4853        allmodexps.append(mod + '/' + exp)
4854
4855print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
4856print '  ** Multi models, multi experiments plots  '
4857print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
4858print '    Model/experiments:', allmodexps
4859
4860owdir = cnf['ofold'] + '/allmodexps'
4861sout = sub.call('mkdir -p ' + owdir, shell=True)
4862
4863# Plotting all model-experiments lines
4864##
4865print "  " + main + ": Plotting all model-experiments variable figures ..."
4866allf = owdir + '/allmodexp.inf'
4867dirfigf = owdir + '/allmodexpfigures.inf'
4868allfigf = owdir+'/all_figures.inf'
4869if figallmodexpscratch:
4870    sout = sub.call('rm ' + allf + ' >& /dev/null', shell=True)
4871    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4872    sout = sub.call('rm ' + allfigf + ' >& /dev/null', shell=True)
4873
4874    objf = open(allfigf,'w')
4875    objf.write("## Drawing of all variables figures for all model-experiments\n")
4876    objf.close()
4877       
4878if addallmodexpfigures:
4879    if dbg:
4880        print '  ' + main + ": adding model-experiment differences removing " +      \
4881          " file '" + allfigf + "' ..."
4882        print '  ' + main + ": adding model-experiment differences' figures " +      \
4883          "removing file '" + dirfigf + "' ..."
4884    sout = sub.call('rm ' + allfigf + ' >& /dev/null', shell=True)
4885    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
4886
4887if not os.path.isfile(dirfigf):
4888    listplots, Nplt = plots_listconstruct(cnf, 'PLOTALLMODEXP', dirfigf, owdir, dbg)
4889else:
4890    print warnmsg
4891    print '  ' + main + ": getting plots to draw already from file !!"
4892    listplots, Nplt = read_plot_file(dirfigf)
4893
4894# End of avoiding to repeat all the plots search
4895
4896print "  Is required to plot:", Nplt, " all model-experiment plots"
4897
4898if dbg:
4899    print '    Plots to draw _______'
4900    gen.printing_dictionary(listplots)
4901
4902if not os.path.isfile(allf):
4903    # Constructing the information of modexp / variable
4904    allmodexpvar, Nallmodexpvar = allmodexps_listconstruct(cnf, allmodexps, owdir, dbg)
4905else:
4906    allmodexpvar, Nallmodexpvar = allmodexps_read(allf)
4907
4908draw_allmodexp_plots(cnf, allmodexpvar, listplots, modexps, owdir, figallmodexpscratch, dbg)
4909
4910print main + ': all files and figures have been properly done !!'
Note: See TracBrowser for help on using the repository browser.