Changeset 468

Timestamp:

Dec 9, 2011, 8:33:29 PM (13 years ago)

Author:

aslmd

Message:

UTIL:PYTHON: added support for Ls and sol indexing through option --axtime (mesoscale model files). also a bit of cleaning, especially with info to user, now more compact.

Location:

trunk/UTIL/PYTHON

Files:

• myplot.py (modified) (9 diffs)
• myscript.py (modified) (1 diff)
• planetoplot.py (modified) (15 diffs)
• pp.py (modified) (2 diffs)

trunk/UTIL/PYTHON/myplot.py

@@ -20 +20 @@
 
 ## Author: AS
-def getname(var=False,winds=False,anomaly=False):
+def getname(var=False,var2=False,winds=False,anomaly=False):
     if var and winds:     basename = var + '_UV'
     elif var:             basename = var
@@ -26 +26 @@
     else:                 errormess("please set at least winds or var",printvar=nc.variables)
     if anomaly:           basename = 'd' + basename
+    if var2:              basename = basename + '_' + var2
     return basename
 
@@ -84 +85 @@
     shape = np.array(input).shape
     #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
-    print 'IN  REDUCEFIELD dim,shape: ',dimension,shape
     if anomaly: print 'ANOMALY ANOMALY'
     output = input
@@ -171 +171 @@
         elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
         elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
-    dimension = np.array(output).ndim
-    shape = np.array(output).shape
-    print 'OUT REDUCEFIELD dim,shape: ',dimension,shape
+    dimension2 = np.array(output).ndim
+    shape2 = np.array(output).shape
+    print 'REDUCEFIELD dim,shape: ',dimension,shape,' >>> ',dimension2,shape2 
     return output, error
 
@@ -264 +264 @@
 
 ## Author: AS
-def getlschar ( namefile ):
+def getlschar ( namefile, getaxis=False ):
     from netCDF4 import Dataset
     from timestuff import sol2ls
@@ -281 +281 @@
         #### strangely enough this does not work for api or ncrcat results!
         zetimestart = getattr(nc, 'START_DATE')
-        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
-        if zeday < 0:    lschar=""  ## might have crossed a month... fix soon
-        else:            lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. )
+        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) - 1 ###?
+        if zeday < 0:    dals = -99    ## might have crossed a month... fix soon
+        else:            dals = int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10.
+        if not getaxis:  lschar = "_Ls"+str(dals)
         ###
         zetime2 = nc.variables['Times'][1]
@@ -290 +291 @@
         zehour    = one
         zehourin  = abs ( next - one )
+        ### A CORRIGER POUR LES SOLS IL FAUT LIRE LE MONTH !!!!
+        if getaxis:
+            zelen = len(nc.variables['Times'][:])
+            yeye = range(zelen) ; lsaxis = range(zelen) ; solaxis = range(zelen) ; ltaxis = range(zelen)
+            for iii in yeye:
+               zetime = nc.variables['Times'][iii] ; zetimestart = getattr(nc, 'START_DATE')
+               zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) - 1
+               ltaxis[iii] = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
+               solaxis[iii] = getattr( nc, 'JULDAY' ) + zeday +  ltaxis[iii] / 24.
+               lsaxis[iii] = sol2ls ( solaxis[iii] )
+               if ltaxis[iii] < ltaxis[iii-1]: ltaxis[iii] = ltaxis[iii] + 24.
+               #print ltaxis[iii], solaxis[iii], lsaxis[iii]
+            lschar = lsaxis ; zehour = solaxis ; zehourin = ltaxis
     else:
         lschar=""
@@ -649 +663 @@
     ###
     if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
-    print "BOUNDS field ", min(fieldcalc), max(fieldcalc)
-    print "BOUNDS adopted ", zevmin, zevmax
+    print "BOUNDS field ", min(fieldcalc), max(fieldcalc), " //// adopted", zevmin, zevmax
     return zevmin, zevmax
 
@@ -1021 +1034 @@
    x = array(x)
    y = array(y)
-   print "CHECK: what_I_plot.shape", what_I_plot.shape
-   print "CHECK: x.shape, y.shape", x.shape, y.shape
+   print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape
    if len(shape) == 1:
       if shape[0] != len(x):

