Changeset 489

Timestamp:

Dec 21, 2011, 3:22:19 PM (13 years ago)

Author:

aslmd

Message:

UTIL PYTHON : 1) no more negative values for --time but use --axtime lt ; this actually fix a stupid bug with negative lat lon. 2) use of -l is now possible with GCM files, try it with -i 2. 3) now correct handling of sol/ls in mesoscale files 4) for mesoscale files with chosen lat+lon now shows a topo map with a cross on it.

Location:

trunk/UTIL/PYTHON

Files:

• README.PP (modified) (1 diff)
• gcm_transformations.py (modified) (3 diffs)
• myplot.py (modified) (5 diffs)
• myscript.py (modified) (2 diffs)
• planetoplot.py (modified) (10 diffs)
• pp.py (modified) (3 diffs)

trunk/UTIL/PYTHON/README.PP

@@ -84 +84 @@
 [only mesoscale for the moment]
 Goal: I want to plot results for two different LOCAL times in the file next to one another
-pp.py -f wrfout**** -v TSURF --time -4 -- time -7
+pp.py -f wrfout**** -v TSURF --time 4 -- time 7 --axtime lt
 
 ***********************************************************************************

trunk/UTIL/PYTHON/gcm_transformations.py

@@ -9 +9 @@
                     input_name      = None, \
                     fields  = 'all', \
+                    limites = None, \
                     predefined = None):
 
@@ -21 +22 @@
     outputfilename=""
     f = open('zrecast.auto.def', 'w')
+
     for zfile in input_name:
         f.write(zfile+"\n")
@@ -52 +54 @@
                 f.write("1"+"\n")
                 f.write("yes"+"\n")
-                f.write("370 0.1"+"\n")
+                if limites[0] != -9999.:  f.write(str(limites[0])+" "+str(limites[-1])+"\n")
+                else:                     f.write("370 0.1"+"\n")
+                #f.write("1000000 100"+"\n")
                 f.write("20"+"\n")
         else:

trunk/UTIL/PYTHON/myplot.py

@@ -284 +284 @@
     nc  = Dataset(namefile)
     zetime = None
+    days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
+    plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
     if 'Times' in nc.variables:
         zetime = nc.variables['Times'][0]
@@ -290 +292 @@
     if zetime is not None \
        and 'vert' not in nc.variables:
-        #### 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]) - 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)
+        ##### strangely enough this does not work for api or ncrcat results!
+        zesol = plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
+        dals = int( 10. * sol2ls ( zesol ) ) / 10.
         ###
         zetime2 = nc.variables['Times'][1]
@@ -302 +301 @@
         zehour    = one
         zehourin  = abs ( next - one )
-        ### A CORRIGER POUR LES SOLS IL FAUT LIRE LE MONTH !!!!
-        if getaxis:
+        if not getaxis:
+            lschar = "_Ls"+str(dals)
+        else:
             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
+               zetime = nc.variables['Times'][iii] 
                ltaxis[iii] = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
-               solaxis[iii] = getattr( nc, 'JULDAY' ) + zeday +  ltaxis[iii] / 24.
+               solaxis[iii] = ltaxis[iii] / 24. + plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
                lsaxis[iii] = sol2ls ( solaxis[iii] )
                if ltaxis[iii] < ltaxis[iii-1]: ltaxis[iii] = ltaxis[iii] + 24.
-               #print ltaxis[iii], solaxis[iii], lsaxis[iii]
+               #print ltaxis[iii], solaxis[iii], lsaxis[iii], getattr( nc, 'JULDAY' )
             lschar = lsaxis ; zehour = solaxis ; zehourin = ltaxis
     else:
@@ -982 +981 @@
       print idx,idy,lat2d[idy,idx],vlat
       var = file.variables["HGT"][:,:,:]
-      mpl.contourf(var[0,:,:],30,cmap = mpl.get_cmap(name="Greys_r") ) ; mpl.axis('off') ; mpl.plot(idx,idy,'bo')
+      mpl.contourf(var[0,:,:],30,cmap = mpl.get_cmap(name="Greys_r") ) ; mpl.axis('off') ; mpl.plot(idx,idy,'mx',mew=4.0,ms=20.0)
       mpl.show()
    return idy,idx
@@ -991 +990 @@
 # output : all indexes to be taken into account for reducing field
   import numpy as np
-  ### added by AS to treat the case of stime = - LT
-  if saxis is not None:
-      if saxis[0][0] < 0: saxis = - saxis
-  ###
   if ( np.array(axis).ndim == 2):
       axis = axis[:,0]

trunk/UTIL/PYTHON/myscript.py

@@ -50 +50 @@
     parser.add_option('--vert',         action='append',dest='svert',  type="string",  default=None, help='slices along vert. 2 comma-separated values: averaging') 
     parser.add_option('--column',       action='store_true',dest='column',             default=False,help='changes --vert z1,z2 from MEAN to INTEGRATE along z')
