source: trunk/UTIL/PYTHON/mcd/mcd.py @ 1076

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

UTIL PYTHON. MCD online. previous changes after progress meeting were not committed.
  • Property svn:executable set to *
File size: 46.6 KB
Line 
1####################################################
2### A Python Class for the Mars Climate Database ###
3### ---------------------------------------------###
4### Aymeric SPIGA 17-21/04/2012                  ###
5### ---------------------------------------------###
6### (see mcdtest.py for examples of use)         ###
7####################################################
8
9import numpy as np
10import matplotlib.pyplot as mpl
11import frozen_myplot as myplot
12
13
14class mcd():
15 
16    def __repr__(self):
17    # print out a help string when help is invoked on the object
18        whatprint = 'MCD object. \"help(mcd)\" for more information\n'
19        return whatprint
20
21########################
22### Default settings ###
23########################
24
25    def __init__(self):
26    # default settings
27        ## 0. general stuff
28        self.name      = "MCD v4.3"
29        self.ack       = "Mars Climate Database (c) LMD/OU/IAA/ESA/CNES"
30        #self.dset      = '/home/aymeric/Science/MCD_v4.3/data/'
31        self.dset      = '/home/marshttp/MCD_v4.3/data/'
32        ## 1. spatio-temporal coordinates
33        self.lat       = 0.
34        self.lats      = None
35        self.late      = None
36        self.lon       = 0.
37        self.lons      = None
38        self.lone      = None
39        self.loct      = 0.
40        self.locts     = None
41        self.locte     = None
42        self.xdate     = 0.  # see datekey
43        self.xdates    = None
44        self.xdatee    = None
45        self.xz        = 10. # see zkey
46        self.xzs       = None
47        self.xze       = None
48        ## 1bis. related settings
49        self.zkey      = 3  # specify that xz is the altitude above surface (m)
50                            # zkey  : <integer>   type of vertical coordinate xz
51                            # 1 = radius from centre of planet (m)
52                            # 2 = height above areoid (m) (MOLA zero datum)
53                            # 3 = height above surface (m)
54                            # 4 = pressure level (Pa)
55                            # 5 = altitude above mean Mars Radius(=3396000m) (m)
56        self.datekey   = 1  # 0 = "Earth time": xdate is given in Julian days (localtime must be set to zero)
57                            # 1 = "Mars date": xdate is the value of Ls
58        ## 2. climatological options
59        self.dust      = 2  #our best guess MY24 scenario, with solar average conditions
60        self.hrkey     = 1  #set high resolution mode on (hrkey=0 to set high resolution off)
61        ## 3. additional settings for advanced use
62        self.extvarkey = 1  #extra output variables (1: yes, 0: no)
63        self.perturkey = 0  #integer perturkey ! perturbation type (0: none)
64        self.seedin    = 1  #random number generator seed (unused if perturkey=0)
65        self.gwlength  = 0. #gravity Wave wavelength (unused if perturkey=0)
66        ## outputs. just to define attributes.
67        ## --> in update
68        self.pres = None ; self.dens = None ; self.temp = None ; self.zonwind = None ; self.merwind = None ; self.meanvar = None ; self.extvar = None
69        self.seedout = None ; self.ierr = None
70        ## --> in prepare
71        self.xcoord = None ; self.ycoord = None
72        self.prestab = None ; self.denstab = None ; self.temptab = None 
73        self.zonwindtab = None ; self.merwindtab = None ; self.meanvartab = None ; self.extvartab = None
74        ## plot stuff
75        self.xlabel = None ; self.ylabel = None ; self.title = ""
76        self.vertplot = False
77        self.fmt = "%.2e" 
78        self.colorm = "jet"
79        self.fixedlt = False
80        self.zonmean = False
81        self.min2d = None
82        self.max2d = None
83        self.dpi = 80.
84        self.islog = False
85
86    def toversion5(self):
87        self.name      = "MCD v5.0"
88        self.dset      = '/home/marshttp/MCD_v5.0/data/'
89        self.extvarkey = np.ones(100)
90
91    def viking1(self): self.name = "Viking 1 site. MCD v4.3 output" ; self.lat = 22.48 ; self.lon = -49.97 ; self.xdate = 97.
92    def viking2(self): self.name = "Viking 2 site. MCD v4.3 output" ; self.lat = 47.97 ; self.lon = -225.74 ; self.xdate = 117.6
93
94    def getdustlabel(self):
95        if self.dust == 1: 
96            self.dustlabel = "MY24 minimum solar scenario"
97            if "v5" in self.name: self.dustlabel = "climatology average solar scenario"
98        elif self.dust == 2: 
99            self.dustlabel = "MY24 average solar scenario"
100            if "v5" in self.name: self.dustlabel = "climatology minimum solar scenario"
101        elif self.dust == 3: 
102            self.dustlabel = "MY24 maximum solar scenario"
103            if "v5" in self.name: self.dustlabel = "climatology maximum solar scenario"
104        elif self.dust == 4: self.dustlabel = "dust storm minimum solar scenario"
105        elif self.dust == 5: self.dustlabel = "dust storm average solar scenario"
106        elif self.dust == 6: self.dustlabel = "dust storm maximum solar scenario"
107        elif self.dust == 7: self.dustlabel = "warm scenario (dusty, maximum solar)"
108        elif self.dust == 8: self.dustlabel = "cold scenario (low dust, minimum solar)"
109
110    def gettitle(self,oneline=False):
111        self.getdustlabel()
112        self.title = self.name + " with " + self.dustlabel + "."
113        if self.datekey == 1:    self.title = self.title + " Ls " + str(self.xdate) + "deg."
114        elif self.datekey == 0:  self.title = self.title + " JD " + str(self.xdate) + "."
115        if not oneline: self.title = self.title + "\n"
116        if self.lats is None:  self.title = self.title + " Latitude " + str(self.lat) + "N"
117        if self.zonmean and self.lats is not None and self.xzs is not None: 
118            self.title = self.title + "Zonal mean over all longitudes."
119        elif self.lons is None: 
120            self.title = self.title + " Longitude " + str(self.lon) + "E"
121        if self.xzs is None:   
122            self.vertunits()
123            self.title = self.title + " Altitude " + str(self.xz) + " " + self.vunits
124        if self.locts is None:
125            self.title = self.title + " Local time " + str(self.loct) + "h"
126            if not self.fixedlt:  self.title = self.title + " (at longitude 0) "
127
128    def getextvarlab(self,num):
129        whichfield = { \
130        91: "Pressure (Pa)", \
131        92: "Density (kg/m3)", \
132        93: "Temperature (K)", \
133        94: "W-E wind component (m/s)", \
134        95: "S-N wind component (m/s)", \
135        1: "Radial distance from planet center (m)",\
136        2: "Altitude above areoid (Mars geoid) (m)",\
137        3: "Altitude above local surface (m)",\
138        4: "orographic height (m) (surf alt above areoid)",\
139        5: "Ls, solar longitude of Mars (deg)",\
140        6: "LST local true solar time (hrs)",\
141        7: "Universal solar time (LST at lon=0) (hrs)",\
142        8: "Air heat capacity Cp (J kg-1 K-1)",\
143        9: "gamma=Cp/Cv Ratio of specific heats",\
144        10: "density RMS day to day variations (kg/m^3)",\
145        11: "[not defined]",\
146        12: "[not defined]",\
147        13: "scale height H(p) (m)",\
148        14: "GCM orography (m)",\
149        15: "surface temperature (K)",\
150        16: "daily max mean surface temperature (K)",\
151        17: "daily min mean surface temperature (K)",\
152        18: "surf. temperature RMS day to day variations (K)",\
153        19: "surface pressure (Pa)",\
154        20: "GCM surface pressure (Pa)",\
155        21: "atmospheric pressure RMS day to day variations (Pa)",\
156        22: "surface pressure RMS day to day variations (Pa)",\
157        23: "temperature RMS day to day variations (K)",\
158        24: "zonal wind RMS day to day variations (m/s)",\
159        25: "meridional wind RMS day to day variations (m/s)",\
160        26: "vertical wind component (m/s) >0 when downwards!",\
161        27: "vertical wind RMS day to day variations (m/s)",\
162        28: "small scale perturbation (gravity wave) (kg/m^3)",\
163        29: "q2: turbulent kinetic energy (m2/s2)",\
164        30: "[not defined]",\
165        31: "thermal IR flux to surface (W/m2)",\
166        32: "solar flux to surface (W/m2)",\
167        33: "thermal IR flux to space (W/m2)",\
168        34: "solar flux reflected to space (W/m2)",\
169        35: "surface CO2 ice layer (kg/m2)",\
170        36: "DOD: Dust column visible optical depth",\
171        37: "Dust mass mixing ratio (kg/kg)",\
172        38: "DOD RMS day to day variations",\
173        39: "DOD total standard deviation over season",\
174        40: "Water vapor column (kg/m2)",\
175        41: "Water vapor vol. mixing ratio (mol/mol)",\
176        42: "Water ice column (kg/m2)",\
177        43: "Water ice mixing ratio (mol/mol)",\
178        44: "O3 ozone vol. mixing ratio (mol/mol)",\
179        45: "[CO2] vol. mixing ratio (mol/mol)",\
180        46: "[O] vol. mixing ratio (mol/mol)",\
181        47: "[N2] vol. mixing ratio (mol/mol)",\
182        48: "[CO] vol. mixing ratio (mol/mol)",\
183        49: "R: Molecular gas constant (J K-1 kg-1)",\
184        50: "Air viscosity estimation (N s m-2)"
185        }
186        ### MCD version 5 new variables. AS 12/2012.
187        if "v5" in self.name:
188          whichfield[29] = "not used (set to zero)"
189          whichfield[30] = "Surface roughness length z0 (m)"
190          whichfield[37] = "DOD RMS day to day variations"
191          whichfield[38] = "Dust mass mixing ratio (kg/kg)"
192          whichfield[39] = "Dust effective radius (m)"
193          whichfield[44] =  whichfield[43]
194          whichfield[43] =  whichfield[42]
195          whichfield[42] =  whichfield[41]
196          whichfield[41] =  whichfield[40]
197          whichfield[40] = "Dust deposition on flat surface (kg m-2 s-1)"
198          whichfield[45] = "Water ice effective radius (m)"
199          whichfield[46] = "Convective PBL height (m)"
200          whichfield[47] = "Max. upward convective wind within the PBL (m/s)"
201          whichfield[48] = "Max. downward convective wind within the PBL (m/s)"
202          whichfield[49] = "Convective vertical wind variance at level z (m2/s2)"
203          whichfield[50] = "Convective eddy vertical heat flux at level z (m/s/K)"
204          whichfield[51] = "Surface wind stress (Kg/m/s2)"
205          whichfield[52] = "Surface sensible heat flux (W/m2) (<0 when flux from surf to atm.)"
206          whichfield[53] = "R: Molecular gas constant (J K-1 kg-1)"
207          whichfield[54] = "Air viscosity estimation (N s m-2)"
208          whichfield[55] = "not used (set to zero)"
209          whichfield[56] = "not used (set to zero)"
210          whichfield[57] = "[CO2] vol. mixing ratio (mol/mol)"
211          whichfield[58] = "[N2] vol. mixing ratio (mol/mol)"
212          whichfield[59] = "[Ar] vol. mixing ratio (mol/mol)"
213          whichfield[60] = "[CO] vol. mixing ratio (mol/mol)"
214          whichfield[61] = "[O] vol. mixing ratio (mol/mol)"
215          whichfield[62] = "[O2] vol. mixing ratio (mol/mol)"
216          whichfield[63] = "[O3] vol. mixing ratio (mol/mol)"
217          whichfield[64] = "[H] vol. mixing ratio (mol/mol)"
218          whichfield[65] = "[H2] vol. mixing ratio (mol/mol)"
219          whichfield[66] = "[electron] vol. mixing ratio (mol/mol)"
220          whichfield[67] = "CO2 column (kg/m2)"
221          whichfield[68] = "N2 column (kg/m2)"
222          whichfield[69] = "Ar column (kg/m2)"
223          whichfield[70] = "CO column (kg/m2)"
224          whichfield[71] = "O column (kg/m2)"
225          whichfield[72] = "O2 column (kg/m2)"
226          whichfield[73] = "O3 column (kg/m2)"
227          whichfield[74] = "H column (kg/m2)"
228          whichfield[75] = "H2 column (kg/m2)"
229          whichfield[76] = "electron column (kg/m2)"
230        if num not in whichfield: myplot.errormess("Incorrect subscript in extvar.")
231        dastuff = whichfield[num]
232        if "(K)" in dastuff:      self.fmt="%.0f"
233        elif "effective radius" in dastuff: self.fmt="%.2e"
234        elif "(Pa)" in dastuff:   self.fmt="%.2e"
235        elif "(W/m2)" in dastuff: self.fmt="%.0f"
236        elif "(m/s)" in dastuff:  self.fmt="%.1f"
237        elif "(m)" in dastuff:    self.fmt="%.0f"
238        else:                     self.fmt="%.2e"
239        return dastuff
240
241    def convertlab(self,num):       
242        ## a conversion from text inquiries to extvar numbers. to be completed.
243        if num == "p": num = 91
244        elif num == "rho": num = 92
245        elif num == "t": num = 93
246        elif num == "u": num = 94
247        elif num == "v": num = 95
248        elif num == "tsurf": num = 15
249        elif num == "topo": num = 4
250        elif num == "h": num = 13
251        elif num == "ps": num = 19
252        elif num == "tau": num = 36
253        elif num == "mtot": 
254            if "v5" in self.name:  num = 41 
255            else:                  num = 40
256        elif num == "icetot": 
257            if "v5" in self.name:  num = 43
258            else:                  num = 42
259        elif num == "h2ovap": 
260            if "v5" in self.name:  num = 42
261            else:                  num = 41
262        elif num == "h2oice": 
263            if "v5" in self.name:  num = 44
264            else:                  num = 43
265        elif num == "cp": num = 8
266        elif num == "rho_ddv": num = 10
267        elif num == "ps_ddv": num = 22
268        elif num == "p_ddv": num = 21
269        elif num == "t_ddv": num = 23
270        elif num == "w": num = 26
271        elif num == "tsurfmx": num = 16
272        elif num == "tsurfmn": num = 17
273        elif num == "lwdown": num = 31
274        elif num == "swdown": num = 32
275        elif num == "lwup": num = 33
276        elif num == "swup": num = 34
277        elif num == "tau": num = 36
278        elif num == "tau_ddv":
279            if "v5" in self.name:  num = 37
280            else:                  num = 38
281        elif num == "qdust":
282            if "v5" in self.name:  num = 38
283            else:                  num = 37
284        elif num == "co2":
285            if "v5" in self.name:  num = 57
286            else:                  num = 45
287        elif num == "o3": 
288            if "v5" in self.name:  num = 63
289            else:                  num = 44
290        elif num == "o": 
291            if "v5" in self.name:  num = 61
292            else:                  num = 46
293        elif num == "co": 
294            if "v5" in self.name:  num = 60
295            else:                  num = 48
296        elif num == "visc": 
297            if "v5" in self.name:  num = 54
298            else:                  num = 50
299        elif num == "co2ice": num = 35
300        elif num == "rdust":
301            if "v5" in self.name:  num = 39
302            else:                  num = 30 # an undefined variable to avoid misleading output
303        elif num == "sdust":
304            if "v5" in self.name:  num = 40
305            else:                  num = 30 # an undefined variable to avoid misleading output
306        elif num == "pbl":
307            if "v5" in self.name:  num = 46
308            else:                  num = 30 # an undefined variable to avoid misleading output
309        elif num == "updraft":
310            if "v5" in self.name:  num = 47
311            else:                  num = 30 # an undefined variable to avoid misleading output
312        elif num == "downdraft":
313            if "v5" in self.name:  num = 48
314            else:                  num = 30 # an undefined variable to avoid misleading output
315        elif num == "pblwvar":
316            if "v5" in self.name:  num = 49
317            else:                  num = 30 # an undefined variable to avoid misleading output
318        elif num == "pblhvar":
319            if "v5" in self.name:  num = 50
320            else:                  num = 30 # an undefined variable to avoid misleading output
321        elif num == "stress":
322            if "v5" in self.name:  num = 51
323            else:                  num = 30 # an undefined variable to avoid misleading output
324        elif num == "ar":
325            if "v5" in self.name:  num = 59
326            else:                  num = 30 # an undefined variable to avoid misleading output
327        elif num == "o2":
328            if "v5" in self.name:  num = 62
329            else:                  num = 30 # an undefined variable to avoid misleading output
330        elif num == "co2col":
331            if "v5" in self.name:  num = 67
332            else:                  num = 30 # an undefined variable to avoid misleading output
333        elif num == "arcol":
334            if "v5" in self.name:  num = 69
335            else:                  num = 30 # an undefined variable to avoid misleading output
336        elif num == "cocol":
337            if "v5" in self.name:  num = 70
338            else:                  num = 30 # an undefined variable to avoid misleading output
339        elif num == "o3col":
340            if "v5" in self.name:  num = 73
341            else:                  num = 30 # an undefined variable to avoid misleading output
342        elif num == "hydro":
343            if "v5" in self.name:  num = 64
344            else:                  num = 30 # an undefined variable to avoid misleading output
345        elif num == "hydro2":
346            if "v5" in self.name:  num = 65
347            else:                  num = 30 # an undefined variable to avoid misleading output
348        elif num == "e":
349            if "v5" in self.name:  num = 66
350            else:                  num = 30 # an undefined variable to avoid misleading output
351        elif num == "ecol":
352            if "v5" in self.name:  num = 76
353            else:                  num = 30 # an undefined variable to avoid misleading output
354        elif not isinstance(num, np.int): myplot.errormess("field reference not found.")
355        return num
356
357###################
358### One request ###
359###################
360
361    def update(self):
362    # retrieve fields from MCD (call_mcd). more info in fmcd.call_mcd.__doc__
363        ## sanity first
364        self.loct = abs(self.loct)%24
365        if self.locts is not None and self.locte is not None: 
366            self.locts = abs(self.locts)%24
367            self.locte = abs(self.locte)%24
368            if self.locts == self.locte: self.locte = self.locts + 24
369        ## now MCD request
370        if "v5" in self.name: from fmcd5 import call_mcd
371        else:                 from fmcd import call_mcd
372        (self.pres, self.dens, self.temp, self.zonwind, self.merwind, \
373         self.meanvar, self.extvar, self.seedout, self.ierr) \
374         = \
375         call_mcd(self.zkey,self.xz,self.lon,self.lat,self.hrkey, \
376             self.datekey,self.xdate,self.loct,self.dset,self.dust, \
377             self.perturkey,self.seedin,self.gwlength,self.extvarkey )
378        ## we use the end of extvar (unused) to store meanvar. this is convenient for getextvar(lab)
379        self.extvar[90] = self.pres ; self.extvar[91] = self.dens
380        self.extvar[92] = self.temp ; self.extvar[93] = self.zonwind ; self.extvar[94] = self.merwind
381        ## treat missing values
382        if self.temp == -999: self.extvar[:] = np.NaN ; self.meanvar[:] = np.NaN
383
384    def printset(self):
385    # print main settings
386        print "zkey",self.zkey,"xz",self.xz,"lon",self.lon,"lat",self.lat,"hrkey",self.hrkey, \
387              "xdate",self.xdate,"loct",self.loct,"dust",self.dust
388
389    def getnameset(self):
390    # set a name referring to settings [convenient for databases]
391        strlat = str(self.lat)+str(self.lats)+str(self.late)
392        strlon = str(self.lon)+str(self.lons)+str(self.lone)
393        strxz = str(self.xz)+str(self.xzs)+str(self.xze)
394        strloct = str(self.loct)+str(self.locts)+str(self.locte)
395        name = str(self.zkey)+strxz+strlon+strlat+str(self.hrkey)+str(self.datekey)+str(self.xdate)+strloct+str(self.dust)
396        if "v5" in self.name: name = "v5_" + name
397        return name
398
399    def printcoord(self):
400    # print requested space-time coordinates
401        print "LAT",self.lat,"LON",self.lon,"LOCT",self.loct,"XDATE",self.xdate
402
403    def printmeanvar(self):
404    # print mean MCD variables
405        print "Pressure = %5.3f pascals. " % (self.pres)
406        print "Density = %5.3f kilograms per cubic meter. " % (self.dens)
407        print "Temperature = %3.0f kelvins (%4.0f degrees celsius)." % (self.temp,self.temp-273.15)
408        print "Zonal wind = %5.3f meters per second." % (self.zonwind)
409        print "Meridional wind = %5.3f meters per second." % (self.merwind)
410        print "Total horizontal wind = %5.3f meters per second." % ( np.sqrt(self.zonwind**2 + self.merwind**2) )
411
412    def printextvar(self,num):
413    # print extra MCD variables
414        num = self.convertlab(num)
415        dastr = str(self.extvar[num-1])
416        if dastr == "nan":   print "!!!! There is a problem, probably a value is requested below the surface !!!!"
417        else:                print self.getextvarlab(num) + " ..... " + dastr
418
419    def printallextvar(self):
420    # print all extra MCD variables   
421        if "v5" in self.name:  limit=76
422        else:                  limit=50
423        for i in range(limit): self.printextvar(i+1)
424
425    def htmlprinttabextvar(self,tabtodo):
426        self.fixedlt = True ## local time is real local time
427        self.gettitle()
428        print "<hr>"
429        print self.title
430        print "<hr>"
431        print "<ul>"
432        for i in range(len(tabtodo)): print "<li>" ; self.printextvar(tabtodo[i]) ; print "</li>"
433        print "</ul>"
434        print "<hr>"
435        print self.ack
436        print "<hr>"
437        #print "SETTINGS<br />"
438        #self.printcoord()
439        #self.printset()
440
441    def printmcd(self):
442    # 1. call MCD 2. print settings 3. print mean vars
443        self.update()
444        self.printcoord()
445        print "-------------------------------------------"
446        self.printmeanvar()
447
448########################
449### Several requests ###
450########################
451
452    def prepare(self,ndx=None,ndy=None):
453    ### prepare I/O arrays for 1d slices
454      if ndx is None:  print "No dimension in prepare. Exit. Set at least ndx." ; exit()
455      else:            self.xcoord = np.ones(ndx)
456      if ndy is None:  dashape = (ndx)     ; dashapemean = (ndx,6)     ; dashapeext = (ndx,101)     ; self.ycoord = None
457      else:            dashape = (ndx,ndy) ; dashapemean = (ndx,ndy,6) ; dashapeext = (ndx,ndy,101) ; self.ycoord = np.ones(ndy)
458      self.prestab = np.ones(dashape) ; self.denstab = np.ones(dashape) ; self.temptab = np.ones(dashape)
459      self.zonwindtab = np.ones(dashape) ; self.merwindtab = np.ones(dashape) 
460      self.meanvartab = np.ones(dashapemean) ; self.extvartab = np.ones(dashapeext)
461
462    def getextvar(self,num):
463    ### get a given var in extvartab
464      try: field=self.extvartab[:,:,num] 
465      except: field=self.extvartab[:,num]
466      return field
467
468    def definefield(self,choice):
469    ### for analysis or plot purposes, set field and field label from user-defined choice
470      choice = self.convertlab(choice)
471      field = self.getextvar(choice); fieldlab = self.getextvarlab(choice)
472      ## fix for possibly slightly negative tracers
473      if "(mol/mol)" in fieldlab or "(kg/kg)" in fieldlab or "(kg/m2)" in fieldlab or "(W/m2)" in fieldlab:
474         ind = np.where(field < 1.e-30)
475         if ind != -1: field[ind] = 0.e0 #1.e-30  ## 0 does not work everywhere.
476      return field,fieldlab
477
478    def ininterv(self,dstart,dend,nd,start=None,end=None,yaxis=False,vertcoord=False):
479    ### user-defined start and end are used to create xcoord (or ycoord) vector
480      if start is not None and end is not None:  first, second = self.correctbounds(start,end,vertcoord)
481      else:                                      first, second = self.correctbounds(dstart,dend,vertcoord) 
482      if self.zkey != 4 or not vertcoord:   tabtab = np.linspace(first,second,nd)
483      else:                                 tabtab = np.logspace(first,second,nd)
484      if not yaxis:      self.xcoord = tabtab
485      else:              self.ycoord = tabtab
486
487    def correctbounds(self,start,end,vertcoord):
488      if self.zkey != 4 or not vertcoord:
489        # regular altitudes
490        if start > end: first = end ; second = start
491        else:           first = start ; second = end
492      else:
493        # pressure: reversed avis
494        if start < end: first = np.log10(end) ; second = np.log10(start)
495        else:           first = np.log10(start) ; second = np.log10(end)
496      return first, second
497
498    def vertlabel(self):
499      if self.zkey == 1:   self.xlabel = "radius from centre of planet (m)"
500      elif self.zkey == 2: self.xlabel = "height above areoid (m) (MOLA zero datum)"
501      elif self.zkey == 3: self.xlabel = "height above surface (m)"
502      elif self.zkey == 4: self.xlabel = "pressure level (Pa)"
503      elif self.zkey == 5: self.xlabel = "altitude above mean Mars Radius(=3396000m) (m)"
504
505    def vertunits(self):
506      if self.zkey == 1:   self.vunits = "m CP"
507      elif self.zkey == 2: self.vunits = "m AMR"
508      elif self.zkey == 3: self.vunits = "m ALS"
509      elif self.zkey == 4: self.vunits = "Pa"
510      elif self.zkey == 5: self.vunits = "m AMMRad"
511
512    def vertaxis(self,number,yaxis=False):
513      if self.zkey == 2:   self.ininterv(-5000.,100000.,number,start=self.xzs,end=self.xze,yaxis=yaxis,vertcoord=True)
514      elif self.zkey == 3: self.ininterv(0.,120000.,number,start=self.xzs,end=self.xze,yaxis=yaxis,vertcoord=True)
515      elif self.zkey == 5: self.ininterv(-5000.,100000.,number,start=self.xzs,end=self.xze,yaxis=yaxis,vertcoord=True)
516      elif self.zkey == 4: self.ininterv(1000.,0.001,number,start=self.xzs,end=self.xze,yaxis=yaxis,vertcoord=True)
517      elif self.zkey == 1: self.ininterv(3396000,3596000,number,start=self.xzs,end=self.xze,yaxis=yaxis,vertcoord=True)
518
519###################
520### 1D analysis ###
521###################
522
523    def put1d(self,i):
524    ## fill in subscript i in output arrays
525    ## (arrays must have been correctly defined through prepare)
526      if self.prestab is None:  myplot.errormess("arrays must be prepared first through self.prepare")
527      self.prestab[i] = self.pres ; self.denstab[i] = self.dens ; self.temptab[i] = self.temp
528      self.zonwindtab[i] = self.zonwind ; self.merwindtab[i] = self.merwind
529      self.meanvartab[i,1:5] = self.meanvar[0:4]  ## note: var numbering according to MCD manual is kept
530      self.extvartab[i,1:100] = self.extvar[0:99] ## note: var numbering according to MCD manual is kept
531
532    def diurnal(self,nd=13):
533    ### retrieve a local time slice
534      self.fixedlt = True  ## local time is real local time
535      save = self.loct
536      self.xlabel = "Local time (Martian hour)"
537      self.prepare(ndx=nd) ; self.ininterv(0.,24.,nd,start=self.locts,end=self.locte) 
538      for i in range(nd): self.loct = self.xcoord[i] ; self.update() ; self.put1d(i)
539      self.loct = save
540
541    def zonal(self,nd=37):
542    ### retrieve a longitude slice
543      save = self.lon
544      self.xlabel = "East longitude (degrees)"
545      self.prepare(ndx=nd) ; self.ininterv(-180.,180.,nd,start=self.lons,end=self.lone)
546      if not self.fixedlt: umst = self.loct
547      for i in range(nd): 
548          self.lon = self.xcoord[i]
549          if not self.fixedlt: self.loct = (umst + self.lon/15.) % 24
550          self.update() ; self.put1d(i)
551      self.lon = save
552
553    def meridional(self,nd=19):
554    ### retrieve a latitude slice
555      self.fixedlt = True  ## local time is real local time
556      save = self.lat
557      self.xlabel = "North latitude (degrees)"
558      self.prepare(ndx=nd) ; self.ininterv(-90.,90.,nd,start=self.lats,end=self.late)
559      for i in range(nd): self.lat = self.xcoord[i] ; self.update() ; self.put1d(i)
560      self.lat = save
561
562    def profile(self,nd=20,tabperso=None):
563    ### retrieve an altitude slice (profile)
564      self.fixedlt = True  ## local time is real local time
565      save = self.xz
566      self.vertlabel()
567      self.vertplot = True
568      if tabperso is not None: nd = len(tabperso)
569      correct = False
570      self.prepare(ndx=nd) ; self.vertaxis(nd)
571      if tabperso is not None: self.xcoord = tabperso
572      for i in range(nd): self.xz = self.xcoord[i] ; self.update() ; self.put1d(i)
573      self.xz = save
574
575    def seasonal(self,nd=12):
576    ### retrieve a seasonal slice
577      save = self.xdate
578      self.xlabel = "Areocentric longitude (degrees)"
579      self.prepare(ndx=nd) ; self.ininterv(0.,360.,nd,start=self.xdates,end=self.xdatee)
580      for i in range(nd): self.xdate = self.xcoord[i] ; self.update() ; self.put1d(i)
581      self.xdate = save
582
583    def getascii(self,tabtodo,filename="output.txt"):
584    ### print out values in an ascii file
585      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
586      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
587      asciifile = open(filename, "w")
588      for i in range(len(tabtodo)): 
589          (field, fieldlab) = self.definefield(tabtodo[i])
590          self.gettitle(oneline=True)
591          asciifile.write("### " + self.title + "\n")
592          asciifile.write("### " + self.ack + "\n")
593          asciifile.write("### Column 1 is " + self.xlabel + "\n")
594          asciifile.write("### Column 2 is " + fieldlab + "\n")
595          for ix in range(len(self.xcoord)):
596              asciifile.write("%15.5e%15.5e\n" % ( self.xcoord[ix], field[ix] ) )
597      asciifile.close()
598      return 
599
600    def makeplot1d(self,choice):
601    ### one 1D plot is created for the user-defined variable in choice.
602      (field, fieldlab) = self.definefield(choice)
603      if not self.vertplot:  absc = self.xcoord ; ordo = field ; ordolab = fieldlab ; absclab = self.xlabel
604      else:                  ordo = self.xcoord ; absc = field ; absclab = fieldlab ; ordolab = self.xlabel
605      mpl.plot(absc,ordo,'-bo') ; mpl.ylabel(ordolab) ; mpl.xlabel(absclab) #; mpl.xticks(query.xcoord)
606      # cases with log axis
607      if self.zkey == 4: mpl.semilogy() ; ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
608      if not self.vertplot and self.islog: mpl.semilogy()
609      if self.vertplot and self.islog: mpl.semilogx()
610      mpl.figtext(0.5, 0.01, self.ack, ha='center')
611
612    def plot1d(self,tabtodo):
613    ### complete 1D figure with possible multiplots
614      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
615      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
616      fig = mpl.figure() ; subv,subh = myplot.definesubplot( len(tabtodo) , fig ) 
617      for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1).grid(True, linestyle=':', color='grey') ; self.makeplot1d(tabtodo[i])
618
619    def htmlplot1d(self,tabtodo,figname="temp.png",title=""):
620    ### complete 1D figure with possible multiplots
621    ### added in 09/2012 for online MCD
622    ### see http://www.dalkescientific.com/writings/diary/archive/2005/04/23/matplotlib_without_gui.html
623      from matplotlib.figure import Figure
624      from matplotlib.backends.backend_agg import FigureCanvasAgg
625      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
626      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
627
628      howmanyplots = len(tabtodo)
629      if howmanyplots == 1: fig = Figure(figsize=(16,8))
630      elif howmanyplots == 2: fig = Figure(figsize=(8,8))
631      elif howmanyplots == 3: fig = Figure(figsize=(8,16))
632      elif howmanyplots == 4: fig = Figure(figsize=(16,8))
633
634      subv,subh = myplot.definesubplot( len(tabtodo) , fig )
635      for i in range(len(tabtodo)):
636        yeah = fig.add_subplot(subv,subh,i+1) #.grid(True, linestyle=':', color='grey')
637        choice = tabtodo[i]
638        (field, fieldlab) = self.definefield(choice)
639
640        # in log plots we do not show the negative values
641        if self.islog: field[np.where(field <= 0.e0)] = np.nan
642
643        if not self.vertplot:  absc = self.xcoord ; ordo = field ; ordolab = fieldlab ; absclab = self.xlabel
644        else:                  ordo = self.xcoord ; absc = field ; absclab = fieldlab ; ordolab = self.xlabel
645
646        yeah.plot(absc,ordo,'-bo') #; mpl.xticks(query.xcoord)
647        ax = fig.gca() ; ax.set_ylabel(ordolab) ; ax.set_xlabel(absclab)
648
649        if self.xzs is not None and self.zkey == 4: ax.set_yscale('log') ; ax.set_ylim(ax.get_ylim()[::-1])
650        if not self.vertplot and self.islog: ax.set_yscale('log')
651        if self.vertplot and self.islog: ax.set_xscale('log')
652
653        if self.lats is not None:      ax.set_xticks(np.arange(-90,91,15)) ; ax.set_xbound(lower=self.lats, upper=self.late)
654        elif self.lons is not None:    ax.set_xticks(np.arange(-360,361,30)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
655        elif self.locts is not None:   ax.set_xticks(np.arange(0,26,2)) ; ax.set_xbound(lower=self.locts, upper=self.locte)
656
657        ax.grid(True, linestyle=':', color='grey')
658
659      self.gettitle()
660      fig.text(0.5, 0.95, self.title, ha='center')
661      fig.text(0.5, 0.01, self.ack, ha='center')
662      canvas = FigureCanvasAgg(fig)
663      # The size * the dpi gives the final image size
664      #   a4"x4" image * 80 dpi ==> 320x320 pixel image
665      canvas.print_figure(figname, dpi=self.dpi)
666
667###################
668### 2D analysis ###
669###################
670
671    def latlon(self,ndx=37,ndy=19):
672    ### retrieve a latitude/longitude slice
673    ### default is: local time is not fixed. user-defined local time is at longitude 0.
674      save1 = self.lon ; save2 = self.lat ; save3 = self.loct
675      self.xlabel = "East longitude (degrees)" ; self.ylabel = "North latitude (degrees)"
676      self.prepare(ndx=ndx,ndy=ndy)
677      self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
678      self.ininterv(-90.,  90.,ndy,start=self.lats,end=self.late,yaxis=True)
679      if not self.fixedlt: umst = self.loct
680      for i in range(ndx):
681       for j in range(ndy):
682         self.lon = self.xcoord[i] ; self.lat = self.ycoord[j]
683         if not self.fixedlt: self.loct = (umst + self.lon/15.) % 24
684         self.update() ; self.put2d(i,j)
685      if not self.fixedlt: self.loct = umst
686      self.lon = save1 ; self.lat = save2 ; self.loct = save3
687
688    def secalt(self,ndx=37,ndy=20,typex="lat"):
689    ### retrieve a coordinate/altitude slice
690      save1 = self.lon ; save2 = self.xz ; save3 = self.loct ; save4 = self.lat
691      self.prepare(ndx=ndx,ndy=ndy)
692      self.vertlabel() ; self.ylabel = self.xlabel
693      self.vertaxis(ndy,yaxis=True)
694      if "lat" in typex:
695          self.xlabel = "North latitude (degrees)"
696          self.ininterv(-90.,90.,ndx,start=self.lats,end=self.late)
697      elif typex == "lon":
698          self.xlabel = "East longitude (degrees)"
699          self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
700      if not self.fixedlt: umst = self.loct
701      for i in range(ndx):
702       for j in range(ndy):
703         if typex == "lat":   self.lat = self.xcoord[i]
704         elif typex == "lon": self.lon = self.xcoord[i]
705         self.xz = self.ycoord[j]
706         if not self.fixedlt: self.loct = (umst + self.lon/15.) % 24
707         self.update() ; self.put2d(i,j)
708      if not self.fixedlt: self.loct = umst
709      self.lon = save1 ; self.xz = save2 ; self.loct = save3 ; self.lat = save4
710
711    def zonalmean(self,ndx=37,ndy=20,ndmean=32):
712    ### retrieve a zonalmean lat/altitude slice
713      self.fixedlt = False
714      save1 = self.lon ; save2 = self.xz ; save3 = self.loct ; save4 = self.lat
715      self.prepare(ndx=ndx,ndy=ndy)
716      self.vertlabel() ; self.ylabel = self.xlabel
717      self.vertaxis(ndy,yaxis=True)
718      self.xlabel = "North latitude (degrees)"
719      self.ininterv(-180.,180.,ndmean)
720      coordmean = self.xcoord
721      self.ininterv(-90.,90.,ndx,start=self.lats,end=self.late)
722      umst = self.loct #fixedlt false for this case
723      for i in range(ndx):
724       self.lat = self.xcoord[i]
725       for j in range(ndy):
726        self.xz = self.ycoord[j]
727        meanpres = 0. ; meandens = 0. ; meantemp = 0. ; meanzonwind = 0. ; meanmerwind = 0. ; meanmeanvar = np.zeros(5) ; meanextvar = np.zeros(100)       
728        for m in range(ndmean):
729           self.lon = coordmean[m]
730           self.loct = (umst + self.lon/15.) % 24 #fixedlt false for this case
731           self.update() 
732           meanpres = meanpres + self.pres/float(ndmean) ; meandens = meandens + self.dens/float(ndmean) ; meantemp = meantemp + self.temp/float(ndmean)
733           meanzonwind = meanzonwind + self.zonwind/float(ndmean) ; meanmerwind = meanmerwind + self.merwind/float(ndmean)
734           meanmeanvar = meanmeanvar + self.meanvar/float(ndmean) ; meanextvar = meanextvar + self.extvar/float(ndmean)
735        self.pres=meanpres ; self.dens=meandens ; self.temp=meantemp ; self.zonwind=meanzonwind ; self.merwind=meanmerwind
736        self.meanvar=meanmeanvar ; self.extvar=meanextvar
737        self.put2d(i,j)
738      self.loct = umst #fixedlt false for this case
739      self.lon = save1 ; self.xz = save2 ; self.loct = save3 ; self.lat = save4
740
741    def hovmoller(self,ndtime=25,ndcoord=20,typex="lat"):
742    ### retrieve a time/other coordinate slice
743      save1 = self.lat ; save2 = self.xz ; save3 = self.loct ; save4 = self.lon
744      if typex == "lat": 
745          ndx = ndcoord ; self.xlabel = "North latitude (degrees)" 
746          ndy = ndtime ; self.ylabel = "Local time (Martian hour)"
747          self.prepare(ndx=ndx,ndy=ndy)
748          self.ininterv(-90.,90.,ndx,start=self.lats,end=self.late)
749          self.ininterv(0.,24.,ndy,start=self.locts,end=self.locte,yaxis=True)
750      elif typex == "lon":
751          ndx = ndcoord ; self.xlabel = "East longitude (degrees)"
752          ndy = ndtime ; self.ylabel = "Local time (Martian hour)"
753          self.prepare(ndx=ndx,ndy=ndy)
754          self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
755          self.ininterv(0.,24.,ndy,start=self.locts,end=self.locte,yaxis=True)
756      elif typex == "alt":
757          ndy = ndcoord ; self.vertlabel() ; self.ylabel = self.xlabel
758          ndx = ndtime ; self.xlabel = "Local time (Martian hour)"
759          self.prepare(ndx=ndx,ndy=ndy)
760          self.vertaxis(ndy,yaxis=True)
761          self.ininterv(0.,24.,ndx,start=self.locts,end=self.locte)
762      for i in range(ndx):
763       for j in range(ndy):
764         if typex == "lat":   self.lat = self.xcoord[i] ; self.loct = self.ycoord[j]
765         elif typex == "lon": self.lon = self.xcoord[i] ; self.loct = self.ycoord[j]
766         elif typex == "alt": self.xz = self.ycoord[j] ; self.loct = self.xcoord[i]
767         self.update() ; self.put2d(i,j)
768      self.lat = save1 ; self.xz = save2 ; self.loct = save3 ; self.lon = save4
769
770    def put2d(self,i,j):
771    ## fill in subscript i,j in output arrays
772    ## (arrays must have been correctly defined through prepare)
773      if self.prestab is None:  myplot.errormess("arrays must be prepared first through self.prepare")
774      self.prestab[i,j] = self.pres ; self.denstab[i,j] = self.dens ; self.temptab[i,j] = self.temp
775      self.zonwindtab[i,j] = self.zonwind ; self.merwindtab[i,j] = self.merwind
776      self.meanvartab[i,j,1:5] = self.meanvar[0:4]  ## note: var numbering according to MCD manual is kept
777      self.extvartab[i,j,1:100] = self.extvar[0:99] ## note: var numbering according to MCD manual is kept
778
779    def makemap2d(self,choice,incwind=False,proj="cyl"):
780    ### one 2D map is created for the user-defined variable in choice.
781      self.latlon() ## a map is implicitely a lat-lon plot. otherwise it is a plot (cf. makeplot2d)
782      if choice == "wind" or incwind:
783          (windx, fieldlabwx) = self.definefield("u")
784          (windy, fieldlabwy) = self.definefield("v")
785      if choice == "wind":
786          field = np.sqrt(windx*windx + windy*windy)
787          fieldlab = "Horizontal wind speed (m/s)"
788      else:   
789          (field, fieldlab) = self.definefield(choice)
790      if incwind:   myplot.maplatlon(self.xcoord,self.ycoord,field,title=fieldlab,proj=proj,vecx=windx,vecy=windy) #,stride=1)
791      else:         myplot.maplatlon(self.xcoord,self.ycoord,field,title=fieldlab,proj=proj)
792      mpl.figtext(0.5, 0.0, self.ack, ha='center')
793
794    def map2d(self,tabtodo,incwind=False,proj="cyl"):
795    ### complete 2D figure with possible multiplots
796      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
797      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
798      fig = mpl.figure()
799      subv,subh = myplot.definesubplot( len(tabtodo) , fig ) 
800      for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1) ; self.makemap2d(tabtodo[i],incwind=incwind,proj=proj)
801
802    def htmlmap2d(self,tabtodo,incwind=False,figname="temp.png",back="zMOL"):
803    ### complete 2D figure with possible multiplots
804    ### added in 09/2012 for online MCD
805    ### see http://www.dalkescientific.com/writings/diary/archive/2005/04/23/matplotlib_without_gui.html
806      from matplotlib.figure import Figure
807      from matplotlib.backends.backend_agg import FigureCanvasAgg
808      from matplotlib.cm import get_cmap
809      from matplotlib import rcParams
810      #from mpl_toolkits.basemap import Basemap # does not work
811      from Scientific.IO import NetCDF
812
813      filename = "/home/marshttp/surface.nc"
814      zefile = NetCDF.NetCDFFile(filename, 'r') 
815      fieldc = zefile.variables[back]
816      yc = zefile.variables['latitude']
817      xc = zefile.variables['longitude']
818
819      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
820      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
821
822      howmanyplots = len(tabtodo)
823      if howmanyplots == 1: fig = Figure(figsize=(16,8)) 
824      elif howmanyplots == 2: fig = Figure(figsize=(8,8)) 
825      elif howmanyplots == 3: fig = Figure(figsize=(8,16)) 
826      elif howmanyplots == 4: fig = Figure(figsize=(16,8)) 
827
828      subv,subh = myplot.definesubplot( len(tabtodo) , fig )
829
830      for i in range(len(tabtodo)):
831        yeah = fig.add_subplot(subv,subh,i+1)
832        choice = tabtodo[i]
833        self.latlon(ndx=64,ndy=48) 
834        ## a map is implicitely a lat-lon plot. otherwise it is a plot (cf. makeplot2d)
835        (field, fieldlab) = self.definefield(choice)
836        if incwind: (windx, fieldlabwx) = self.definefield("u") ; (windy, fieldlabwy) = self.definefield("v")
837
838        proj="moll" ; colorb= self.colorm ; ndiv=20 ; zeback="molabw" ; trans=1.0 #0.6
839        vecx=None ; vecy=None ; stride=2
840        lon = self.xcoord
841        lat = self.ycoord
842       
843        #[lon2d,lat2d] = np.meshgrid(lon,lat)
844        ##### define projection and background. define x and y given the projection
845        ##[wlon,wlat] = myplot.latinterv()
846        ##yeahm = myplot.define_proj(proj,wlon,wlat,back=zeback,blat=None,blon=None)
847        ##x, y = yeahm(lon2d, lat2d)
848        #map = Basemap(projection='ortho',lat_0=45,lon_0=-100)
849        #x, y = map(lon2d, lat2d)
850
851        #### TEMP
852        x = lon ; y = lat
853
854        ## define field. bound field.
855        what_I_plot = np.transpose(field)
856        zevmin, zevmax = myplot.calculate_bounds(what_I_plot,vmin=self.min2d,vmax=self.max2d) 
857        what_I_plot = myplot.bounds(what_I_plot,zevmin,zevmax)
858        ## define contour field levels. define color palette
859        ticks = ndiv + 1
860        zelevels = np.linspace(zevmin,zevmax,ticks)
861        palette = get_cmap(name=colorb)
862
863        # You can set negative contours to be solid instead of dashed:
864        rcParams['contour.negative_linestyle'] = 'solid'
865        ## contours topo
866        zelevc = np.linspace(-9.,20.,11)
867        yeah.contour( xc, yc, fieldc, zelevc, colors='black',linewidths = 0.4)
868        yeah.contour( np.array(xc) + 360., yc, fieldc, zelevc, colors='black',linewidths = 0.4)
869        yeah.contour( np.array(xc) - 360., yc, fieldc, zelevc, colors='black',linewidths = 0.4)
870        # contour field
871        c = yeah.contourf( x, y, what_I_plot, zelevels, cmap = palette, alpha = trans )
872        clb = Figure.colorbar(fig,c,orientation='vertical',format=self.fmt,ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])))
873        clb.set_label(fieldlab)
874        if incwind:
875          [x2d,y2d] = np.meshgrid(x,y)
876          yeah.quiver(x2d,y2d,np.transpose(windx),np.transpose(windy))
877        ax = fig.gca() ; ax.set_ylabel("Latitude") ; ax.set_xlabel("Longitude")
878        ax.set_xticks(np.arange(-360,361,45)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
879        ax.set_yticks(np.arange(-90,91,30)) ; ax.set_ybound(lower=self.lats, upper=self.late)
880      self.gettitle()
881      fig.text(0.5, 0.95, self.title, ha='center')
882      fig.text(0.5, 0.01, self.ack, ha='center')
883      canvas = FigureCanvasAgg(fig)
884      # The size * the dpi gives the final image size
885      #   a4"x4" image * 80 dpi ==> 320x320 pixel image
886      canvas.print_figure(figname, dpi=self.dpi)
887
888    def htmlplot2d(self,tabtodo,figname="temp.png"):
889    ### complete 2D figure with possible multiplots
890    ### added in 10/2012 for online MCD
891    ### see http://www.dalkescientific.com/writings/diary/archive/2005/04/23/matplotlib_without_gui.html
892      from matplotlib.figure import Figure
893      from matplotlib.backends.backend_agg import FigureCanvasAgg
894      from matplotlib.cm import get_cmap
895      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
896      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
897
898      howmanyplots = len(tabtodo)
899      if howmanyplots == 1: fig = Figure(figsize=(16,8))
900      elif howmanyplots == 2: fig = Figure(figsize=(8,8))
901      elif howmanyplots == 3: fig = Figure(figsize=(8,16))
902      elif howmanyplots == 4: fig = Figure(figsize=(16,8))
903
904      subv,subh = myplot.definesubplot( len(tabtodo) , fig )
905
906      for i in range(len(tabtodo)):
907        yeah = fig.add_subplot(subv,subh,i+1)
908        choice = tabtodo[i]
909
910        if self.lons is not None:   
911           if self.locts is None:  self.secalt(ndx=64,ndy=35,typex="lon")
912           else:                   self.hovmoller(ndcoord=64,typex="lon")
913        elif self.lats is not None: 
914           if self.locts is None: 
915               if self.zonmean:   self.zonalmean()
916               else:         self.secalt(ndx=48,ndy=35,typex="lat")
917           else:                   self.hovmoller(ndcoord=48,typex="lat")
918        else:
919           self.hovmoller(ndcoord=35,typex="alt")
920
921        (field, fieldlab) = self.definefield(choice)
922
923        colorb=self.colorm ; ndiv=20 
924
925        ## define field. bound field.
926        what_I_plot = np.transpose(field)
927        zevmin, zevmax = myplot.calculate_bounds(what_I_plot,vmin=self.min2d,vmax=self.max2d) 
928        what_I_plot = myplot.bounds(what_I_plot,zevmin,zevmax)
929        ## define contour field levels. define color palette
930        ticks = ndiv + 1
931        zelevels = np.linspace(zevmin,zevmax,ticks)
932        palette = get_cmap(name=colorb)
933        # contour field
934        c = yeah.contourf( self.xcoord, self.ycoord, what_I_plot, zelevels, cmap = palette )
935        clb = Figure.colorbar(fig,c,orientation='vertical',format=self.fmt,ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])))
936        clb.set_label(fieldlab)
937        ax = fig.gca() ; ax.set_ylabel(self.ylabel) ; ax.set_xlabel(self.xlabel)
938
939        if self.lons is not None:   ax.set_xticks(np.arange(-360,361,45)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
940        elif self.lats is not None: ax.set_xticks(np.arange(-90,91,30)) ; ax.set_xbound(lower=self.lats, upper=self.late)
941
942        if self.locts is not None: 
943            if self.xzs is not None: ax.set_xticks(np.arange(0,26,2)) ; ax.set_xbound(lower=self.locts, upper=self.locte)
944            else:                    ax.set_yticks(np.arange(0,26,2)) ; ax.set_ybound(lower=self.locts, upper=self.locte)
945
946        if self.zkey == 4 and self.xzs is not None: 
947            ax.set_yscale('log') ; ax.set_ylim(ax.get_ylim()[::-1])
948        else:
949            #ax.set_yticks(np.arange(self.xzs,self.xze,10000.)) ;
950            ax.set_ybound(lower=self.xzs, upper=self.xze)
951
952      self.gettitle()
953      fig.text(0.5, 0.95, self.title, ha='center')
954      fig.text(0.5, 0.01, self.ack, ha='center')
955      canvas = FigureCanvasAgg(fig)
956      # The size * the dpi gives the final image size
957      #   a4"x4" image * 80 dpi ==> 320x320 pixel image
958      canvas.print_figure(figname, dpi=self.dpi)
959
960    ### TODO: makeplot2d, plot2d, passer plot settings
961
Note: See TracBrowser for help on using the repository browser.