trunk/UTIL/PYTHON/myscript.py

@@ -53 +53 @@
     parser.add_option('--ymin',         action='store',dest='ymin',    type="float",   default=None, help='min value for y-axis in contour-plots [min(yaxis)]')
     parser.add_option('--inverty',      action='store_true',dest='inverty',            default=False,help='force decreasing values along y-axis (e.g. p-levels)')
-    parser.add_option('--logy',         action='store_true',dest='logy',               default=False,help='set y-axis to logarithmic.')   
+    parser.add_option('--logy',         action='store_true',dest='logy',               default=False,help='set y-axis to logarithmic')   
+    parser.add_option('--axtime',       action='store',dest='axtime',  type="string",  default=None, help='choose "ls" or "sol" for time reference')
 
     ### OPERATIONS BETWEEN FILES

trunk/UTIL/PYTHON/planetoplot.py

@@ -57 +57 @@
            mquality=False,\
            trans=1,\
-           zarea=None):
+           zarea=None,\
+           axtime=None):
 
 
@@ -72 +73 @@
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
-    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend
+    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel
     from matplotlib.cm import get_cmap
     import numpy as np
@@ -113 +114 @@
     ### Loop over the files + vars initially separated by commas to be plotted on the same figure ###
     #################################################################################################
-    k = 0 ; firstfile = True
+    k = 0 ; firstfile = True ; count = 0
     for nnn in range(len(namefiles)):
      for vvv in range(len(var)): 
 
-      print "********** LOOP..... THIS IS SUBPLOT NUMBER.....",k
-
       ######################
       ### Load NETCDF object
-      namefile = namefiles[nnn] ; print "********** THE NAMEFILE IS....", namefile 
+      namefile = namefiles[nnn] 
       nc  = Dataset(namefile)
 
@@ -134 +133 @@
       ### ... VAR
       varname=var[vvv]
-      print "********** THE VAR IS....",varname, var2
       if varname not in nc.variables: varname = False
       ### ... WINDS
@@ -153 +151 @@
           else:                           errormess("no time axis found.")
           vert = nc.variables["altitude"][:]
+          if axtime in ["ls","sol"]:   errormess("not supported. should not be too difficult though.")
       elif typefile in ['meso','mesoapi','geo','mesoideal']:
           ### the following lines are kind of dirty... not possible to ask for several lats and lons
-          if vlon is not None or vlat is not None:   indices = bidimfind(lon2d,lat2d,vlon,vlat) ; print '********** INDICES: ', indices
+          if vlon is not None or vlat is not None:   indices = bidimfind(lon2d,lat2d,vlon,vlat) 
           if slon is not None: slon[0][0] = indices[0] ; slon[0][1] = indices[0]
           if slat is not None: slat[0][0] = indices[1] ; slat[0][1] = indices[1]
@@ -169 +168 @@
           else:                 londim='WEST-EAST_PATCH_END_UNSTAG'
           lon = np.arange(0,getattr(nc,londim),1) ; lat = np.arange(0,getattr(nc,latdim),1)
-          if "Times" in nc.variables:time = np.arange(0,len(nc.variables["Times"]),1)
-          elif "Time" in nc.variables:time = np.arange(0,len(nc.variables["Time"]),1)
+          ###
+          if axtime in ["ls","sol"]:
+              lstab, soltab, lttab = getlschar ( namefile, getaxis = True )
+              if axtime == "ls":      time = lstab 
+              elif axtime == "sol":   time = soltab
+          else:
+              if "Times" in nc.variables:   time = count + np.arange(0,len(nc.variables["Times"]),1)
+              elif "Time" in nc.variables:  time = count + np.arange(0,len(nc.variables["Time"]),1)
+              count = time[-1] + 1  ## so that a cat is possible with simple subscripts
+          ###
           if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
           else:
@@ -209 +216 @@
       ##### GENERAL STUFF HERE
           all_var[k] = getfield(nc,varname)
-      print "********** all_var[k].shape", all_var[k].shape
-      k += 1
-      firstfile = False
+      
+      print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
       #### End of for namefile in namefiles
 
@@ -217 +223 @@
     ### Operation on files
     if ope is not None:
