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

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

UTIL PYTHON : a more definitive version of what could make a new and simple MCD web interface. added a lot of capabilities: variable setting, beginner mode, wind vector, fixed or not local time, etc etc etc. some modifications were made to the python MCD interface too (this makes the basis for the web interface).

  • Property svn:executable set to *
File size: 16.3 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 fmcd
11import matplotlib.pyplot as mpl
12import myplot
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 output"
29        self.dset      = '/home/aymeric/Science/MCD_v4.3/data/'
30        ## 1. spatio-temporal coordinates
31        self.lat       = 0.
32        self.lon       = 0.
33        self.loct      = 0.
34        self.xdate     = 0.  # see datekey
35        self.xz        = 10. # see zkey
36        ## 1bis. related settings
37        self.zkey      = 3  # specify that xz is the altitude above surface (m)
38        self.datekey   = 1  # 0 = "Earth time": xdate is given in Julian days (localtime must be set to zero)
39                            # 1 = "Mars date": xdate is the value of Ls
40        ## 2. climatological options
41        self.dust      = 2  #our best guess MY24 scenario, with solar average conditions
42        self.hrkey     = 1  #set high resolution mode on (hrkey=0 to set high resolution off)
43        ## 3. additional settings for advanced use
44        self.extvarkey = 1  #extra output variables (1: yes, 0: no)
45        self.perturkey = 0  #integer perturkey ! perturbation type (0: none)
46        self.seedin    = 1  #random number generator seed (unused if perturkey=0)
47        self.gwlength  = 0. #gravity Wave wavelength (unused if perturkey=0)
48        ## outputs. just to define attributes.
49        ## --> in update
50        self.pres = None ; self.dens = None ; self.temp = None ; self.zonwind = None ; self.merwind = None ; self.meanvar = None ; self.extvar = None
51        self.seedout = None ; self.ierr = None
52        ## --> in prepare
53        self.xcoord = None ; self.ycoord = None
54        self.prestab = None ; self.denstab = None ; self.temptab = None 
55        self.zonwindtab = None ; self.merwindtab = None ; self.meanvartab = None ; self.extvartab = None
56
57    def viking1(self): self.name = "Viking 1 site. MCD v4.3 output" ; self.lat = 22.48 ; self.lon = -49.97 ; self.xdate = 97.
58    def viking2(self): self.name = "Viking 2 site. MCD v4.3 output" ; self.lat = 47.97 ; self.lon = -225.74 ; self.xdate = 117.6
59
60    def getextvarlab(self,num):
61        whichfield = { \
62        91: "Pressure (Pa)", \
63        92: "Density (kg/m3)", \
64        93: "Temperature (K)", \
65        94: "W-E wind component (m/s)", \
66        95: "S-N wind component (m/s)", \
67        1: "Radial distance from planet center (m)",\
68        2: "Altitude above areoid (Mars geoid) (m)",\
69        3: "Altitude above local surface (m)",\
70        4: "orographic height (m) (surface altitude above areoid)",\
71        5: "Ls, solar longitude of Mars (deg)",\
72        6: "LST local true solar time (hrs)",\
73        7: "Universal solar time (LST at lon=0) (hrs)",\
74        8: "Air heat capacity Cp (J kg-1 K-1)",\
75        9: "gamma=Cp/Cv Ratio of specific heats",\
76        10: "density RMS day to day variations (kg/m^3)",\
77        11: "[not defined]",\
78        12: "[not defined]",\
79        13: "scale height H(p) (m)",\
80        14: "GCM orography (m)",\
81        15: "surface temperature (K)",\
82        16: "daily maximum mean surface temperature (K)",\
83        17: "daily minimum mean surface temperature (K)",\
84        18: "surf. temperature RMS day to day variations (K)",\
85        19: "surface pressure (high resolution if hireskey=1)",\
86        20: "GCM surface pressure (Pa)",\
87        21: "atmospheric pressure RMS day to day variations (Pa)",\
88        22: "surface pressure RMS day to day variations (Pa)",\
89        23: "temperature RMS day to day variations (K)",\
90        24: "zonal wind RMS day to day variations (m/s)",\
91        25: "meridional wind RMS day to day variations (m/s)",\
92        26: "vertical wind component (m/s) >0 when downwards!",\
93        27: "vertical wind RMS day to day variations (m/s)",\
94        28: "small scale perturbation (gravity wave) (kg/m^3)",\
95        29: "q2: turbulent kinetic energy (m2/s2)",\
96        30: "[not defined]",\
97        31: "thermal IR flux to surface (W/m2)",\
98        32: "solar flux to surface (W/m2)",\
99        33: "thermal IR flux to space (W/m2)",\
100        34: "solar flux reflected to space (W/m2)",\
101        35: "surface CO2 ice layer (kg/m2)",\
102        36: "DOD: Dust column visible optical depth",\
103        37: "Dust mass mixing ratio (kg/kg)",\
104        38: "DOD RMS day to day variations",\
105        39: "DOD total standard deviation over season",\
106        40: "Water vapor column (kg/m2)",\
107        41: "Water vapor vol. mixing ratio (mol/mol)",\
108        42: "Water ice column (kg/m2)",\
109        43: "Water ice mixing ratio (mol/mol)",\
110        44: "O3 ozone vol. mixing ratio (mol/mol)",\
111        45: "[CO2] vol. mixing ratio (mol/mol)",\
112        46: "[O] vol. mixing ratio (mol/mol)",\
113        47: "[N2] vol. mixing ratio (mol/mol)",\
114        48: "[CO] vol. mixing ratio (mol/mol)",\
115        49: "R: Molecular gas constant (J K-1 kg-1)",\
116        50: "Air viscosity estimation (N s m-2)"
117        }
118        if num not in whichfield: myplot.errormess("Incorrect subscript in extvar.")
119        return whichfield[num]
120
121    def convertlab(self,num):       
122        ## a conversion from text inquiries to extvar numbers. to be completed.
123        if num == "p": num = 91
124        elif num == "rho": num = 92
125        elif num == "t": num = 93
126        elif num == "u": num = 94
127        elif num == "v": num = 95
128        elif num == "tsurf": num = 15
129        elif num == "topo": num = 4
130        elif num == "h": num = 13
131        elif num == "ps": num = 19
132        elif num == "tau": num = 36
133        elif num == "mtot": num = 40
134        elif num == "icetot": num = 42
135        elif num == "ps_ddv": num = 22
136        elif num == "h2ovap": num = 41
137        elif num == "h2oice": num = 43
138        elif num == "cp": num = 8
139        elif num == "rho_ddv": num = 10
140        elif num == "tsurfmx": num = 16
141        elif num == "tsurfmn": num = 17
142        elif num == "lwdown": num = 31
143        elif num == "swdown": num = 32
144        elif num == "lwup": num = 33
145        elif num == "swup": num = 34
146        elif num == "o3": num = 44
147        elif num == "o": num = 46
148        elif num == "co": num = 48
149        elif num == "visc": num = 50
150        elif num == "co2ice": num = 35
151        elif not isinstance(num, np.int): myplot.errormess("field reference not found.")
152        return num
153
154###################
155### One request ###
156###################
157
158    def update(self):
159    # retrieve fields from MCD (call_mcd). more info in fmcd.call_mcd.__doc__
160        (self.pres, self.dens, self.temp, self.zonwind, self.merwind, \
161         self.meanvar, self.extvar, self.seedout, self.ierr) \
162         = \
163         fmcd.call_mcd(self.zkey,self.xz,self.lon,self.lat,self.hrkey, \
164             self.datekey,self.xdate,self.loct,self.dset,self.dust, \
165             self.perturkey,self.seedin,self.gwlength,self.extvarkey )
166        ## we use the end of extvar (unused) to store meanvar. this is convenient for getextvar(lab)
167        self.extvar[90] = self.pres ; self.extvar[91] = self.dens
168        self.extvar[92] = self.temp ; self.extvar[93] = self.zonwind ; self.extvar[94] = self.merwind
169
170    def printset(self):
171    # print main settings
172        print "zkey",self.zkey,"xz",self.xz,"lon",self.lon,"lat",self.lat,"hrkey",self.hrkey, \
173              "xdate",self.xdate,"loct",self.loct,"dust",self.dust
174
175    def getnameset(self):
176    # set a name referring to settings [convenient for databases]
177        name = str(self.zkey)+str(self.xz)+str(self.lon)+str(self.lat)+str(self.hrkey)+str(self.datekey)+str(self.xdate)+str(self.loct)+str(self.dust)
178        return name
179
180    def printcoord(self):
181    # print requested space-time coordinates
182        print "LAT",self.lat,"LON",self.lon,"LOCT",self.loct,"XDATE",self.xdate
183
184    def printmeanvar(self):
185    # print mean MCD variables
186        print "Pressure = %5.3f pascals. " % (self.pres)
187        print "Density = %5.3f kilograms per cubic meter. " % (self.dens)
188        print "Temperature = %3.0f kelvins (%4.0f degrees celsius)." % (self.temp,self.temp-273.15)
189        print "Zonal wind = %5.3f meters per second." % (self.zonwind)
190        print "Meridional wind = %5.3f meters per second." % (self.merwind)
191
192    def printextvar(self,num):
193    # print extra MCD variables
194        num = self.convertlab(num)
195        print self.getextvarlab(num) + " ..... " + str(self.extvar[num-1])
196
197    def printallextvar(self):
198    # print all extra MCD variables   
199        for i in range(50): self.printextvar(i+1)
200
201    def htmlprinttabextvar(self,tabtodo):
202        print "Results from the Mars Climate Database"
203        print "<ul>"
204        for i in range(len(tabtodo)): print "<li>" ; self.printextvar(tabtodo[i]) ; print "</li>"
205        print "</ul>"
206        print "<hr>"
207        print "SETTINGS<br />"
208        self.printcoord()
209        self.printset()
210
211    def printmcd(self):
212    # 1. call MCD 2. print settings 3. print mean vars
213        self.update()
214        self.printcoord()
215        print "-------------------------------------------"
216        self.printmeanvar()
217
218########################
219### Several requests ###
220########################
221
222    def prepare(self,ndx=None,ndy=None):
223    ### prepare I/O arrays for 1d slices
224      if ndx is None:  print "No dimension in prepare. Exit. Set at least ndx." ; exit()
225      else:            self.xcoord = np.ones(ndx)
226      if ndy is None:  dashape = (ndx)     ; dashapemean = (ndx,6)     ; dashapeext = (ndx,101)     ; self.ycoord = None
227      else:            dashape = (ndx,ndy) ; dashapemean = (ndx,ndy,6) ; dashapeext = (ndx,ndy,101) ; self.ycoord = np.ones(ndy)
228      self.prestab = np.ones(dashape) ; self.denstab = np.ones(dashape) ; self.temptab = np.ones(dashape)
229      self.zonwindtab = np.ones(dashape) ; self.merwindtab = np.ones(dashape) 
230      self.meanvartab = np.ones(dashapemean) ; self.extvartab = np.ones(dashapeext)
231
232    def getextvar(self,num):
233    ### get a given var in extvartab
234      try: field=self.extvartab[:,:,num] 
235      except: field=self.extvartab[:,num]
236      return field
237
238    def definefield(self,choice):
239    ### for analysis or plot purposes, set field and field label from user-defined choice
240      choice = self.convertlab(choice)
241      field = self.getextvar(choice); fieldlab = self.getextvarlab(choice)
242      return field,fieldlab
243
244###################
245### 1D analysis ###
246###################
247
248    def put1d(self,i):
249    ## fill in subscript i in output arrays
250    ## (arrays must have been correctly defined through prepare)
251      if self.prestab is None:  myplot.errormess("arrays must be prepared first through self.prepare")
252      self.prestab[i] = self.pres ; self.denstab[i] = self.dens ; self.temptab[i] = self.temp
253      self.zonwindtab[i] = self.zonwind ; self.merwindtab[i] = self.merwind
254      self.meanvartab[i,1:5] = self.meanvar[0:4]  ## note: var numbering according to MCD manual is kept
255      self.extvartab[i,1:100] = self.extvar[0:99] ## note: var numbering according to MCD manual is kept
256
257    def diurnal(self,nd=13,start=0.,end=24.):
258    ### retrieve a local time slice
259      self.xlabel = "Local time (Martian hour)"
260      self.prepare(ndx=nd) ; self.xcoord = np.linspace(start,end,nd)
261      for i in range(nd): self.loct = self.xcoord[i] ; self.update() ; self.put1d(i)
262
263    def zonal(self,nd=37,start=-180.,end=180.):
264    ### retrieve a longitude slice
265      self.xlabel = "East longitude (degrees)"
266      self.prepare(ndx=nd) ; self.xcoord = np.linspace(start,end,nd)
267      for i in range(nd): self.lon = self.xcoord[i] ; self.update() ; self.put1d(i)
268
269    def meridional(self,nd=19,start=-90.,end=90.):
270    ### retrieve a latitude slice
271      self.xlabel = "North latitude (degrees)"
272      self.prepare(ndx=nd) ; self.xcoord = np.linspace(start,end,nd)
273      for i in range(nd): self.lat = self.xcoord[i] ; self.update() ; self.put1d(i)
274
275    def profile(self,nd=20,start=0.,end=120000.,tabperso=None):
276    ### retrieve an altitude slice (profile)
277      self.xlabel = "Altitude (m)"
278      if tabperso is not None: nd = len(tabperso)
279      self.prepare(ndx=nd) 
280      if tabperso is None:  self.xcoord = np.linspace(start,end,nd)
281      else:                 self.xcoord = tabperso
282      for i in range(nd): self.xz = self.xcoord[i] ; self.update() ; self.put1d(i)
283
284    def seasonal(self,nd=12,start=0.,end=360.):
285    ### retrieve a seasonal slice
286      self.xlabel = "Areocentric longitude (degrees)"
287      self.prepare(ndx=nd) ; self.xcoord = np.linspace(start,end,nd)
288      for i in range(nd): self.xdate = self.xcoord[i] ; self.update() ; self.put1d(i)
289
290    def makeplot1d(self,choice,vertplot=0):
291    ### one 1D plot is created for the user-defined variable in choice.
292      (field, fieldlab) = self.definefield(choice)
293      if vertplot != 1:  absc = self.xcoord ; ordo = field ; ordolab = fieldlab ; absclab = self.xlabel
294      else:              ordo = self.xcoord ; absc = field ; absclab = fieldlab ; ordolab = self.xlabel
295      mpl.plot(absc,ordo,'-bo') ; mpl.ylabel(ordolab) ; mpl.xlabel(absclab) #; mpl.xticks(query.xcoord)
296      mpl.figtext(0.5, 0.01, "Mars Climate Database (c) LMD/OU/IAA/ESA/CNES", ha='center')
297
298    def plot1d(self,tabtodo,vertplot=0):
299    ### complete 1D figure with possible multiplots
300      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
301      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
302      fig = mpl.figure() ; subv,subh = myplot.definesubplot( len(tabtodo) , fig ) 
303      for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1).grid(True, linestyle=':', color='grey') ; self.makeplot1d(tabtodo[i],vertplot)
304
305###################
306### 2D analysis ###
307###################
308
309    def latlon(self,ndx=37,startx=-180.,endx=180.,ndy=19,starty=-90.,endy=90.,fixedlt=False):
310    ### retrieve a latitude/longitude slice
311    ### default is: local time is not fixed. user-defined local time is at longitude 0.
312      self.xlabel = "East longitude (degrees)" ; self.ylabel = "North latitude (degrees)"
313      self.prepare(ndx=ndx,ndy=ndy)
314      self.xcoord = np.linspace(startx,endx,ndx) ; self.ycoord = np.linspace(starty,endy,ndy)
315      if not fixedlt: umst = self.loct
316      for i in range(ndx):
317       for j in range(ndy):
318         self.lon = self.xcoord[i] ; self.lat = self.ycoord[j]
319         if not fixedlt: self.loct = (umst + self.lon/15.) % 24
320         self.update() ; self.put2d(i,j)
321      if not fixedlt: self.loct = umst
322
323    def put2d(self,i,j):
324    ## fill in subscript i,j in output arrays
325    ## (arrays must have been correctly defined through prepare)
326      if self.prestab is None:  myplot.errormess("arrays must be prepared first through self.prepare")
327      self.prestab[i,j] = self.pres ; self.denstab[i,j] = self.dens ; self.temptab[i,j] = self.temp
328      self.zonwindtab[i,j] = self.zonwind ; self.merwindtab[i,j] = self.merwind
329      self.meanvartab[i,j,1:5] = self.meanvar[0:4]  ## note: var numbering according to MCD manual is kept
330      self.extvartab[i,j,1:100] = self.extvar[0:99] ## note: var numbering according to MCD manual is kept
331
332    def makemap2d(self,choice,incwind=False,fixedlt=False):
333    ### one 2D map is created for the user-defined variable in choice.
334      self.latlon(fixedlt=fixedlt) ## a map is implicitely a lat-lon plot. otherwise it is a plot (cf. makeplot2d)
335      (field, fieldlab) = self.definefield(choice)
336      if incwind:
337          (windx, fieldlabwx) = self.definefield("u")
338          (windy, fieldlabwy) = self.definefield("v")
339          myplot.maplatlon(self.xcoord,self.ycoord,field,title=fieldlab,proj="cyl",vecx=windx,vecy=windy)
340      else:
341          myplot.maplatlon(self.xcoord,self.ycoord,field,title=fieldlab,proj="moll")
342      mpl.figtext(0.5, 0.0, "Mars Climate Database (c) LMD/OU/IAA/ESA/CNES", ha='center')
343
344    def map2d(self,tabtodo,incwind=False,fixedlt=False):
345    ### complete 2D figure with possible multiplots
346      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
347      if isinstance(tabtodo,np.int): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
348      fig = mpl.figure() ; subv,subh = myplot.definesubplot( len(tabtodo) , fig ) 
349      for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1) ; self.makemap2d(tabtodo[i],incwind=incwind,fixedlt=fixedlt)
350
351    ### TODO: makeplot2d, plot2d, passer plot settings
352
Note: See TracBrowser for help on using the repository browser.