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

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

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