-        print ope
+        print "********** OPERATION: ",ope
         if "var" not in ope:
              if len(var) > 1: errormess("for this operation... please set only one var !")
@@ -228 +234 @@
                 all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; numplot = numplot+2
              elif ope in ["cat"]:
-                tab = all_var[0];k = 1
+                tabtime = all_time[0];tab = all_var[0];k = 1
                 if var2: tab2 = all_var2[0]
                 while k != len(namefiles):
                     if var2: tab2 = np.append(tab2,all_var2[k],axis=0) 
-                    tab = np.append(tab,all_var[k],axis=0) ; k += 1
-                all_time[0] = np.arange(0,len(tab),1) ### AS: time reference is too simplistic, should be better
-                all_var[0] = np.array(tab) ; numplot = 1
+                    tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
+                #all_time[0] = np.arange(0,len(tab),1) ### AS: time reference is too simplistic, should be better
+                all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
                 if var2: all_var2[0] = np.array(tab2)
              else: errormess(ope+" : non-implemented operation. Check pp.py --help")
@@ -284 +290 @@
        ltst = None 
        if typefile in ['mesoapi','meso'] and indextime is not None:  ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) )  
-       print "********** index lon, lat, vert, time ",indexlon,indexlat,indexvert,indextime
+       print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
        ####################################################################
        error = False
@@ -380 +386 @@
                         if colorb != 'nobar':        
                             if (fileref is not None) and (index_f is numplot-1):   daformat = "%.3f"
+                            elif mult != 1:                                        daformat = "%.1f"
                             else:                                                  daformat = fmtvar(fvar.upper())
                             colorbar( fraction=0.05,pad=0.03,format=daformat,\
-                                      ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),extend='neither',spacing='proportional' ) 
+                                      ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' ) 
                         if winds:
                             if typefile in ['mesoapi','meso']:
@@ -443 +450 @@
                  if indextime is not None: lbl = lbl + " it" + str(indextime)
 
-                 if indexvert is not None:    plot(x,what_I_plot,label=lbl)  ## regular plot
-                 else:                        plot(what_I_plot,x,label=lbl)  ## vertical profile
-                 legend(loc='best')
+                 if indexvert is not None or indextime is None:    plot(x,what_I_plot,label=lbl)  ## regular plot
+                 else:                                             plot(what_I_plot,x,label=lbl)  ## vertical profile
+                 if nplot > 1: legend(loc='best')
+                 if indextime is None and axtime is not None:      xlabel(axtime.upper()) ## define the right label
                  if save == 'txt':  writeascii(np.transpose(what_I_plot),'profile'+str(nplot)+'.txt')
 
@@ -456 +464 @@
 
        ### Next subplot
-       basename = getname(var=varname,winds=winds,anomaly=anomaly)
+       basename = getname(var=varname,var2=var2,winds=winds,anomaly=anomaly)
        if len(what_I_plot.shape) > 3:
            basename = basename + getstralt(nc,level) 
@@ -464 +472 @@
             if slat is not None: basename = basename + "_lat_" + str(int(lat2d[indices[1],indices[0]]))
             plottitle = basename
-            if addchar:   [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] )  ;  plottitle = plottitle + addchar
+            ### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus)
+            if addchar and indextime is not None:   [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] )  ;  plottitle = plottitle + addchar
             if ltst is not None and ( (mapmode == 0) or (proj in ["lcc","laea","merc","nsper"]) ):  plottitle = plottitle + "_LT" + str(ltst)
        else:
@@ -472 +481 @@
                 else:                        plottitle = basename+' '+namefiles[0]#index_f]
             else:                            plottitle = basename+' '+namefiles[0]#index_f]
-       if mult != 1:                         plottitle = str(mult) + "*" + plottitle
+       if mult != 1:                         plottitle = '{:.0e}'.format(mult) + "*" + plottitle
        if zetitle != "fill":                 
           plottitle = zetitle

trunk/UTIL/PYTHON/pp.py

@@ -141 +141 @@
                 invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
                 blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,\
-                mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area)
+                mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime)
         print 'DONE: '+name
         system("rm -f to_be_erased")
@@ -156 +156 @@
     f.write(command)
 
-    print "********** OPTIONS: ", opt
+    #print "********** OPTIONS: ", opt
     print "********************************************************** END"