-    parser.add_option('--time',         action='append',dest='stime',  type="string",  default=None, help='slices along time. 2 comma-separated values: averaging. negative: local time [meso].') 
+    parser.add_option('--time',         action='append',dest='stime',  type="string",  default=None, help='slices along time. 2 comma-separated values: averaging') 
     parser.add_option('--xmax',         action='store',dest='xmax',    type="float",   default=None, help='max value for x-axis in contour-plots [max(xaxis)]')
     parser.add_option('--ymax',         action='store',dest='ymax',    type="float",   default=None, help='max value for y-axis in contour-plots [max(yaxis)]')
@@ -57 +57 @@
     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('--axtime',       action='store',dest='axtime',  type="string",  default=None, help='choose "ls" or "sol" for time reference')
+    parser.add_option('--axtime',       action='store',dest='axtime',  type="string",  default=None, help='choose "ls","sol","lt" for time ref (1D or --time)')
 
     ### OPERATIONS BETWEEN FILES

trunk/UTIL/PYTHON/planetoplot.py

@@ -74 +74 @@
     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, xlabel
+    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel, axis
     from matplotlib.cm import get_cmap
     import numpy as np
@@ -147 +147 @@
       if ((proj == None) and (typefile not in ['mesoideal'])):   proj = getproj(nc)                 
 
-##########################################################
+      if firstfile:
+         ##########################
+         ### Define plot boundaries
+         ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
+         if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
+         elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
+         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
+         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom)
+         elif zarea is not None: [wlon,wlat] = latinterv(area=zarea)
+
+########################################################## LOAD 4D DIMENSIONS : x, y, z, t
       if typefile == "gcm":
           lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:]
@@ -158 +168 @@
           ### 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) #,file=nc)  #placer un point sur carte
+              if firstfile:  iwantawhereplot = nc
+              else:          iwantawhereplot = None
+              indices = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)  #placer un point sur carte
               lonp,latp = ( lon2d[indices[0],indices[1]] , lat2d[indices[0],indices[1]] )
           if typefile in ['mesoapi','meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  ### important to do that now and not before
@@ -183 +195 @@
               elif "Time" in nc.variables:  time = count + np.arange(0,len(nc.variables["Time"]),1)
               else:                         time = count + np.arange(0,1,1)
-              count = time[-1] + 1  ## so that a cat is possible with simple subscripts
+              if nnn > 0:  count = time[-1] + 1  ## so that a cat is possible with simple subscripts
+              else:        count = 0
+          if axtime in ["lt"]:
+              for i in range(len(time)):  time[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
+              print "LOCAL TIMES.... ", time
           ###
           if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
@@ -198 +214 @@
        ## Faire d'autre checks sur les compatibilites entre fichiers!!
 ##########################################################
-
-      if firstfile:
-         ##########################
-         ### Define plot boundaries
-         ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
-         if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
-         elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
-         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
-         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom)
-         elif zarea is not None: [wlon,wlat] = latinterv(area=zarea) 
 
       all_varname[k] = varname
@@ -246 +252 @@
                     if var2: tab2 = np.append(tab2,all_var2[k],axis=0) 
                     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)
@@ -287 +292 @@
        else:                            yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
        time = all_time[index_f]
-       if stime is not None:
-           if stime[0][0] < 0:
-               if typefile in ['mesoapi','meso']:
-                   for i in range(len(time)):  time[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
-                   print "OK... WORKING WITH LOCAL TIMES"
-               else: errormess("local times not supported. not too hard to modify the code though.")
        if mrate is not None:                 indextime = None 
        else:                                 indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
@@ -298 +297 @@
        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:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
+       #var = nc.variables["phisinit"][:,:]
+       #contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)
+       #show()
+       #exit()
        ####################################################################
        ########## REDUCE FIELDS
@@ -337 +340 @@
                if colorb in ["def","nobar"]:                           palette = get_cmap(name=defcolorb(fvar.upper()))
                else:                                                   palette = get_cmap(name=colorb)
-               ##### 1. ELIMINATE >3D CASES
-               if len(what_I_plot.shape) >= 4:
+               ##### 1. ELIMINATE 0D or >3D CASES
+               if len(what_I_plot.shape) == 0:   
+                 print "VALUE VALUE VALUE VALUE ::: ",what_I_plot
+                 save = 'donothing'
+               elif len(what_I_plot.shape) >= 4:
                  print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
                  errormess("Are you sure you did not forget to prescribe a dimension ?")
@@ -542 +548 @@
         elif save in ['eps','svg','pdf']:     makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)
         elif save == 'gui':                   show()
+        elif save == 'donothing':             pass
         elif save == 'txt':                   print "Saved results in txt file." 
         else: 

trunk/UTIL/PYTHON/pp.py

@@ -44 +44 @@
         yeah = glob.glob(opt.file[0]) ; yeah.sort()
         opt.file[0] = yeah[0] 
-        for file in yeah: opt.file[0] = opt.file[0] + "," + file
+        for file in yeah[1:]: opt.file[0] = opt.file[0] + "," + file
 
     #############################
@@ -79 +79 @@
             zelevel = float(inputnvert[0])
             ze_interp_levels = [-9999.]
+        elif np.array(inputnvert).size == 2:
+            zelevel = -99.
+            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),20)
         else:
             zelevel = -99.
@@ -122 +125 @@
                     input_name=zefiles,\
                     fields=inputvar,\
+                    limites = ze_interp_levels,\
                     predefined=opt.intas)