source: trunk/UTIL/PYTHON/planetoplot_v2/ppclass.py @ 934

Last change on this file since 934 was 934, checked in by aslmd, 12 years ago

UTIL PYTHON planetoplot_v2. thanks to Thomas Navarro added the possibility to change time axis through attribute self.changetime; for now only mars_sol2ls is available, but easy to add one: put calculations in ppcompute, and add the option in the method performtimechange.

File size: 81.0 KB
Line 
1###############################################
2## PLANETOPLOT                               ##
3## --> PPCLASS                               ##
4## A generic and versatile Python module     ##
5## ... to read netCDF files and plot         ##
6###############################################
7## Author: Aymeric Spiga. 02-03/2013         ##
8###############################################
9# python built-in librairies
10import os
11import time as timelib
12import pickle
13# added librairies
14import numpy as np
15import netCDF4
16import matplotlib.pyplot as mpl
17# personal librairies
18import ppplot
19import ppcompute
20###############################################
21
22###################################
23#### HEADER                      ##
24#### ... executed when imported  ##
25###################################
26# where settings files are located...
27whereset = None
28whereset = ppcompute.findset(whereset)
29# ... we load user-defined automatic settings from set_ppclass.txt
30zefile = "set_ppclass.txt"
31glob_listx = [] ; glob_listy = [] ; glob_listz = [] ; glob_listt = []
32glob_listarea = []
33try: 
34    f = open(whereset+zefile, 'r') ; lines = f.readlines()
35    for stuff in lines[5].strip().split(';'): glob_listx.append(stuff)
36    for stuff in lines[8].strip().split(';'): glob_listy.append(stuff)
37    for stuff in lines[11].strip().split(';'): glob_listz.append(stuff)
38    for stuff in lines[14].strip().split(';'): glob_listt.append(stuff)
39    for stuff in lines[17].strip().split(';'): glob_listarea.append(stuff)
40except IOError: 
41    print "warning: "+zefile+" not in "+whereset+" ; no presets."
42
43##################################
44#### USEFUL GENERIC FUNCTIONS ####
45##################################
46
47# inspect variables and dimensions in a netCDF file
48def inspect(filename):
49    print "**** INSPECT FILE",filename
50    test = netCDF4.Dataset(filename)
51    print "**** VARIABLES: ",test.variables.keys()
52    for dim in test.dimensions.keys():
53        output = "**** DIMENSION: "+str(dim)+" "+str(len(test.dimensions[dim]))
54        try: output = output + " ----> "+str(test.variables[dim][0])+"  "+str(test.variables[dim][-1])
55        except: pass
56        print output ; output = ""
57
58# check a tab and exit if wrong. if just one string make it a list.
59# (if allownumber, convert this into a string).
60def checktab(tab,mess="",allownone=False,allownumber=False):
61    if tab is None: 
62      if not allownone:  print "pp.define: no "+mess ; exit()
63      else: pass
64    else:
65      if not isinstance(tab, list):
66        if isinstance(tab, str): 
67            tab = [tab]
68        elif (isinstance(tab, int) or isinstance(tab, float)) and allownumber: 
69            tab = [str(tab)] 
70        else: 
71            print "pp.define: "+mess+" should be either a string or a list of strings!" ; exit()
72      elif isinstance(tab, list):
73        if isinstance(tab[0],str): 
74            pass
75        elif (isinstance(tab[0], int) or isinstance(tab[0], float)) and allownumber:
76            for iii in range(len(tab)): tab[iii] = str(tab[iii])
77        else: 
78            print "pp.define: "+mess+" should be either a string or a list of strings!" ; exit()
79    return tab
80
81# determine which method is to be applied to a given dimension
82def findmethod(tab):
83    if tab is None:              output = "free"
84    elif tab[0,0] != tab[0,1]:   output = "comp"
85    else:                        output = "fixed"
86    return output
87
88# read what is given by the user (version of T. Navarro)
89def readslices(saxis):
90    if saxis == None:
91        zesaxis = None
92    else:
93        zesaxis = np.empty((len(saxis),2))
94        for i in range(len(saxis)):
95            a = separatenames(saxis[i])
96            if len(a) == 1:
97                zesaxis[i,:] = float(a[0])
98            else:
99                zesaxis[i,0] = float(a[0])
100                zesaxis[i,1] = float(a[1])         
101    return zesaxis
102
103# look for comas in the input name to separate different names (files, variables,etc ..)
104# (needed by readslices)
105def separatenames (name):
106    if name is None: names = None
107    else:
108      names = [] ; stop = 0 ; currentname = name
109      while stop == 0:
110        indexvir = currentname.find(',')
111        if indexvir == -1: stop = 1 ; name1 = currentname
112        else: name1 = currentname[0:indexvir]
113        names = np.concatenate((names,[name1]))
114        currentname = currentname[indexvir+1:len(currentname)]
115    return names
116
117#######################
118### THE MAIN OBJECT ###
119#######################
120class pp():
121
122    # print out a help string when help is invoked on the object
123    def __repr__(self):
124        whatprint = 'pp object. \"help(pp)\" for more information\n'
125        return whatprint
126
127    # default settings
128    # -- user can define settings by two methods.
129    # -- 1. yeah = pp(file="file.nc")
130    # -- 2. yeah = pp() ; yeah.file = "file.nc"
131    def __init__(self,file=None,var="notset",\
132                      filegoal=None,vargoal=None,\
133                      x=None,y=None,z=None,t=None,\
134                      stridex=1,stridey=1,\
135                      stridez=1,stridet=1,\
136                      compute="mean",\
137                      verbose=False,noproj=False,\
138                      superpose=False,\
139                      plotin=None,\
140                      forcedimplot=-1,\
141                      out="gui",\
142                      filename="myplot",\
143                      folder="./",\
144                      includedate=True,\
145                      xlabel=None,ylabel=None,\
146                      xcoeff=None,ycoeff=None,\
147                      proj=None,\
148                      vmin=None,vmax=None,\
149                      div=None,\
150                      colorb=None,\
151                      lstyle=None,\
152                      marker=None,\
153                      color=None,\
154                      label=None,\
155                      changetime=None,\
156                      title=None):
157        self.request = None
158        self.nfin = 0 ; self.nvin = 0
159        self.nplotx = None ; self.nploty = None
160        self.nplotz = None ; self.nplott = None
161        self.status = "init"
162        self.fig = None ; self.subv = None ; self.subh = None 
163        self.n = 0 ; self.howmanyplots = 0
164        self.nplot = 0
165        self.p = None
166        self.customplot = False
167        self.allfield = None
168        ## what could be defined by the user
169        self.file = file
170        self.var = var
171        self.filegoal = filegoal
172        self.vargoal = vargoal
173        self.x = x ; self.y = y   ## if None, free dimension
174        self.z = z ; self.t = t   ## if None, free dimension
175        self.stridex = stridex ; self.stridey = stridey
176        self.stridez = stridez ; self.stridet = stridet
177        self.compute = compute
178        self.verbose = verbose
179        self.noproj = noproj
180        self.plotin = plotin
181        self.superpose = superpose
182        self.forcedimplot = forcedimplot
183        self.out = out
184        self.filename = filename
185        self.folder = folder
186        self.includedate = includedate
187        self.changetime = changetime
188        ## here are user-defined plot settings
189        ## -- if not None, valid on all plots in the pp() objects
190        self.xlabel = xlabel ; self.xcoeff = xcoeff
191        self.ylabel = ylabel ; self.ycoeff = ycoeff
192        self.proj = proj
193        self.vmin = vmin ; self.vmax = vmax
194        self.div = div
195        self.colorb = colorb
196        self.lstyle = lstyle
197        self.marker = marker
198        self.color = color
199        self.label = label
200        self.title = title
201
202    # print status
203    def printstatus(self):
204        if self.filename == "THIS_IS_A_CLONE":
205            pass
206        else:
207            print "**** PPCLASS. Done step: " + self.status
208
209    #####################################################
210    # EMULATE OPERATORS + - * / ** << FOR PP() OBJECTS  #
211    #####################################################
212
213    # define the operation <<
214    # ... e.g. obj2 << obj1
215    # ... means: get init for pp object obj2 from another pp object obj1
216    # ... (this should solve the affectation trap obj2 = obj1)
217    def __lshift__(self,other):
218        if other.__class__.__name__ == "pp":
219            self.file = other.file
220            self.var = other.var
221            self.filegoal = other.filegoal
222            self.vargoal = other.vargoal
223            self.x = other.x ; self.y = other.y   ## if None, free dimension
224            self.z = other.z ; self.t = other.t   ## if None, free dimension
225            self.stridex = other.stridex ; self.stridey = other.stridey
226            self.stridez = other.stridez ; self.stridet = other.stridet
227            self.verbose = other.verbose
228            self.noproj = other.noproj
229            self.plotin = other.plotin
230            self.superpose = other.superpose
231            self.forcedimplot = other.forcedimplot
232            self.out = other.out
233            self.filename = other.filename
234            self.folder = other.folder
235            self.xlabel = other.xlabel ; self.xcoeff = other.xcoeff
236            self.ylabel = other.ylabel ; self.ycoeff = other.ycoeff
237            self.proj = other.proj
238            self.vmin = other.vmin ; self.vmax = other.vmax
239            self.div = other.div
240            self.colorb = other.colorb
241            self.lstyle = other.lstyle
242            self.marker = other.marker
243            self.color = other.color
244            self.label = other.label
245            self.title = other.title
246            self.includedate = other.includedate
247            self.changetime = other.changetime
248        else:
249            print "!! ERROR !! argument must be a pp object." ; exit()
250
251    # check the compatibility of two objects for operations
252    # --> if other is a pp class, test sizes and return isnum = False
253    # --> if other is an int or a float, return isnum = True
254    # --> otherwise, just print an error and exit
255    def checktwo(self,other):
256        if other.__class__.__name__ == "pp":
257          isnum = False
258          if self.status in ["init","defined"] or other.status in ["init","define"]: 
259             print "!! ERROR !! Please use .retrieve to get fields for plots with one of your pp operands." ; exit()
260          if self.nfin   != other.nfin   or \
261             self.nvin   != other.nvin   or \
262             self.nplott != other.nplott or \
263             self.nplotz != other.nploty or \
264             self.nploty != other.nploty or \
265             self.nplotx != other.nplotx :
266               print "!! ERROR !! The two operands do not have the same number of files, variables, t z y x requests."
267               exit()
268        elif isinstance(other,int) or isinstance(other,float):
269          isnum = True
270        else:
271          print "!! ERROR !! The operand is neither a pp class nor an integer or a float." ; exit()
272        return isnum
273
274    # define a selective copy of a pp() object for operations
275    # ... copy.copy() is not conservative (still acts like a pointer)
276    # ... copy.deepcopy() does not work with netCDF objects
277    # so what is done here is a copy of everything except
278    # (to avoid sharing with self and therefore modifying self through operations)
279    # - request attribute of pp() object
280    # - field attribute of the onerequest() objects
281    def selective_copy(self):
282        if self.status in ["init","defined"]:
283            print "!! ERROR !! Please use .retrieve to get fields for the object you want to copy from." ; exit()
284        the_clone = pp()
285        for k, v in vars(self).items():
286           if k != "request":
287               setattr(the_clone,k,v)
288        the_clone.verbose = False
289        the_clone.filename = "THIS_IS_A_CLONE" # trick to avoid additional outputs
290        the_clone.define()
291        for i in range(self.nfin):
292         for j in range(self.nvin):
293          for t in range(self.nplott):
294           for z in range(self.nplotz):
295            for y in range(self.nploty):
296             for x in range(self.nplotx):
297              obj_ref = self.request[i][j][t][z][y][x]
298              obj = the_clone.request[i][j][t][z][y][x]
299              for k, v in vars(obj_ref).items():
300               if k != "field":
301                setattr(obj,k,v)
302        the_clone.status = "retrieved"
303        the_clone.filename = self.filename
304        return the_clone
305
306    # define the operation + on two objects. or with an int/float.
307    # ... with selective_copy the self object is not modified.
308    def __add__(self,other):
309        isnum = self.checktwo(other)
310        the_clone = self.selective_copy()
311        for i in range(self.nfin):
312         for j in range(self.nvin):
313          for t in range(self.nplott):
314           for z in range(self.nplotz):
315            for y in range(self.nploty):
316             for x in range(self.nplotx):
317              obj = the_clone.request[i][j][t][z][y][x]
318              obj_ref = self.request[i][j][t][z][y][x]
319              if not isnum:   
320                  ope = other.request[i][j][t][z][y][x].field
321                  if obj_ref.field.shape != ope.shape:
322                    print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
323                    exit()
324              else:           
325                  ope = other
326              goal = self.vargoal[j] + self.filegoal[i]
327              if ("vector" in goal) or ("contour" in goal):
328                  if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
329                  obj.field = obj_ref.field
330              else:
331                  obj.field = obj_ref.field + ope
332        return the_clone
333
334    # define the operation - on two objects. or with an int/float.
335    # ... with selective_copy the self object is not modified.
336    def __sub__(self,other):
337        isnum = self.checktwo(other)
338        the_clone = self.selective_copy()
339        for i in range(self.nfin):
340         for j in range(self.nvin):
341          for t in range(self.nplott):
342           for z in range(self.nplotz):
343            for y in range(self.nploty):
344             for x in range(self.nplotx):
345              obj = the_clone.request[i][j][t][z][y][x]
346              obj_ref = self.request[i][j][t][z][y][x]
347              if not isnum:
348                  ope = other.request[i][j][t][z][y][x].field
349                  if obj_ref.field.shape != ope.shape:
350                    print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
351                    exit()
352              else:
353                  ope = other
354              goal = self.vargoal[j] + self.filegoal[i]
355              if ("vector" in goal) or ("contour" in goal):
356                  if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
357                  obj.field = obj_ref.field
358              else:
359                  obj.field = obj_ref.field - ope
360        return the_clone
361
362    # define the operation * on two objects. or with an int/float.
363    # ... with selective_copy the self object is not modified.
364    def __mul__(self,other):
365        isnum = self.checktwo(other)
366        the_clone = self.selective_copy()
367        for i in range(self.nfin):
368         for j in range(self.nvin):
369          for t in range(self.nplott):
370           for z in range(self.nplotz):
371            for y in range(self.nploty):
372             for x in range(self.nplotx):
373              obj = the_clone.request[i][j][t][z][y][x]
374              obj_ref = self.request[i][j][t][z][y][x]
375              if not isnum:
376                  ope = other.request[i][j][t][z][y][x].field
377                  if obj_ref.field.shape != ope.shape:
378                    print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
379                    exit()
380              else:
381                  ope = other
382              goal = self.vargoal[j] + self.filegoal[i]
383              if ("vector" in goal) or ("contour" in goal):
384                  if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
385                  obj.field = obj_ref.field
386              else:
387                  obj.field = obj_ref.field * ope
388        return the_clone
389
390    # define the operation / on two objects. or with an int/float.
391    # ... with selective_copy the self object is not modified.
392    def __div__(self,other):
393        isnum = self.checktwo(other)
394        the_clone = self.selective_copy()
395        for i in range(self.nfin):
396         for j in range(self.nvin):
397          for t in range(self.nplott):
398           for z in range(self.nplotz):
399            for y in range(self.nploty):
400             for x in range(self.nplotx):
401              obj = the_clone.request[i][j][t][z][y][x]
402              obj_ref = self.request[i][j][t][z][y][x]
403              if not isnum:
404                  ope = other.request[i][j][t][z][y][x].field
405                  if obj_ref.field.shape != ope.shape:
406                    print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
407                    exit()
408              else:
409                  ope = other
410              goal = self.vargoal[j] + self.filegoal[i]
411              if ("vector" in goal) or ("contour" in goal):
412                  if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
413                  obj.field = obj_ref.field
414              else:
415                  obj.field = obj_ref.field / ope
416        return the_clone
417
418    # define the reverse operation float/int + object
419    def __radd__(self,other):
420        isnum = self.checktwo(other)
421        if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
422        return self.__add__(other)
423
424    # define the reverse operation float/int - object
425    def __rsub__(self,other):
426        isnum = self.checktwo(other)
427        if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
428        return self.__sub__(other)
429
430    # define the reverse operation float/int * object
431    def __rmul__(self,other):
432        isnum = self.checktwo(other)
433        if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
434        return self.__mul__(other)
435
436    # define the reverse operation float/int / object
437    def __rdiv__(self,other):
438        isnum = self.checktwo(other)
439        if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
440        return self.__div__(other)
441
442    # define the operation ** on one object.
443    # ... with selective_copy the self object is not modified.
444    def __pow__(self,num):
445        the_clone = self.selective_copy()
446        if isinstance(num,int) or isinstance(num,float):
447            for i in range(self.nfin):
448             for j in range(self.nvin):
449              for t in range(self.nplott):
450               for z in range(self.nplotz):
451                for y in range(self.nploty):
452                 for x in range(self.nplotx):
453                  obj  = the_clone.request[i][j][t][z][y][x]
454                  obj_ref = self.request[i][j][t][z][y][x]
455                  goal = self.vargoal[j] + self.filegoal[i]
456                  if ("vector" in goal) or ("contour" in goal):
457                      if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
458                      obj.field = obj_ref.field
459                  else:
460                      obj.field = obj_ref.field ** num
461        else:
462            print "!! ERROR !! To define a power, either an int or a float is needed." ; exit()
463        return the_clone
464
465    ### TBD: reverse power? for exponentials?
466
467    ##############################################################################################
468    # define method
469    # ---------
470    # ... (file and var are either one string or a vector of strings)
471    # ... the goal of define is to define a 2D array of onerequest() objects (see class below)
472    #     given the number of file, var, x, y, z, t asked by the user
473    # ... objectives for file or var are given through filegoal and vargoal
474    #     --> possible values: main contour vector
475    # ---------
476    # ... then onerequest() objects are being defined more precisely
477    #     by getting index_x index_y index_z index_t
478    #     and setting method_x method_y method_z method_t to either
479    #      - "free" for free dimensions (plot dimensions)
480    #      - "comp" for averages, max, min
481    #      - "fixed" for fixed dimensions (possibly several i.e. multislice)
482    ##############################################################################################
483    def define(self):
484        self.printstatus()
485        # initial check and get dimensions
486        self.file = checktab(self.file,mess="file")
487        self.nfin = len(self.file)
488        if self.verbose:
489            for i in range(self.nfin): inspect(self.file[i])
490        self.var = checktab(self.var,mess="var")
491        self.nvin = len(self.var)
492        # check goal tabs for files and variables
493        # ... default is to plot everything
494        if self.filegoal is None: self.filegoal = ["main"]*self.nfin
495        if self.vargoal is None:  self.vargoal  = ["main"]*self.nvin
496        self.filegoal = checktab(self.filegoal, mess="filegoal")
497        self.vargoal  = checktab(self.vargoal,  mess="vargoal")
498        if len(self.filegoal) != self.nfin:  print "!! ERROR !! filegoal must be the same size as file." ; exit()
499        if len(self.vargoal)  != self.nvin:  print "!! ERROR !! vargoal must be the same size as var." ; exit()
500        # variables: initial check
501        self.x = checktab(self.x,mess="x",allownone=True,allownumber=True)
502        self.y = checktab(self.y,mess="y",allownone=True,allownumber=True)
503        self.z = checktab(self.z,mess="z",allownone=True,allownumber=True)
504        self.t = checktab(self.t,mess="t",allownone=True,allownumber=True)
505        # for the moment not var- nor file- dependent.
506        # but this could be the case.
507        sx = readslices(self.x) ; sy = readslices(self.y)
508        sz = readslices(self.z) ; st = readslices(self.t)
509        # get methods
510        mx = findmethod(sx) ; my = findmethod(sy)
511        mz = findmethod(sz) ; mt = findmethod(st)
512        # get number of plots to be done
513        if mx in ["fixed","comp"]: self.nplotx = sx.size/2
514        else:                      self.nplotx = 1
515        if my in ["fixed","comp"]: self.nploty = sy.size/2
516        else:                      self.nploty = 1
517        if mz in ["fixed","comp"]: self.nplotz = sz.size/2
518        else:                      self.nplotz = 1
519        if mt in ["fixed","comp"]: self.nplott = st.size/2
520        else:                      self.nplott = 1
521        if self.verbose:  print "**** OK. Plots over x,y,z,t -->",self.nplotx,self.nploty,self.nplotz,self.nplott
522        # create the list of onerequest() objects
523        self.request = [[[[[[ \
524                       onerequest() \
525                       for x in range(self.nplotx)] for y in range(self.nploty)] \
526                       for z in range(self.nplotz)] for t in range(self.nplott)] \
527                       for j in range(self.nvin)]   for i in range(self.nfin)] 
528        # loop on onerequest() objects
529        for i in range(self.nfin):
530         for j in range(self.nvin):
531          for t in range(self.nplott):
532           for z in range(self.nplotz):
533            for y in range(self.nploty):
534             for x in range(self.nplotx):
535              obj = self.request[i][j][t][z][y][x]
536              # fill in names for files and variables
537              obj.verbose = self.verbose
538              obj.file = self.file[i]
539              obj.var = self.var[j]
540              # indicate the computation method
541              obj.compute = self.compute
542              # open the files (the same file might be opened several times but this is cheap)
543              obj.openfile()
544              ### get x,y,z,t dimensions from file
545              obj.getdim()
546              ### possible time axis change
547              obj.changetime = self.changetime
548              obj.performtimechange()
549              ### get methods
550              obj.method_x = mx ; obj.method_y = my
551              obj.method_z = mz ; obj.method_t = mt           
552              ### get index
553              obj.getindextime(dalist=st,ind=t,stride=self.stridet)
554              obj.getindexvert(dalist=sz,ind=z,stride=self.stridez)
555              obj.getindexhori(dalistx=sx,dalisty=sy,indx=x,indy=y,stridex=self.stridex,stridey=self.stridey)
556        # change status
557        self.status = "defined"
558        return self
559
560    ##############################################################################################
561    # retrieve method
562    # --> for each element onerequest() in the array, get field .var from .f file
563    # --> see below the onerequest() class:
564    #        - only get what is needed for computing and plotting
565    #        - averages etc... are computed here
566    # --> RESULT: each onerequest() object has now its attribute .field filled
567    # --> if one wants to perform operations on fields, this should be done after retrieve()
568    ##############################################################################################
569    def retrieve(self):
570        self.printstatus()
571        # check if things were done OK before
572        if self.status != "defined": print "!! ERROR !! Please use .define() to define your pp object." ; exit()
573        # create the list of allfield() objects
574        # --> so that the user can easily access values
575        self.allfield = [[[[[[ \
576                        [] \
577                        for x in range(self.nplotx)] for y in range(self.nploty)] \
578                        for z in range(self.nplotz)] for t in range(self.nplott)] \
579                        for j in range(self.nvin)]   for i in range(self.nfin)]
580        ## first get fields
581        ## ... only what is needed is extracted from the files
582        ## ... and computations are performed
583        for i in range(self.nfin):
584         for j in range(self.nvin):
585          for t in range(self.nplott):
586           for z in range(self.nplotz):
587            for y in range(self.nploty):
588             for x in range(self.nplotx):
589              obj = self.request[i][j][t][z][y][x]
590              obj.getfield()
591              obj.computations()
592              self.allfield[i][j][t][z][y][x] = obj.field
593        # change status
594        self.status = "retrieved"
595        return self
596
597    ##########################################################
598    # get: a shortcut method for the define + retrieve chain #
599    ##########################################################
600    def get(self):
601        self.define()
602        self.retrieve()
603        return self 
604
605    ########################################
606    # smooth: smooth the field in 1D or 2D #
607    ########################################
608    ## TBD: smooth not OK with masked array in the end of retrieve()
609    def smooth(self,window):
610        if self.verbose: 
611            print "!! WARNING !! Performing a smoothing with a window size",window
612            print "!! WARNING !! To come back to unsmoothed file, use .get() again"
613        for i in range(self.nfin):
614         for j in range(self.nvin):
615          for t in range(self.nplott):
616           for z in range(self.nplotz):
617            for y in range(self.nploty):
618             for x in range(self.nplotx):
619              obj = self.request[i][j][t][z][y][x]
620              if obj.field.ndim == 1:
621                  print "!! ERROR !! 1D smoothing not supported yet because reduces array sizes."
622                  exit()
623                  # TBD TBD TBD
624                  #obj.field = ppcompute.smooth1d(obj.field,window=window)
625              elif obj.field.ndim == 2:
626                  obj.field = ppcompute.smooth2d(obj.field,window=window)
627
628    ############################################################################################## 
629    # defineplot method
630    # --> defineplot first defines arrays of plot objects and set each of them
631    #     ... simple looping except cases where goal is not main (e.g. contour or vector)
632    # --> principle: each onerequest() object gives birth to a subplot
633    # --> defineplot vs. makeplot: defining plot and actually plotting it are clearly separated
634    # --> THE KEY OUPUT OF defineplot IS AN ARRAY self.p OF PLOT OBJECTS
635    # optional arguments
636    # --> extraplot: to indicate a number of plots to be added afterwards (use self.plotin)
637    # --> loadfile: to use self.p from a previously saved file
638    ##############################################################################################
639    def defineplot(self,extraplot=0,loadfile=None):
640        # -----------------------------------------------------
641        # LOAD MODE: load a self.p object. count plots from it.
642        # -----------------------------------------------------
643        if loadfile is not None:
644            try: filehandler = open(loadfile, 'r') ; self.p = pickle.load(filehandler) 
645            except IOError: print "!! ERROR !! Cannot find object file to load." ; exit()
646            self.status = "definedplot" ; self.plotin = None
647            self.nplot = len(self.p) ; self.howmanyplots = self.nplot
648            return
649        # -----------------------------------------------------
650        # REGULAR MODE
651        # -----------------------------------------------------
652        self.printstatus()
653        # check if things were done OK before
654        if self.status in ["init","defined"]: 
655            print "!! ERROR !! Please use .retrieve() to get fields for plots with your pp object." ; exit()
656        # check self.plotin (an existing fig on which to add plots afterwards)
657        if self.plotin.__class__.__name__ == "pp":
658            if self.plotin.fig is None:
659                self.plotin = None # this is an additional security in case
660                                   #   a pp object is given without figure opened yet.
661        elif self.plotin is not None:
662            print "!! ERROR !! plotin argument must be a pp object." ; exit()
663        # initialize the array of subplot objects
664        # either something new or attributes coming from plotin object
665        if self.plotin is None:  self.p = [ ]
666        else:                    self.p = self.plotin.p
667        # create an array of subplot objects
668        # ... in theory the order of looping can be changed without any harm
669        # ... the only important thing is to keep i,j,t,z,y,x resp. for file,var,t,z,y,x
670        count = 0
671        for i in range(self.nfin):
672         if self.filegoal[i] == "main": 
673          for j in range(self.nvin):
674           if self.vargoal[j] == "main":
675            for t in range(self.nplott):
676             for z in range(self.nplotz):
677              for y in range(self.nploty):
678               for x in range(self.nplotx):
679                # look at dimension and append the right plot object
680                obj = self.request[i][j][t][z][y][x]
681                dp = obj.dimplot
682                if dp == 1 or self.forcedimplot == 1:    plobj = ppplot.plot1d()
683                elif dp == 2 or self.forcedimplot == 2:  plobj = ppplot.plot2d()
684                elif dp == 0: print "**** OK. VALUES VALUES VALUES",obj.field
685                else:         print "!! ERROR !! 3D or 4D plots not supported" ; exit()
686                # load abscissa and ordinate in obj
687                obj.definecoord()
688                # we start to define things here before appending
689                # (convenient: could be overridden by the user before makeplot)
690                # ... the if loop is necessary so that we can loop above on the dp=0 case
691                if dp in [1,2]:
692                    # and define what to do in plobj
693                    plobj.invert = obj.invert_axes
694                    plobj.swap = obj.swap_axes
695                    # axis labels
696                    plobj.xlabel = obj.absclab ; plobj.ylabel = obj.ordilab
697                    # superpose or not (this is mostly for saving purpose)
698                    plobj.superpose = self.superpose
699                    # get title, colormaps, labels, etc.. from var
700                    plobj.var = obj.var
701                    plobj.define_from_var()
702                    # generic 1D/2D: load field and coord in plot object
703                    plobj.field = obj.field    # field to be plotted
704                    plobj.absc = obj.absc      # abscissa (or longitude)
705                    plobj.ordi = obj.ordi      # ordinate (or latitude)
706                                               # -- useless in 1D but not used anyway
707                    # specific 1D plot stuff
708                    if dp == 1:
709                        # -- a default label
710                        plobj.label = ""
711                        if self.nfin > 1: plobj.label = plobj.label + " file #"+str(i+1)
712                        if self.nvin > 1: plobj.label = plobj.label + " var #"+str(j+1)
713                        if self.nplott > 1: plobj.label = plobj.label + " t #"+str(t+1)
714                        if self.nplotz > 1: plobj.label = plobj.label + " z #"+str(z+1)
715                        if self.nploty > 1: plobj.label = plobj.label + " y #"+str(y+1)
716                        if self.nplotx > 1: plobj.label = plobj.label + " x #"+str(x+1)
717                    # specific 2d plot stuff
718                    if dp == 2:
719                        # -- light grey background for missing values
720                        if type(plobj.field).__name__ in 'MaskedArray': plobj.axisbg = '0.75'
721                        # -- define if it is a map or a plot
722                        plobj.mapmode = ( obj.method_x+obj.method_y == "freefree" \
723                                       and "grid points" not in obj.name_x \
724                                       and not self.noproj )
725                    # possible user-defined plot settings shared by all plots
726                    if self.div is not None: plobj.div = self.div
727                    if self.xlabel is not None: plobj.xlabel = self.xlabel
728                    if self.xcoeff is not None: plobj.xcoeff = self.xcoeff
729                    if self.ylabel is not None: plobj.ylabel = self.ylabel
730                    if self.ycoeff is not None: plobj.ycoeff = self.ycoeff
731                    if self.title is not None: plobj.title = self.title
732                    # -- 1D specific
733                    if dp == 1:
734                        if self.lstyle is not None: plobj.lstyle = self.lstyle
735                        if self.marker is not None: plobj.marker = self.marker
736                        if self.color is not None: plobj.color = self.color
737                        if self.label is not None: plobj.label = self.label
738                    # -- 2D specific
739                    elif dp == 2:
740                        if self.proj is not None and not self.noproj: plobj.proj = self.proj
741                        if self.vmin is not None: plobj.vmin = self.vmin
742                        if self.vmax is not None: plobj.vmax = self.vmax
743                        if self.colorb is not None: plobj.colorb = self.colorb
744                    # finally append plot object
745                    self.p.append(plobj)
746                    count = count + 1
747        # self.nplot is number of plot to be defined in this call to defineplot()
748        # (because of self.plotin this might less than length of self.p)
749        self.nplot = count
750        # --- superimposed contours and vectors ---
751        # we have to start another loop because we need forward information
752        # TBD: there is probably a more flexible way to do that
753        count = 0
754        for i in range(self.nfin):
755         for j in range(self.nvin):
756          for t in range(self.nplott):
757           for z in range(self.nplotz):
758            for y in range(self.nploty):
759             for x in range(self.nplotx):
760              goal = self.vargoal[j] + self.filegoal[i]
761              obj = self.request[i][j][t][z][y][x]
762              if "mainmain" in goal and obj.dimplot == 2:
763                  # the plot object we consider in the loop
764                  pl = self.p[count]
765                  # -- see if there is a contour requested...
766                  # (we use try because we might be at the end of the list)
767                  found = 0
768                  try:    condvar = self.vargoal[j+1]
769                  except: condvar = "itisok"
770                  try:    condfile = self.filegoal[i+1]
771                  except: condfile = "itisok"
772                  # ... get contour
773                  ##########################################
774                  # NB: contour is expected to be right after main otherwise it is not displayed
775                  ##########################################
776                  if condvar == "contour":
777                      plobj.addcontour = self.request[i][j+1][t][z][y][x].field ; found += 1
778                  if condfile == "contour":
779                      plobj.addcontour = self.request[i+1][j][t][z][y][x].field ; found += 1
780                  # see if there is a vector requested...
781                  # (we use try because we might be at the end of the list)
782                  try:    condvar = self.vargoal[j+found+1]+self.vargoal[j+found+2]
783                  except: condvar = "itisok"
784                  try:    condfile = self.filegoal[i+found+1]+self.filegoal[i+found+2]
785                  except: condfile = "itisok"
786                  # ... get vector and go directly to the next iteration
787                  # (in some cases we would do this twice but this is cheap)
788                  if "vector" in condvar:
789                      plobj.addvecx = self.request[i][j+found+1][t][z][y][x].field
790                      plobj.addvecy = self.request[i][j+found+2][t][z][y][x].field
791                  if "vector" in condfile:
792                      plobj.addvecx = self.request[i+found+1][j][t][z][y][x].field
793                      plobj.addvecy = self.request[i+found+2][j][t][z][y][x].field
794                  count = count + 1
795        # COUNT PLOTS. if 0 just exit.
796        # self.howmanyplots is self.nplot + possible extraplots
797        self.howmanyplots = self.nplot + extraplot
798        if self.howmanyplots > 0:
799            if self.verbose: print "**** OK. expect %i plots" % (self.howmanyplots)
800        else:
801            exit() # because this means that we only had 0D values !
802        # final status
803        self.status = "definedplot"
804        return self
805
806    ##############################################################################################
807    # makeplot method
808    # --> after defineplot and before makeplot, user-defined plot settings can be easily given
809    #     simply by modifying the attributes of each elements of self.p
810    # --> to change only one plot setting, no need to call defineplot again before makeplot
811    # --> in the end, the array self.p of plot objects is saved for easy and convenient replotting
812    # --> see practical examples in the folder 'examples'
813    ##############################################################################################
814    def makeplot(self):
815        self.printstatus()
816        # a few initial operations
817        # ------------------------
818        if "definedplot" not in self.status: 
819            print "!! ERROR !! Please use .defineplot() before .makeplot() can be used with your pp object." ; exit()
820        # open a figure and define subplots         
821        # ---------------------------------
822        if self.plotin is None: 
823            # start from scratch
824            self.fig = mpl.figure(figsize=(16,8))
825            self.subv,self.subh = ppplot.definesubplot(self.howmanyplots,self.fig) 
826            self.n = 0
827            ## adapted space for labels etc
828            ## ... except for ortho because there is no label anyway
829            self.customplot = self.p[0].field.ndim == 2 \
830                        and self.p[0].mapmode == True \
831                        and self.p[0].proj not in ["ortho"]
832            if self.customplot:
833                margin = 0.07
834                self.fig.subplots_adjust(left=margin,right=1-margin,bottom=margin,top=1-margin)
835        else:
836            # start from an existing figure.
837            # extraplot must have been set in the call to the previous figure.
838            self.fig = self.plotin.fig
839            self.subv,self.subh = self.plotin.subv,self.plotin.subh
840            self.n = self.plotin.n
841            self.howmanyplots = self.plotin.howmanyplots
842            self.customplot = self.plotin.customplot
843        # LOOP on all subplots
844        # NB: cannot use 'for pl in self.p' if self.plotin not None
845        # --------------------
846        for count in range(self.nplot):
847            # the plot object we consider in the loop
848            pl = self.p[self.n]
849            # before making the plot, create a subplot. the first one is numbered 1 not 0.
850            # ... if pl.superpose, we use only one and only figure
851            # ... (and we have to be careful with not doing things several times)
852            if pl.superpose:
853                if self.n == 0: 
854                    self.fig.add_subplot(1,1,1,axisbg=pl.axisbg) # define one subplot (still needed for user-defined font sizes)
855                    sav = pl.xlabel,pl.ylabel,pl.xcoeff,pl.ycoeff,pl.title,pl.swaplab # save titles and labels
856                else: 
857                    pl.invert = False ; pl.lstyle = None # don't invert again axis
858                    # set saved titles and labels
859                    if self.plotin is None:
860                        pl.xlabel,pl.ylabel,pl.xcoeff,pl.ycoeff,pl.title,pl.swaplab = sav
861                    else:
862                        prev_plot = self.plotin.p[self.n-1]
863                        pl.xlabel = prev_plot.xlabel
864                        pl.ylabel = prev_plot.ylabel
865                        pl.xcoeff = prev_plot.xcoeff
866                        pl.ycoeff = prev_plot.ycoeff
867                        pl.title = prev_plot.title
868                        pl.swaplab = prev_plot.swaplab
869            else:
870                self.fig.add_subplot(self.subv,self.subh,self.n+1,axisbg=pl.axisbg)
871            if self.verbose: print "**** Done subplot %i / %i " %( self.n+1,self.howmanyplots ) 
872            # finally make the plot
873            pl.make()
874            # increment plot count (and propagate this in plotin)
875            self.n = self.n+1
876            if self.plotin is not None: self.plotin.n = self.n
877        # once completed show the plot (cannot show intermediate plotin)
878        # ... added a fix (customplot=True) for the label problem in basemap
879        print "**** PPCLASS. Done step: makeplot"
880        if (self.n == self.howmanyplots):
881            ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=self.customplot,includedate=self.includedate)
882            mpl.close()
883        # SAVE A PICKLE FILE WITH THE self.p ARRAY OF OBJECTS
884        if self.verbose: print "**** Saving session in "+self.filename + ".ppobj"
885        savfile = self.folder + "/" + self.filename + ".ppobj"
886        try: 
887            filehandler = open(savfile, 'w')
888            pickle.dump(self.p, filehandler)
889        except IOError: 
890            print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
891        return self
892
893    ###########################################################
894    # plot: a shortcut method for the defineplot + plot chain #
895    ###########################################################
896    def plot(self,extraplot=0):
897        self.defineplot(extraplot=extraplot)
898        self.makeplot()
899        return self
900
901    #######################################################
902    # getplot: a shortcut method for the get + plot chain #
903    #######################################################
904    def getplot(self,extraplot=0):
905        self.get()
906        self.plot(extraplot=extraplot)
907        return self
908
909    ###################################################################
910    # getdefineplot: a shortcut method for the get + defineplot chain #
911    ###################################################################
912    def getdefineplot(self,extraplot=0):
913        self.get()
914        self.defineplot(extraplot=extraplot)
915        return self
916
917    ##############################################################
918    # f: operation on two pp objects being on status 'definedplot'
919    # this allows for one field being function of another one
920    # e.g. u.f(v) means u will be displayed as a function of v
921    # ... no need to do defineplot after u.f(v), makeplot directly
922    ##############################################################
923    def f(self,other):
924        # preamble: for this operation to work, defineplot() must have been done
925        if self.status != "definedplot":
926            if self.verbose: print "!! WARNING !! performing defineplot on operand"
927            self.defineplot()
928        if other.status != "definedplot":
929            if self.verbose: print "!! WARNING !! performing defineplot on operand"
930            other.defineplot()
931        # check total number of plots
932        if self.howmanyplots != other.howmanyplots:
933               print "!! ERROR !! The two operands do not have the same number of subplots."
934               exit()
935        # and now operation.
936        count = 0
937        while count < self.howmanyplots:
938           sobj = self.p[count] ; oobj = other.p[count]
939           if sobj.field.ndim !=1 or oobj.field.ndim !=1:
940               if self.verbose: print "!! WARNING !! Flattening arrays because more than one-dimensional."
941               sobj.field = np.ravel(sobj.field)
942               oobj.field = np.ravel(oobj.field)
943           sobj.absc = oobj.field
944           sobj.xlabel = oobj.ylabel
945           if sobj.absc.size > sobj.field.size:
946               if self.verbose:
947                   print "!! WARNING !! Trying to define y=f(x) with x and y not at the same size.",sobj.absc.size,sobj.field.size
948                   print "!! WARNING !! Modifying x to fit y size but please check." 
949               sobj.absc = sobj.absc[0:sobj.field.size]
950           count = count + 1
951        return self
952
953    ###########################################################
954    # copyopt: get options from e.g. a parser
955    # ... allow for simple scripting and user-defined settings
956    # ... must be called between defineplot and makeplot
957    # REQUIRED: attributes of opt must be the same as in the pp object
958    ###########################################################
959    def getopt(self,opt):
960        # -- if only one, or less than the number of plots --> we take the first one
961        # -- if as many as number of plots --> OK, each plot has its own setting
962        # (except a few cases such as trans)
963        for iii in range(self.howmanyplots):
964            ###
965            try: self.p[iii].trans = opt.trans
966            except: pass
967            ###
968            try: self.p[iii].div = opt.div
969            except: pass
970            ###
971            try: self.p[iii].logy = opt.logy
972            except: pass
973            ###
974            try: self.p[iii].colorb = opt.colorb[iii]
975            except: 
976                try: self.p[iii].colorb = opt.colorb[0]
977                except: pass
978            ###
979            try: self.p[iii].title = opt.title[iii]
980            except: 
981                try: self.p[iii].title = opt.title[0]
982                except: pass
983            ###
984            try: self.p[iii].xlabel = opt.xlabel[iii]
985            except: 
986                try: self.p[iii].xlabel = opt.xlabel[0]
987                except: pass
988            ###
989            try: self.p[iii].ylabel = opt.ylabel[iii]
990            except: 
991                try: self.p[iii].ylabel = opt.ylabel[0]
992                except: pass
993            ###
994            try: self.p[iii].lstyle = opt.lstyle[iii]
995            except: 
996                try: self.p[iii].lstyle = opt.lstyle[0]
997                except: pass
998            ###
999            try: self.p[iii].color = opt.color[iii]
1000            except: 
1001                try: self.p[iii].color = opt.color[0]
1002                except: pass
1003            ###
1004            try: self.p[iii].marker = opt.marker[iii]
1005            except: 
1006                try: self.p[iii].marker = opt.marker[0]
1007                except: pass
1008            ###
1009            try: self.p[iii].label = opt.label[iii]
1010            except:
1011                try: self.p[iii].label = opt.label[0]
1012                except: pass
1013            ###
1014            try: self.p[iii].proj = opt.proj[iii]
1015            except: 
1016                try: self.p[iii].proj = opt.proj[0]
1017                except: pass
1018            ###
1019            try: self.p[iii].back = opt.back[iii]
1020            except: 
1021                try: self.p[iii].back = opt.back[0]
1022                except: pass
1023            ###
1024            try: self.p[iii].area = opt.area[iii]
1025            except: 
1026                try: self.p[iii].area = opt.area[0]
1027                except: pass
1028            ###
1029            try: self.p[iii].blon = opt.blon[iii]
1030            except: 
1031                try: self.p[iii].blon = opt.blon[0]
1032                except: pass
1033            ###
1034            try: self.p[iii].blat = opt.blat[iii]
1035            except: 
1036                try: self.p[iii].blat = opt.blat[0]
1037                except: pass
1038            ###
1039            try: self.p[iii].vmin = opt.vmin[iii]
1040            except: 
1041                try: self.p[iii].vmin = opt.vmin[0]
1042                except: pass
1043            ###
1044            try: self.p[iii].vmax = opt.vmax[iii]
1045            except: 
1046                try: self.p[iii].vmax = opt.vmax[0]
1047                except: pass
1048
1049##########################################################
1050### THE ONEREQUEST SUBOBJECT TO PP (ON WHICH IT LOOPS) ###
1051##########################################################
1052class onerequest():
1053
1054    # default settings. mostly initialized to diagnose problem, except dimplot, nplot, verbose, swap_axes, invert_axes
1055    # -------------------------------
1056    def __init__(self):
1057        self.file  = '!! file: I am not set, damned !!'
1058        self.f     = None
1059        self.dim   = None
1060        self.var   = '!! var: I am not set, damned !!'
1061        self.index_x = [] ; self.index_y = [] ; self.index_z = [] ; self.index_t = []
1062        self.index_x2d = [] ; self.index_y2d = []
1063        self.method_x = '!! method_x: I am not set, damned !!'
1064        self.method_y = '!! method_y: I am not set, damned !!'
1065        self.method_z = '!! method_z: I am not set, damned !!'
1066        self.method_t = '!! method_t: I am not set, damned !!'
1067        self.field = None
1068        self.name_x = None ; self.name_y = None ; self.name_z = None ; self.name_t = None
1069        self.dim_x = None ; self.dim_y = None ; self.dim_z = None ; self.dim_t = None
1070        self.field_x = None ; self.field_y = None ; self.field_z = None ; self.field_t = None
1071        self.dimplot = 0
1072        self.nplot = 1
1073        self.absc = None ; self.ordi = None ; self.absclab = None ; self.ordilab = None
1074        self.verbose = True
1075        self.swap_axes = False ; self.invert_axes = False
1076        self.compute = None
1077        self.changetime = None
1078
1079    # open a file. for now it is netcdf. TBD for other formats.
1080    # check that self.var is inside.
1081    # -------------------------------
1082    def openfile(self):
1083        if not os.path.exists(self.file): print '!! ERROR !! I could not find the following file: '+self.file ; exit()
1084        if not os.path.isfile(self.file): print '!! ERROR !! This does not appear to be a file: '+self.file ; exit()
1085        self.f = netCDF4.Dataset(self.file)
1086        if self.verbose: print "**** OK. Opened file "+self.file
1087        if self.var not in self.f.variables.keys(): 
1088            print '!! ERROR !! File '+self.file+' does not contain variable: '+self.var
1089            print '..... try instead with ',self.f.variables.keys() ; exit()
1090
1091    # copy attributes from another existing object
1092    # --------------------------------------------
1093    def copy(self,source):
1094        for k, v in vars(source).items():
1095            setattr(self,k,v)
1096
1097    # get x,y,z,t dimensions from NETCDF file
1098    # TBD: user could request for a specific altitude dimension
1099    # TBD: staggered variables could request specific dimensions
1100    # -------------------------------
1101    def getdim(self):
1102          # GET SIZES OF EACH DIMENSION
1103          if self.verbose: print "**** OK. Found variable "+self.var
1104          shape = self.f.variables[self.var].shape
1105          self.dim = len(shape)
1106          if self.dim == 1:
1107              if self.verbose: print "**** OK. 1D field. I assume this varies with time."
1108              self.dim_x = 1 ; self.dim_y = 1 ; self.dim_z = 1 ; self.dim_t = shape[0]
1109          elif self.dim == 2:
1110              if self.verbose: print "**** OK. 2D field. I assume this is not-time-varying lat-lon map."
1111              self.dim_x = shape[1] ; self.dim_y = shape[0] ; self.dim_z = 1 ; self.dim_t = 1
1112          elif self.dim == 3:
1113              if self.verbose: print "**** OK. 3D field. I assume this is time-varying lat-lon map."
1114              self.dim_x = shape[2] ; self.dim_y = shape[1] ; self.dim_z = 1 ; self.dim_t = shape[0]
1115          elif self.dim == 4:
1116              if self.verbose: print "**** OK. 4D field."
1117              self.dim_x = shape[3] ; self.dim_y = shape[2] ; self.dim_z = shape[1] ; self.dim_t = shape[0]
1118          # LONGITUDE. Try preset fields. If not present set grid points axis.
1119          self.name_x = "nothing"
1120          for c in glob_listx:
1121            if c in self.f.variables.keys():
1122             self.name_x = c
1123          if self.name_x == "nothing":
1124            self.field_x = np.array(range(self.dim_x))
1125            self.name_x = "x grid points"
1126          else:
1127            self.field_x = self.f.variables[self.name_x]
1128          # LATITUDE. Try preset fields. If not present set grid points axis.
1129          self.name_y = "nothing"
1130          for c in glob_listy:
1131            if c in self.f.variables.keys():
1132             self.name_y = c
1133          if self.name_y == "nothing":
1134            self.field_y = np.array(range(self.dim_y))
1135            self.name_y = "y grid points"
1136          else:
1137            self.field_y = self.f.variables[self.name_y]
1138          # ensure that lon and lat are 2D fields
1139          # 1. simple 1D case (not time-varying)
1140          if len(self.field_x.shape)*len(self.field_y.shape) == 1:
1141               if self.verbose: print "**** OK. recasting lon and lat as 2D fields." 
1142               [self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
1143          # 2. complex 3D case (time-varying, actually just copied over time axis)
1144          elif len(self.field_x.shape)*len(self.field_y.shape) == 9:
1145               if self.verbose: print "**** OK. reducing lon and lat as 2D fields. get rid of time."
1146               self.field_x = self.field_x[0,:,:]
1147               self.field_y = self.field_y[0,:,:]
1148          # if xy axis are apparently undefined, set 2D grid points axis.
1149          if "grid points" not in self.name_x:
1150            if self.field_x.all() == self.field_x[0,0]:
1151               print "!! WARNING !! xy axis look undefined. creating non-dummy ones."
1152               self.field_x = np.array(range(self.dim_x)) ; self.name_x = "x grid points"
1153               self.field_y = np.array(range(self.dim_y)) ; self.name_y = "y grid points"
1154               [self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
1155          if self.dim_x > 1: 
1156               if self.verbose: print "**** OK. x axis %4.0f values [%5.1f,%5.1f]" % (self.dim_x,self.field_x.min(),self.field_x.max())
1157          if self.dim_y > 1: 
1158               if self.verbose: print "**** OK. y axis %4.0f values [%5.1f,%5.1f]" % (self.dim_y,self.field_y.min(),self.field_y.max())
1159          # ALTITUDE. Try preset fields. If not present set grid points axis.
1160          # WARNING: how do we do if several are available?
1161          self.name_z = "nothing"
1162          for c in glob_listz:
1163            if c in self.f.variables.keys():
1164             self.name_z = c
1165          if self.name_z == "nothing":
1166            self.field_z = np.array(range(self.dim_z))
1167            self.name_z = "z grid points"
1168          else:
1169            self.field_z = self.f.variables[self.name_z][:] # specify dimension
1170                                                            # TBD: have to check that this is not a 3D field
1171          if self.dim_z > 1: 
1172               if self.verbose: print "**** OK. z axis %4.0f values [%5.1f,%5.1f]" % (self.dim_z,self.field_z.min(),self.field_z.max())
1173          # TIME. Try preset fields.
1174          self.name_t = "nothing"
1175          for c in glob_listt:
1176            if c in self.f.dimensions.keys():
1177             self.name_t = c
1178          try:
1179            # speed up: only get first value, last one.
1180            dafirst = self.f.variables[self.name_t][0]
1181            dalast = self.f.variables[self.name_t][self.dim_t-1]
1182            self.field_t = np.linspace(dafirst,dalast,num=self.dim_t)
1183            if dafirst == dalast: self.field_t = np.array([dafirst])
1184          except:
1185            # ... or if a problem encountered, define a simple time axis
1186            self.field_t = np.array(range(self.dim_t))
1187            self.name_t = "t grid points"
1188          if self.dim_t > 1: 
1189               if self.verbose: print "**** OK. t axis %4.0f values [%5.1f,%5.1f]" % (self.dim_t,self.field_t.min(),self.field_t.max())     
1190
1191    # change time axis
1192    # ... add your options here!
1193    # --------------------------
1194    def performtimechange(self):
1195        if self.changetime is not None \
1196          and self.name_t != "t grid points":
1197            if self.verbose: print "**** OK. Converting time axis:",self.changetime
1198            if self.changetime == "mars_sol2ls": 
1199                self.field_t = ppcompute.mars_sol2ls(self.field_t)
1200            else:
1201                print "!! WARNING !! This time change is not implemented. Nothing is done."
1202
1203    # get list of index to be retrieved for time axis
1204    ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1205    # -------------------------------
1206    def getindextime(self,dalist=None,ind=None,stride=1):
1207        if self.method_t == "free": 
1208            self.index_t = np.arange(0,self.dim_t,stride)
1209            if self.dim_t > 1: 
1210                self.dimplot = self.dimplot + 1 
1211                if self.verbose: print "**** OK. t values. all."
1212            else:               
1213                self.method_t = "fixed"
1214                if self.verbose: print "**** OK. no t dimension."
1215        elif self.method_t == "comp":
1216            start = np.argmin( np.abs( self.field_t - dalist[ind][0] ) )
1217            stop = np.argmin( np.abs( self.field_t - dalist[ind][1] ) )
1218            self.index_t = np.arange(start,stop,stride)
1219            if self.verbose: print "**** OK. t values. comp over interval ",self.field_t[start],self.field_t[stop]," nvalues=",self.index_t.size
1220        elif self.method_t == "fixed":
1221            self.index_t.append( np.argmin( np.abs( self.field_t - dalist[ind][0] ) ))
1222            if self.verbose: print "**** OK. t values",self.field_t[self.index_t]
1223        else:
1224            print "!! ERROR !! method "+self.method_t+" not supported"
1225        self.index_t = np.array(self.index_t)
1226             
1227    # get list of index to be retrieved for vertical axis
1228    ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1229    # -------------------------------
1230    def getindexvert(self,dalist=None,ind=None,stride=1):
1231        if self.method_z == "free": 
1232            self.index_z = np.arange(0,self.dim_z,stride)
1233            if self.dim_z > 1: 
1234                self.dimplot = self.dimplot + 1
1235                if self.verbose: print "**** OK. z values. all."
1236            else:               
1237                self.method_z = "fixed"
1238                if self.verbose: print "**** OK. no z dimension."
1239        elif self.method_z == "comp":
1240            start = np.argmin( np.abs( self.field_z - dalist[ind][0] ) )
1241            stop = np.argmin( np.abs( self.field_z - dalist[ind][1] ) )
1242            self.index_z = np.arange(start,stop,stride)
1243            if self.verbose: print "**** OK. z values. comp over interval",self.field_z[start],self.field_z[stop]," nvalues=",self.index_z.size
1244        elif self.method_z == "fixed":
1245            self.index_z.append( np.argmin( np.abs( self.field_z - dalist[ind][0] ) ))
1246            if self.verbose: print "**** OK. z values",self.field_z[self.index_z]
1247        else:
1248            if self.verbose: print "!! ERROR !! method "+self.method_z+" not supported"
1249        self.index_z = np.array(self.index_z)
1250
1251    # get list of index to be retrieved for horizontal grid
1252    # --> index_x and index_y are slices to be retrieved from NETCDF files
1253    # --> index_x2d and index_y2d are the actual (x,y) coordinates corresponding to each relevant point
1254    # [this is slightly more complicated because 2D arrays for lat-lon projection possibly irregular]
1255    # NB: to append index we use lists (the most convenient) then we convert into a numpy.array
1256    ### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1257    # -------------------------------
1258    def getindexhori(self,dalistx=None,dalisty=None,indx=None,indy=None,stridex=1,stridey=1):
1259        ## get what is the method over x and y axis
1260        test = self.method_x+self.method_y
1261        ## CASE 0, EASY CASES:
1262        ## - LAT IS FREE (we do here what must be done whatever LON case is)
1263        ## - LON IS FREE (we do here what must be done whatever LAT case is)
1264        ## - LAT IS COMP AND LON IS FREE
1265        ## - LON IS COMP AND LAT IS FREE
1266        if self.method_x == "free" or test in ["compfree","compcomp"]:
1267            self.index_x = range(0,self.dim_x,stridex)
1268            if self.dim_x > 1: 
1269                if self.method_x == "free": self.dimplot = self.dimplot + 1
1270                if self.verbose: print "**** OK. x values. all."
1271            else:               
1272                self.method_x = "fixed"
1273                if self.verbose: print "**** OK. no x dimension."
1274        if self.method_y == "free" or test in ["freecomp","compcomp"]:
1275            self.index_y = range(0,self.dim_y,stridey)
1276            if self.dim_y > 1: 
1277                if self.method_y == "free": self.dimplot = self.dimplot + 1
1278                if self.verbose: print "**** OK. y values. all."
1279            else:               
1280                self.method_y = "fixed"
1281                if self.verbose: print "**** OK. no y dimension."
1282        ## CASE 0 above, this is just for continuity for free.
1283        ## ... for comp we have to select bounds.
1284        ## ... TBD: take the bool array strategy for what follows!
1285        if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
1286            ### ref1_dirty_hack
1287            ### ... for the moment this is a hack. but actually this is more powerful.
1288            if self.method_x == "comp":
1289                yeah = (self.field_x >= dalistx[indx][0])*(self.field_x <= dalistx[indx][1])
1290                self.index_x = yeah[0,:]
1291            if self.method_y == "comp":
1292                yeah = (self.field_y >= dalisty[indy][0]) * (self.field_y <= dalisty[indy][1])
1293                self.index_y = yeah[:,0]
1294            self.index_x2d = self.index_x
1295            self.index_y2d = self.index_y
1296        ## AND NOW THE LITTLE BIT MORE COMPLICATED CASES
1297        ## CASE 1 LAT AND LON ARE FIXED
1298        elif test == "fixedfixed":
1299            idy,idx = np.unravel_index( np.argmin( ( self.field_x - dalistx[indx][0])**2 + (self.field_y - dalisty[indy][0])**2 ), self.field_x.shape ) 
1300                          #TBD: pb with staggered coord
1301            if idx not in self.index_x:  self.index_x.append(idx)
1302            if idy not in self.index_y:  self.index_y.append(idy)
1303            self.index_x2d.append(idx)
1304            self.index_y2d.append(idy)
1305        ## CASE 2 LON IS FIXED BUT NOT LAT
1306        elif test in ["fixedfree","fixedcomp"]:
1307            # find where are requested x values for each y on the free dimension
1308            # NB: this does not work for non-bijective cases e.g. polar stereographic
1309            for iy in range(self.dim_y):
1310              idx = np.argmin( np.abs( self.field_x[iy,:] - dalistx[indx][0] ) )
1311              # if comp is requested we select only indexes which yield values between requested min and max
1312              storeval = (self.method_y == "comp") and (self.field_y[iy,idx] > dalisty[indy][0]) and (self.field_y[iy,idx] < dalisty[indy][1])
1313              storeval = storeval or (self.method_y == "free")
1314              if storeval:
1315                  if idx not in self.index_x:  self.index_x.append(idx)
1316                  if iy not in self.index_y and self.method_y == "comp": self.index_y.append(iy)
1317                  if idx not in self.index_x2d or iy not in self.index_y2d:
1318                    self.index_x2d.append(idx)
1319                    self.index_y2d.append(iy)
1320        ## CASE 3 LAT IS FIXED BUT NOT LON
1321        elif test in ["freefixed","compfixed"]:
1322            # find where are requested y values for each x on the free dimension
1323            # NB: this does not work for non-bijective cases e.g. polar stereographic
1324            for ix in range(self.dim_x):
1325              idy = np.argmin( np.abs( self.field_y[:,ix] - dalisty[indy][0] ) )
1326              # if comp is requested we select only indexes which yield values between requested min and max
1327              storeval = (self.method_x == "comp") and (self.field_x[idy,ix] > dalistx[indx][0]) and (self.field_x[idy,ix] < dalistx[indx][1])
1328              storeval = storeval or (self.method_x == "free")
1329              if storeval:
1330                  if idy not in self.index_y:  self.index_y.append(idy)
1331                  if ix not in self.index_x and self.method_x == "comp": self.index_x.append(ix)
1332                  if ix not in self.index_x2d or idy not in self.index_y2d:
1333                    self.index_x2d.append(ix)
1334                    self.index_y2d.append(idy)
1335        ## check index tab
1336        if len(self.index_x) == 0 or len(self.index_y) == 0:
1337            print "!! ERROR !! no indices found. check prescribed latitudes or longitudes" ; exit()
1338        ## ensure the array is a numpy array for getfield to work
1339        self.index_x = np.array(self.index_x)
1340        self.index_y = np.array(self.index_y)
1341        self.index_x2d = np.array(self.index_x2d)
1342        self.index_y2d = np.array(self.index_y2d)
1343        ### print extrema
1344        printx = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
1345        printy = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
1346        if self.verbose: 
1347            print "**** OK. x values (min,max).", printx.min(),printx.max()
1348            print "**** OK. y values (min,max).", printy.min(),printy.max()
1349
1350    # get the field from the NETCDF file and perform averages
1351    # -------------------------------
1352    def getfield(self):
1353        ## first tell what is to be done
1354        if self.dimplot > 2:                       print "**** !! ERROR !! "+str(self.dimplot)+"D plots not supported!" ; exit()
1355        elif self.dimplot == 0 and self.verbose:   print "**** OK. 0D value requested."
1356        elif self.dimplot == 1 and self.verbose:   print "**** OK. 1D plot requested."
1357        elif self.verbose:                         print "**** OK. 2D section requested."
1358        # well, now get field from netcdf file
1359        # part below is necessary otherwise there is an index error below
1360        if self.index_x.size == 1: self.index_x = self.index_x[0]
1361        if self.index_y.size == 1: self.index_y = self.index_y[0]
1362        if self.index_z.size == 1: self.index_z = self.index_z[0]
1363        if self.index_t.size == 1: self.index_t = self.index_t[0]
1364        # then retrieve what is requested by user
1365        # each self.dim case corresponds to tests in the beginning of getdim.
1366        time0 = timelib.time()
1367        if self.verbose: print "**** OK. I am getting values from files. Please wait."
1368        if self.dim == 1: 
1369            nt = self.index_t.size ; nz = 1 ; ny = 1 ; nx = 1
1370            self.field = self.f.variables[self.var][self.index_t]
1371        elif self.dim == 2:
1372            nt = 1 ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
1373            self.field = self.f.variables[self.var][self.index_y,self.index_x]
1374        elif self.dim == 3:
1375            nt = self.index_t.size ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
1376            self.field = self.f.variables[self.var][self.index_t,self.index_y,self.index_x]
1377            # this is far faster than retrieving each term with a loop
1378        elif self.dim == 4:
1379            nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = self.index_x.size
1380            self.field = self.f.variables[self.var][self.index_t,self.index_z,self.index_y,self.index_x]
1381        else:
1382            print "!! ERROR !! field would have more than four dimensions ?" ; exit()
1383        # dirty hack (AS) ref1_dirty_hack
1384        # waiting for more fundamental modifications. case when self.index_y is a bool array.
1385        # ... be careful if no point...
1386        try:
1387            if type(self.index_x[0]) == np.bool_: nx = np.sum(self.index_x) ## gives the size of the True part!
1388            if type(self.index_y[0]) == np.bool_: ny = np.sum(self.index_y) ## gives the size of the True part!
1389        except:
1390            pass
1391        # NB: ... always 4D array but possibly with "size 1" dimensions
1392        #     ... if one dimension is missing because 1D 2D or 3D requests, make it appear again
1393        self.field = np.reshape(self.field,(nt,nz,ny,nx))
1394        if self.verbose: print "**** OK. I got %7.1e values. This took me %6.4f seconds" % (nx*ny*nz*nt,timelib.time() - time0)
1395        if self.verbose: print "**** OK. I got var "+self.var+" with shape",self.field.shape
1396        # reduce coordinates to useful points
1397        # ... TBD: this should be ordered in the case of non-regular projections
1398        if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
1399          # we need 2D coordinates (free) or we get broadcast problem (comp) so we use np.ix
1400          self.field_x = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
1401          self.field_y = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
1402        else:
1403          # we are OK with 1D coordinates
1404          self.field_x = self.field_x[self.index_y2d, self.index_x2d]
1405          self.field_y = self.field_y[self.index_y2d, self.index_x2d]
1406        self.field_z = self.field_z[self.index_z]
1407        self.field_t = self.field_t[self.index_t]
1408        # extract relevant horizontal points
1409        # TBD: is compfree OK with computing on irregular grid?
1410        test = self.method_x + self.method_y
1411        if test in ["fixedfixed","freefree"]:
1412          pass
1413        elif test in ["fixedfree","fixedcomp"] or test in ["freefixed","compfixed"]: 
1414          time0 = timelib.time()
1415          # now have to obtain the new indexes which correspond to the extracted self.field
1416          # for it to work with unique index, ensure that any index_* is a numpy array
1417          if not isinstance(self.index_x, np.ndarray): self.index_x = np.array([self.index_x])
1418          if not isinstance(self.index_y, np.ndarray): self.index_y = np.array([self.index_y])
1419          if not isinstance(self.index_z, np.ndarray): self.index_z = np.array([self.index_z])
1420          if not isinstance(self.index_t, np.ndarray): self.index_t = np.array([self.index_t])
1421          for val in self.index_x: self.index_x2d[np.where(self.index_x2d == val)] = np.where(self.index_x == val)[0]
1422          for val in self.index_y: self.index_y2d[np.where(self.index_y2d == val)] = np.where(self.index_y == val)[0]
1423               ##### VERY EXPENSIVE
1424               ## recast self.field with 2D horizontal arrays because we might have extracted
1425               ## more than what is to be actually plot or computed, in particular for comps on 2D lat,lon coordinates
1426               #self.field = self.field[np.ix_(self.index_t,self.index_z,self.index_y2d,self.index_x2d)]
1427               #(nt,nz,ny,nx) = self.field.shape       
1428          # prepare the loop on all relevant horizontal points
1429          if self.method_x in ["comp","free"]:   
1430              nnn = self.index_x2d.shape[0] ; what_I_am_supposed_to_do = "keepx"
1431          elif self.method_y in ["comp","free"]: 
1432              nnn = self.index_y2d.shape[0] ; what_I_am_supposed_to_do = "keepy" 
1433          # LOOP to extract only useful values over horizontal dimensions
1434          # only take diagonal terms, do not loop on all self.index_x2d*self.index_y2d
1435          # ... this method is fast enough, perhaps there is a faster way though
1436          # ... (for sure the method with np.diag is much slower)
1437          for iii in range(nnn):
1438           ix = self.index_x2d[iii] ; iy = self.index_y2d[iii]
1439           for iz in range(self.index_z.size):
1440            for it in range(self.index_t.size):
1441              if what_I_am_supposed_to_do == "keepx":    self.field[it,iz,0,ix] = self.field[it,iz,iy,ix]
1442              elif what_I_am_supposed_to_do == "keepy":  self.field[it,iz,iy,0] = self.field[it,iz,iy,ix]
1443          if self.verbose: print "**** OK. I got to pick the right values for your request. This took me %6.4f seconds" % (timelib.time() - time0)
1444          # we only keep the one value that was modified on the dimension which is not free
1445          if what_I_am_supposed_to_do == "keepx":     self.field = self.field[:,:,0,:] ; ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1446          elif what_I_am_supposed_to_do == "keepy":   self.field = self.field[:,:,:,0] ; nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1447        # make a mask in case there are non-NaN missing values. (what about NaN missing values?)
1448        # ... this is important for computations below (see ppcompute)
1449        masked = np.ma.masked_where(np.abs(self.field) > 1e25,self.field)
1450        if masked.mask.any() == True:
1451             if self.verbose: print "!! WARNING !! Values over +-1e25 are considered missing values."
1452             self.field = masked
1453             self.field.set_fill_value([np.NaN])
1454
1455    # perform computations
1456    # -------------------------------
1457    # available: mean, max, min, meanarea
1458    # TB: integrals? for derivatives, define a function self.dx()
1459    def computations(self): 
1460        nt,nz,ny,nx = self.field.shape
1461        # treat the case of mean on fields normalized with grid mesh area
1462        # ... this is done in the .area() method.
1463        # after that self.field contains field*area/totarea
1464        if "area" in self.compute: 
1465            if "comp" in self.method_x+self.method_y: 
1466                self.area()
1467            else:
1468                if self.verbose: print "!! WARNING !! No area accounted for (computing on t and/or z axis)."
1469        # now ready to compute [TBD: we would like to have e.g. mean over x,y and min over t ??]
1470        if self.method_t == "comp":
1471            if self.verbose: print "**** OK. Computing over t axis."
1472            if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=0)
1473            elif self.compute == "min": self.field = ppcompute.min(self.field,axis=0)
1474            elif self.compute == "max": self.field = ppcompute.max(self.field,axis=0)
1475            else: print "!! ERROR !! operation not supported." ; exit()
1476            nt = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1477        if self.method_z == "comp": 
1478            if self.verbose: print "**** OK. Computing over z axis."
1479            if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=1)
1480            elif self.compute == "min": self.field = ppcompute.min(self.field,axis=1)
1481            elif self.compute == "max": self.field = ppcompute.max(self.field,axis=1)
1482            else: print "!! ERROR !! operation not supported." ; exit()
1483            nz = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1484        if self.method_y == "comp": 
1485            if self.verbose: print "**** OK. Computing over y axis."
1486            if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=2)
1487            elif self.compute == "min": self.field = ppcompute.min(self.field,axis=2)
1488            elif self.compute == "max": self.field = ppcompute.max(self.field,axis=2)
1489            elif self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=2)
1490            else: print "!! ERROR !! operation not supported." ; exit()
1491            ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1492            if self.field_x.ndim == 2: self.field_x = self.field_x[0,:] # TBD: this is OK for regular grid but not for irregular
1493        if self.method_x == "comp":
1494            if self.verbose: print "**** OK. Computing over x axis."
1495            if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=3)
1496            elif self.compute == "min": self.field = ppcompute.min(self.field,axis=3)
1497            elif self.compute == "max": self.field = ppcompute.max(self.field,axis=3)
1498            elif self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=3)
1499            else: print "!! ERROR !! operation not supported." ; exit()
1500            nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1501            if self.field_y.ndim == 2: self.field_y = self.field_y[:,0] # TBD: this is OK for regular grid but not for irregular
1502        # remove all dimensions with size 1 to prepare plot (and check the resulting dimension with dimplot)
1503        self.field = np.squeeze(self.field)
1504        if self.field.ndim != self.dimplot: 
1505            print "!! ERROR !! Problem: self.field is different than plot dimensions", self.field.ndim, self.dimplot ; exit()
1506        if self.verbose: 
1507            print "**** OK. Final shape for "+self.var+" after averaging and squeezing",self.field.shape
1508   
1509    # get areas for computations and ponderate self.field by area/totarea
1510    # -------------------------------------------------------------------
1511    def area(self):
1512        if self.verbose: print "**** OK. Get area array for computations."
1513        # create a request object for area
1514        # ... and copy known attributes from self
1515        aire = onerequest()
1516        aire.copy(self)
1517        # get area field name
1518        aire.var = "nothing"
1519        for c in glob_listarea:
1520         if c in aire.f.variables.keys():
1521            aire.var = c
1522        # do not try to calculate areas automatically
1523        if aire.var == "nothing":
1524            print "!! ERROR !! area variable not found... needs to be added in set_ppclass.txt?"
1525            exit()
1526        # define area request dimensions
1527        aire.getdim()
1528        # ensure this is a 2D horizontal request and define indexes
1529        #    ... areas are not supposed to vary with time and height
1530        aire.method_x = "free" ; aire.method_y = "free"
1531        aire.getindexhori() ; aire.dimplot = 2
1532        aire.method_z = "fixed" ; aire.field_z = np.array([0]) ; aire.index_z = np.array([0])
1533        aire.method_t = "fixed" ; aire.field_t = np.array([0]) ; aire.index_t = np.array([0])
1534        # read the 2D area array in netCDF file
1535        aire.getfield()
1536        aire.field = np.squeeze(aire.field)
1537        # reduce with self horizontal indexes
1538        if "fixed" in self.method_x+self.method_y:
1539            aire.field = aire.field[self.index_y,self.index_x]
1540        # calculate total area
1541        # ... 2D comp is easy. 1D comp is a bit less easy but simple array manipulation.
1542        if "free" in self.method_x+self.method_y:
1543            if self.method_x == "free":
1544                totarea = ppcompute.sum(aire.field,axis=0)
1545                totarea = np.reshape(totarea,(1,totarea.size))
1546                totarea = np.tile(totarea,(1,self.index_x))
1547            elif self.method_y == "free":
1548                totarea = ppcompute.sum(aire.field,axis=1)
1549                totarea = np.reshape(totarea,(totarea.size,1))
1550                totarea = np.tile(totarea,(1,self.index_x.size))
1551        elif self.method_x == "comp" and self.method_y == "comp":
1552            aire.field = aire.field[np.ix_(self.index_y, self.index_x)] # reduce to requested indexes only
1553            totarea = ppcompute.sum(ppcompute.sum(aire.field,axis=1),axis=0)
1554        else:
1555            if self.verbose: print "!! WARNING !! Not account for areas. Only averaging over z and/or t axis."
1556        # normalize by total area
1557        print "**** OK. I can now normalize field by areas."
1558        aire.field = aire.field / totarea
1559        # tile area array over self t and z axis so that area field could be multiplied with self.field
1560        aire.field = np.tile(aire.field,(self.index_t.size,self.index_z.size,1,1))
1561        if self.field.shape != aire.field.shape:
1562            print "!! ERROR !! Problem in area(). Check array shapes."
1563            print "Field vs. aire:",self.field.shape,aire.field.shape ; exit()
1564        else:
1565            self.field = self.field*aire.field
1566
1567    # define coordinates for plot
1568    # -------------------------------
1569    def definecoord(self):
1570        I_got_abs = False ; I_got_ord = False
1571        # here is the thing. time is usually taken as an abscissa so we start with time.
1572        if self.method_t ==  "free": 
1573            self.absc = self.field_t ; self.absclab = self.name_t
1574            I_got_abs = True
1575        # then we usually have x as an abscissa.
1576        if self.method_x == "free":
1577            if I_got_abs: 
1578                self.ordi = self.field_x ; self.ordilab = self.name_x
1579                I_got_ord = True
1580            else:         
1581                self.absc = self.field_x ; self.absclab = self.name_x
1582                I_got_abs = True
1583        # ... or we have y
1584        if self.method_y == "free":
1585            if I_got_abs:   
1586                self.ordi = self.field_y ; self.ordilab = self.name_y
1587                I_got_ord = True
1588            else:         
1589                self.absc = self.field_y ; self.absclab = self.name_y
1590                I_got_abs = True
1591        # ... and we end with z because it is usually not an abscissa (profiles).
1592        if self.method_z == "free":
1593            if self.field_z[0] > self.field_z[1]:
1594                self.invert_axes = True # the axis will be turned upside-down
1595            if I_got_abs: 
1596                self.ordi = self.field_z ; self.ordilab = self.name_z
1597                I_got_ord = True
1598            else:
1599                self.absc = self.field_z ; self.absclab = self.name_z
1600                I_got_abs = True
1601                self.swap_axes = True # says that altitude is not supposed to remain as an abscissa
1602        if I_got_abs and self.verbose: print "**** OK. abscissa:",self.absclab, self.absc.shape
1603        if I_got_ord and self.verbose: print "**** OK. ordinate:",self.ordilab, self.ordi.shape
Note: See TracBrowser for help on using the repository browser.