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

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

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