Ignore:
Timestamp:
Oct 29, 2012, 11:06:27 AM (12 years ago)
Author:
aslmd
Message:

UTIL PYTHON. MCD online interface. addressed the TODO list of rev 812, plus other small improvements. this version was sent to beta-testers in the euromars list.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/mcd/mcd.py

    r813 r821  
    7474        self.zonwindtab = None ; self.merwindtab = None ; self.meanvartab = None ; self.extvartab = None
    7575        ## plot stuff
    76         self.xlabel = None ; self.ylabel = None
     76        self.xlabel = None ; self.ylabel = None ; self.title = ""
    7777        self.vertplot = False
    78         self.fmt = "%.2e"
     78        self.fmt = "%.2e"
     79        self.colorm = "jet"
     80        self.fixedlt = False
     81        self.zonmean = False
    7982
    8083    def viking1(self): self.name = "Viking 1 site. MCD v4.3 output" ; self.lat = 22.48 ; self.lon = -49.97 ; self.xdate = 97.
    … …  
    97100        elif self.datekey == 0:  self.title = self.title + " JD " + str(self.xdate) + "."
    98101        if not oneline: self.title = self.title + "\n"
    99         if self.lats is None:  self.title = self.title + " Latitude " + str(self.lat) + "E"
    100         if self.lons is None:  self.title = self.title + " Longitude " + str(self.lon) + "N"
     102        if self.lats is None:  self.title = self.title + " Latitude " + str(self.lat) + "N"
     103        if self.zonmean: self.title = self.title + "Zonal mean over all longitudes."
     104        elif self.lons is None: self.title = self.title + " Longitude " + str(self.lon) + "E"
    101105        if self.xzs is None:   
    102106            self.vertunits()
    103107            self.title = self.title + " Altitude " + str(self.xz) + " " + self.vunits
    104         if self.locts is None: self.title = self.title + " Local time " + str(self.loct) + "h"
     108        if self.locts is None:
     109            self.title = self.title + " Local time " + str(self.loct) + "h"
     110            if not self.fixedlt:  self.title = self.title + " (at longitude 0) "
    105111
    106112    def getextvarlab(self,num):
    … …  
    269275    # print extra MCD variables
    270276        num = self.convertlab(num)
    271         print self.getextvarlab(num) + " ..... " + str(self.extvar[num-1])
     277        dastr = str(self.extvar[num-1])
     278        if dastr == "nan":   print "!!!! There is a problem, probably a value is requested below the surface !!!!"
     279        else:                print self.getextvarlab(num) + " ..... " + dastr
    272280
    273281    def printallextvar(self):
    … …  
    276284
    277285    def htmlprinttabextvar(self,tabtodo):
     286        self.fixedlt = True ## local time is real local time
    278287        self.gettitle()
    279288        print "<hr>"
    … …  
    383392    def diurnal(self,nd=13):
    384393    ### retrieve a local time slice
     394      self.fixedlt = True  ## local time is real local time
    385395      save = self.loct
    386396      self.xlabel = "Local time (Martian hour)"
    … …  
    394404      self.xlabel = "East longitude (degrees)"
    395405      self.prepare(ndx=nd) ; self.ininterv(-180.,180.,nd,start=self.lons,end=self.lone)
    396       for i in range(nd): self.lon = self.xcoord[i] ; self.update() ; self.put1d(i)
     406      if not self.fixedlt: umst = self.loct
     407      for i in range(nd):
     408          self.lon = self.xcoord[i]
     409          if not self.fixedlt: self.loct = (umst + self.lon/15.) % 24
     410          self.update() ; self.put1d(i)
    397411      self.lon = save
    398412
    399413    def meridional(self,nd=19):
    400414    ### retrieve a latitude slice
     415      self.fixedlt = True  ## local time is real local time
    401416      save = self.lat
    402417      self.xlabel = "North latitude (degrees)"
    … …  
    407422    def profile(self,nd=20,tabperso=None):
    408423    ### retrieve an altitude slice (profile)
     424      self.fixedlt = True  ## local time is real local time
    409425      save = self.xz
    410426      self.vertlabel()
    … …  
    482498        yeah.plot(absc,ordo,'-bo') #; mpl.xticks(query.xcoord)
    483499        ax = fig.gca() ; ax.set_ylabel(ordolab) ; ax.set_xlabel(absclab)
    484         if self.zkey == 4: ax.set_yscale('log') ; ax.set_ylim(ax.get_ylim()[::-1])
     500
     501        if self.xzs is not None and self.zkey == 4: ax.set_yscale('log') ; ax.set_ylim(ax.get_ylim()[::-1])
    485502
    486503        if self.lats is not None:      ax.set_xticks(np.arange(-90,91,15)) ; ax.set_xbound(lower=self.lats, upper=self.late)
    487504        elif self.lons is not None:    ax.set_xticks(np.arange(-360,361,30)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
    488505        elif self.locts is not None:   ax.set_xticks(np.arange(0,26,2)) ; ax.set_xbound(lower=self.locts, upper=self.locte)
     506
     507        ax.grid(True, linestyle=':', color='grey')
    489508
    490509      self.gettitle()
    … …  
    500519###################
    501520
    502     def latlon(self,ndx=37,ndy=19,fixedlt=False):
     521    def latlon(self,ndx=37,ndy=19):
    503522    ### retrieve a latitude/longitude slice
    504523    ### default is: local time is not fixed. user-defined local time is at longitude 0.
    … …  
    508527      self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
    509528      self.ininterv(-90.,  90.,ndy,start=self.lats,end=self.late,yaxis=True)
    510       if not fixedlt: umst = self.loct
     529      if not self.fixedlt: umst = self.loct
    511530      for i in range(ndx):
    512531       for j in range(ndy):
    513532         self.lon = self.xcoord[i] ; self.lat = self.ycoord[j]
    514          if not fixedlt: self.loct = (umst + self.lon/15.) % 24
     533         if not self.fixedlt: self.loct = (umst + self.lon/15.) % 24
    515534         self.update() ; self.put2d(i,j)
    516       if not fixedlt: self.loct = umst
     535      if not self.fixedlt: self.loct = umst
    517536      self.lon = save1 ; self.lat = save2 ; self.loct = save3
    518537
    519     def lonalt(self,ndx=37,ndy=20,fixedlt=False):
    520     ### retrieve a longitude/altitude slice
    521       save1 = self.lon ; save2 = self.xz ; save3 = self.loct
     538    def secalt(self,ndx=37,ndy=20,typex="lat"):
     539    ### retrieve a coordinate/altitude slice
     540      save1 = self.lon ; save2 = self.xz ; save3 = self.loct ; save4 = self.lat
     541      self.prepare(ndx=ndx,ndy=ndy)
    522542      self.vertlabel() ; self.ylabel = self.xlabel
    523       self.xlabel = "East longitude (degrees)"
    524       self.prepare(ndx=ndx,ndy=ndy)
    525       self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
    526543      self.vertaxis(ndy,yaxis=True)
    527       if not fixedlt: umst = self.loct
     544      if "lat" in typex:
     545          self.xlabel = "North latitude (degrees)"
     546          self.ininterv(-90.,90.,ndx,start=self.lats,end=self.late)
     547      elif typex == "lon":
     548          self.xlabel = "East longitude (degrees)"
     549          self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
     550      if not self.fixedlt: umst = self.loct
    528551      for i in range(ndx):
    529552       for j in range(ndy):
    530          self.lon = self.xcoord[i] ; self.xz = self.ycoord[j]
    531          if not fixedlt: self.loct = (umst + self.lon/15.) % 24
     553         if typex == "lat":   self.lat = self.xcoord[i]
     554         elif typex == "lon": self.lon = self.xcoord[i]
     555         self.xz = self.ycoord[j]
     556         if not self.fixedlt: self.loct = (umst + self.lon/15.) % 24
    532557         self.update() ; self.put2d(i,j)
    533       if not fixedlt: self.loct = umst
    534       self.lon = save1 ; self.xz = save2 ; self.loct = save3
    535 
    536     def latalt(self,ndx=19,ndy=20,fixedlt=False):
    537     ### retrieve a latitude/altitude slice
    538       save1 = self.lat ; save2 = self.xz ; save3 = self.loct
     558      if not self.fixedlt: self.loct = umst
     559      self.lon = save1 ; self.xz = save2 ; self.loct = save3 ; self.lat = save4
     560
     561    def zonalmean(self,ndx=37,ndy=20,ndmean=32):
     562    ### retrieve a zonalmean lat/altitude slice
     563      self.fixedlt = False
     564      save1 = self.lon ; save2 = self.xz ; save3 = self.loct ; save4 = self.lat
     565      self.prepare(ndx=ndx,ndy=ndy)
    539566      self.vertlabel() ; self.ylabel = self.xlabel
     567      self.vertaxis(ndy,yaxis=True)
    540568      self.xlabel = "North latitude (degrees)"
    541       self.prepare(ndx=ndx,ndy=ndy)
     569      self.ininterv(-180.,180.,ndmean)
     570      coordmean = self.xcoord
    542571      self.ininterv(-90.,90.,ndx,start=self.lats,end=self.late)
    543       self.vertaxis(ndy,yaxis=True)
    544       if not fixedlt: umst = self.loct
     572      umst = self.loct #fixedlt false for this case
     573      for i in range(ndx):
     574       self.lat = self.xcoord[i]
     575       for j in range(ndy):
     576        self.xz = self.ycoord[j]
     577        meanpres = 0. ; meandens = 0. ; meantemp = 0. ; meanzonwind = 0. ; meanmerwind = 0. ; meanmeanvar = np.zeros(5) ; meanextvar = np.zeros(100)       
     578        for m in range(ndmean):
     579           self.lon = coordmean[m]
     580           self.loct = (umst + self.lon/15.) % 24 #fixedlt false for this case
     581           self.update()
     582           meanpres = meanpres + self.pres/float(ndmean) ; meandens = meandens + self.dens/float(ndmean) ; meantemp = meantemp + self.temp/float(ndmean)
     583           meanzonwind = meanzonwind + self.zonwind/float(ndmean) ; meanmerwind = meanmerwind + self.merwind/float(ndmean)
     584           meanmeanvar = meanmeanvar + self.meanvar/float(ndmean) ; meanextvar = meanextvar + self.extvar/float(ndmean)
     585        self.pres=meanpres ; self.dens=meandens ; self.temp=meantemp ; self.zonwind=meanzonwind ; self.merwind=meanmerwind
     586        self.meanvar=meanmeanvar ; self.extvar=meanextvar
     587        self.put2d(i,j)
     588      self.loct = umst #fixedlt false for this case
     589      self.lon = save1 ; self.xz = save2 ; self.loct = save3 ; self.lat = save4
     590
     591    def hovmoller(self,ndtime=25,ndcoord=20,typex="lat"):
     592    ### retrieve a time/other coordinate slice
     593      save1 = self.lat ; save2 = self.xz ; save3 = self.loct ; save4 = self.lon
     594      if typex == "lat":
     595          ndx = ndcoord ; self.xlabel = "North latitude (degrees)"
     596          ndy = ndtime ; self.ylabel = "Local time (Martian hour)"
     597          self.prepare(ndx=ndx,ndy=ndy)
     598          self.ininterv(-90.,90.,ndx,start=self.lats,end=self.late)
     599          self.ininterv(0.,24.,ndy,start=self.locts,end=self.locte,yaxis=True)
     600      elif typex == "lon":
     601          ndx = ndcoord ; self.xlabel = "East longitude (degrees)"
     602          ndy = ndtime ; self.ylabel = "Local time (Martian hour)"
     603          self.prepare(ndx=ndx,ndy=ndy)
     604          self.ininterv(-180.,180.,ndx,start=self.lons,end=self.lone)
     605          self.ininterv(0.,24.,ndy,start=self.locts,end=self.locte,yaxis=True)
     606      elif typex == "alt":
     607          ndy = ndcoord ; self.vertlabel() ; self.ylabel = self.xlabel
     608          ndx = ndtime ; self.xlabel = "Local time (Martian hour)"
     609          self.prepare(ndx=ndx,ndy=ndy)
     610          self.vertaxis(ndy,yaxis=True)
     611          self.ininterv(0.,24.,ndx,start=self.locts,end=self.locte)
    545612      for i in range(ndx):
    546613       for j in range(ndy):
    547          self.lat = self.xcoord[i] ; self.xz = self.ycoord[j]
    548          if not fixedlt: self.loct = (umst + self.lon/15.) % 24
     614         if typex == "lat":   self.lat = self.xcoord[i] ; self.loct = self.ycoord[j]
     615         elif typex == "lon": self.lon = self.xcoord[i] ; self.loct = self.ycoord[j]
     616         elif typex == "alt": self.xz = self.ycoord[j] ; self.loct = self.xcoord[i]
    549617         self.update() ; self.put2d(i,j)
    550       if not fixedlt: self.loct = umst
    551       self.lat = save1 ; self.xz = save2 ; self.loct = save3
     618      self.lat = save1 ; self.xz = save2 ; self.loct = save3 ; self.lon = save4
    552619
    553620    def put2d(self,i,j):
    … …  
    560627      self.extvartab[i,j,1:100] = self.extvar[0:99] ## note: var numbering according to MCD manual is kept
    561628
    562     def makemap2d(self,choice,incwind=False,fixedlt=False,proj="cyl"):
     629    def makemap2d(self,choice,incwind=False,proj="cyl"):
    563630    ### one 2D map is created for the user-defined variable in choice.
    564       self.latlon(fixedlt=fixedlt) ## a map is implicitely a lat-lon plot. otherwise it is a plot (cf. makeplot2d)
     631      self.latlon() ## a map is implicitely a lat-lon plot. otherwise it is a plot (cf. makeplot2d)
    565632      if choice == "wind" or incwind:
    566633          (windx, fieldlabwx) = self.definefield("u")
    … …  
    575642      mpl.figtext(0.5, 0.0, self.ack, ha='center')
    576643
    577     def map2d(self,tabtodo,incwind=False,fixedlt=False,proj="cyl"):
     644    def map2d(self,tabtodo,incwind=False,proj="cyl"):
    578645    ### complete 2D figure with possible multiplots
    579646      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
    … …  
    581648      fig = mpl.figure()
    582649      subv,subh = myplot.definesubplot( len(tabtodo) , fig )
    583       for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1) ; self.makemap2d(tabtodo[i],incwind=incwind,fixedlt=fixedlt,proj=proj)
    584 
    585     def htmlmap2d(self,tabtodo,incwind=False,fixedlt=False,figname="temp.png",title="",back="zMOL"):
     650      for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1) ; self.makemap2d(tabtodo[i],incwind=incwind,proj=proj)
     651
     652    def htmlmap2d(self,tabtodo,incwind=False,figname="temp.png",back="zMOL"):
    586653    ### complete 2D figure with possible multiplots
    587654    ### added in 09/2012 for online MCD
    … …  
    591658      from matplotlib.cm import get_cmap
    592659      from matplotlib import rcParams
    593 
    594       #from mpl_toolkits.basemap import Basemap
    595 
     660      #from mpl_toolkits.basemap import Basemap # does not work
    596661      from Scientific.IO import NetCDF
     662
    597663      filename = "/home/marshttp/surface.nc"
    598664      zefile = NetCDF.NetCDFFile(filename, 'r')
    … …  
    600666      yc = zefile.variables['latitude']
    601667      xc = zefile.variables['longitude']
    602       ## plutot que fieldc = self.getextvar(self.convertlab("topo"))
    603668
    604669      if isinstance(tabtodo,np.str): tabtodo=[tabtodo] ## so that asking one element without [] is possible.
    … …  
    616681        yeah = fig.add_subplot(subv,subh,i+1)
    617682        choice = tabtodo[i]
    618         self.latlon(fixedlt=fixedlt,ndx=64,ndy=48)
     683        self.latlon(ndx=64,ndy=48)
    619684        ## a map is implicitely a lat-lon plot. otherwise it is a plot (cf. makeplot2d)
    620685        (field, fieldlab) = self.definefield(choice)
    621686        if incwind: (windx, fieldlabwx) = self.definefield("u") ; (windy, fieldlabwy) = self.definefield("v")
    622687
    623         proj="moll" ; colorb="jet" ; ndiv=20 ; zeback="molabw" ; trans=1.0 #0.6
    624         title="" ; vecx=None ; vecy=None ; stride=2
     688        proj="moll" ; colorb= self.colorm ; ndiv=20 ; zeback="molabw" ; trans=1.0 #0.6
     689        vecx=None ; vecy=None ; stride=2
    625690        lon = self.xcoord
    626691        lat = self.ycoord
    … …  
    671736      canvas.print_figure(figname, dpi=80)
    672737
    673     def htmlplot2d(self,tabtodo,fixedlt=False,figname="temp.png",title=""):
     738    def htmlplot2d(self,tabtodo,figname="temp.png"):
    674739    ### complete 2D figure with possible multiplots
    675740    ### added in 10/2012 for online MCD
    … …  
    693758        choice = tabtodo[i]
    694759
    695         if self.lons is not None:    self.lonalt(fixedlt=fixedlt,ndx=64,ndy=35)
    696         elif self.lats is not None:  self.latalt(fixedlt=fixedlt,ndx=48,ndy=35)
     760        if self.lons is not None:   
     761           if self.locts is None:  self.secalt(ndx=64,ndy=35,typex="lon")
     762           else:                   self.hovmoller(ndcoord=64,typex="lon")
     763        elif self.lats is not None: 
     764           if self.locts is None: 
     765               if self.zonmean:   self.zonalmean()
     766               else:         self.secalt(ndx=48,ndy=35,typex="lat")
     767           else:                   self.hovmoller(ndcoord=48,typex="lat")
     768        else:
     769           self.hovmoller(ndcoord=35,typex="alt")
    697770
    698771        (field, fieldlab) = self.definefield(choice)
    699772
    700         colorb="jet" ; ndiv=20 ; title=""
     773        colorb=self.colorm ; ndiv=20
    701774
    702775        ## define field. bound field.
    … …  
    714787        ax = fig.gca() ; ax.set_ylabel(self.ylabel) ; ax.set_xlabel(self.xlabel)
    715788
    716         if self.zkey == 4:
     789        if self.lons is not None:   ax.set_xticks(np.arange(-360,361,45)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
     790        elif self.lats is not None: ax.set_xticks(np.arange(-90,91,30)) ; ax.set_xbound(lower=self.lats, upper=self.late)
     791
     792        if self.locts is not None:
     793            if self.xzs is not None: ax.set_xticks(np.arange(0,26,2)) ; ax.set_xbound(lower=self.locts, upper=self.locte)
     794            else:                    ax.set_yticks(np.arange(0,26,2)) ; ax.set_ybound(lower=self.locts, upper=self.locte)
     795
     796        if self.zkey == 4 and self.xzs is not None:
    717797            ax.set_yscale('log') ; ax.set_ylim(ax.get_ylim()[::-1])
    718798        else:
    719799            #ax.set_yticks(np.arange(self.xzs,self.xze,10000.)) ;
    720800            ax.set_ybound(lower=self.xzs, upper=self.xze)
    721 
    722         if self.lons is not None: ax.set_xticks(np.arange(-360,361,45)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
    723         elif self.lats is not None: ax.set_xticks(np.arange(-90,91,30)) ; ax.set_xbound(lower=self.lats, upper=self.late)
    724801
    725802      self.gettitle()
Note: See TracChangeset for help on using the changeset viewer.