Changeset 475 for trunk/UTIL


Ignore:
Timestamp:
Dec 15, 2011, 4:05:49 PM (13 years ago)
Author:
aslmd
Message:

UTIL: Python graphics: 1. corrected locations lat/lon for meso. 2. fixed yaxis settings which did not work in some 1D and 2D cases. 3. moved 1D stuff in imov loop, so the loop is now consistent and generic. Also updated a bit farm_tour to avoid displaying e.g. grep as running job.

Location:
trunk/UTIL
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/UTIL/PYTHON/README.PP

    r453 r475  
    7171Goal: I want to plot averaged results in the file from one time to another time
    7272pp.py -f diagfi.nc -v tsurf --time 4,7
     73
     74Goal: I want to plot a globally-averaged 1D temperature profile
     75pp.py -f diagfi.nc -v temp --time 4 --lat -90,90 --lon -180,180
     76
     77Goal: I want to overplot few globally-averaged 1D temperature profiles at different times
     78pp.py -f diagfi.nc -v temp --time 4 --time 7 --lat -90,90 --lon -180,180
    7379
    7480[only mesoscale for the moment]
  • trunk/UTIL/PYTHON/myplot.py

    r469 r475  
    414414    if typefile in ['mesoapi','meso']:
    415415        [lon2d,lat2d] = getcoord2d(nc)
    416         lon2d = dumpbdy(lon2d,6)
    417         lat2d = dumpbdy(lat2d,6)
    418416    elif typefile in ['gcm']:
    419417        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
     
    958956
    959957## Author: AS
    960 def bidimfind(lon2d,lat2d,vlon,vlat):
     958def bidimfind(lon2d,lat2d,vlon,vlat,file=None):
    961959   import numpy as np
     960   import matplotlib.pyplot as mpl
    962961   if vlat is None:    array = (lon2d - vlon)**2
    963962   elif vlon is None:  array = (lat2d - vlat)**2
     
    965964   idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
    966965   if vlon is not None:
    967        #print lon2d[idy,idx],vlon
    968        if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
     966      if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
    969967   if vlat is not None:
    970        #print lat2d[idy,idx],vlat
    971        if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
    972    return idx,idy
     968      if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
     969   if file is not None:
     970      print idx,idy,lon2d[idy,idx],vlon
     971      print idx,idy,lat2d[idy,idx],vlat
     972      var = file.variables["HGT"][:,:,:]
     973      mpl.contourf(var[0,:,:],30,cmap = mpl.get_cmap(name="Greys_r") ) ; mpl.axis('off') ; mpl.plot(idx,idy,'bo')
     974      mpl.show()
     975   return idy,idx
    973976
    974977## Author: TN
  • trunk/UTIL/PYTHON/planetoplot.py

    r468 r475  
    88### A. Spiga     -- LMD -- 11~12/2011 -- Extended multivar subplot capabilities + cosmetic changes + general cleaning and tests
    99### A. Colaitis  -- LMD --    12/2011 -- Added movie capability [mencoder must be installed]
    10 ### A. Spiga     -- LMD --    12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...]
     10### A. Spiga     -- LMD --    12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...] + consistent generic movie loop
    1111
    1212def planetoplot (namefiles,\
     
    154154      elif typefile in ['meso','mesoapi','geo','mesoideal']:
    155155          ### the following lines are kind of dirty... not possible to ask for several lats and lons
    156           if vlon is not None or vlat is not None:   indices = bidimfind(lon2d,lat2d,vlon,vlat)
    157           if slon is not None: slon[0][0] = indices[0] ; slon[0][1] = indices[0]
    158           if slat is not None: slat[0][0] = indices[1] ; slat[0][1] = indices[1]
     156          if vlon is not None or vlat is not None:   
     157              indices = bidimfind(lon2d,lat2d,vlon,vlat) #,file=nc)  #placer un point sur carte
     158              lonp,latp = ( lon2d[indices[0],indices[1]] , lat2d[indices[0],indices[1]] )
     159          if typefile in ['mesoapi','meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  ### important to do that now and not before
     160          if slon is not None: slon[0][0] = indices[1] ; slon[0][1] = indices[1]  #...this is idx
     161          if slat is not None: slat[0][0] = indices[0] ; slat[0][1] = indices[0]  #...this is idy
    159162          if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
    160163          else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
     
    236239                tabtime = all_time[0];tab = all_var[0];k = 1
    237240                if var2: tab2 = all_var2[0]
    238                 while k != len(namefiles):
     241                while k != len(namefiles) and len(all_time[k]) != 0:
    239242                    if var2: tab2 = np.append(tab2,all_var2[k],axis=0)
    240243                    tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
     
    292295       print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
    293296       ####################################################################
     297       ########## REDUCE FIELDS
     298       ####################################################################
    294299       error = False
    295300       varname = all_varname[index_f]
     
    308313       changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1)  ## default for 1D plots: superimposed. to be reworked for better flexibility.
    309314       if changesubplot: subplot(subv,subh,nplot)
    310        ### Map projection                   
    311        if mapmode == 1:     m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ; x, y = m(lon2d, lat2d)
    312        elif mapmode ==0:    m = None ; x = None ; y = None
    313315       ####################################################################
    314 
    315        if not error:
     316       if error:
     317               errormess("There is an error in reducing field !")
     318       else:
    316319               ticks = ndiv + 1
    317320               fvar = varname
    318321               if anomaly: fvar = 'anomaly'
    319                if mapmode == 0:
     322               ###
     323               if mapmode == 0:    ### could this be moved inside imov loop ?
    320324                   itime=indextime
    321325                   if len(what_I_plot.shape) is 3:itime=[0]
     326                   m = None ; x = None ; y = None
    322327                   what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    323328                         itime,what_I_plot, len(all_var[index_f].shape),vertmode)
    324                    zxmin, zxmax = xaxis ; zymin, zymax = yaxis
    325                    if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
    326                    if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
    327                    if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
    328                    if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
    329                    if invert_y:     lima,limb = mpl.pyplot.ylim() ; mpl.pyplot.ylim(limb,lima)
    330                    if ylog:         mpl.pyplot.semilogy()
    331 
    332                if (fileref is not None) and (index_f is numplot-1):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
     329               ###
     330               if (fileref is not None) and (index_f == numplot-1):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
    333331               else:                                                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
     332               if (fileref is not None) and (index_f == numplot-1):    colorb = "RdBu_r"
    334333               if colorb in ["def","nobar"]:                           palette = get_cmap(name=defcolorb(fvar.upper()))
    335                elif (fileref is not None) and (index_f is numplot-1):  palette = get_cmap(name="RdBu_r")
    336334               else:                                                   palette = get_cmap(name=colorb)
    337 
    338                ##### simple 2D field and movies of 2D fields
    339                if len(what_I_plot.shape) >= 2:
    340                  if (len(what_I_plot.shape) is 3 and mrate is None):  errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
     335               ##### 1. ELIMINATE >3D CASES
     336               if len(what_I_plot.shape) >= 4:
     337                 print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
     338                 errormess("Are you sure you did not forget to prescribe a dimension ?")
     339               ##### 2. HANDLE simple 1D/2D field and movies of 1D/2D fields
     340               else:
    341341                 if mrate is not None: iend=len(time)-1
    342342                 else:                 iend=0
    343343                 imov = 0
    344                  if var2:  which = "contour" ## have to start with contours rather than shading
    345                  else:     which = "regular"
     344                 if len(what_I_plot.shape) == 3:
     345                    if var2:               which = "contour" ## have to start with contours rather than shading
     346                    else:                  which = "regular"
     347                    if mrate is None:      errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
     348                 elif len(what_I_plot.shape) == 2:
     349                    if var2:               which = "contour" ## have to start with contours rather than shading
     350                    else:                  which = "regular"
     351                    if mrate is not None:  which = "unidim"
     352                 elif len(what_I_plot.shape) == 1:
     353                    which = "unidim"
     354                 ##### IMOV LOOP #### IMOV LOOP
    346355                 while imov <= iend:
    347356                    print "-> frame ",imov+1, which
     
    355364                        if mrate is None or what_I_plot_contour.ndim < 3:   what_I_plot_frame = what_I_plot_contour
    356365                        else:                                               what_I_plot_frame = what_I_plot_contour[imov,:,:]
    357                     if mrate is not None:     
    358                         if mapmode == 1:
     366                    elif which == "unidim":
     367                        if mrate is None:                                   what_I_plot_frame = what_I_plot
     368                        else:                                               what_I_plot_frame = what_I_plot[:,imov]  ## because swapaxes
     369                    #if mrate is not None:     
     370                    if mapmode == 1:
    359371                            m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)  ## this is dirty, defined above but out of imov loop
    360372                            x, y = m(lon2d, lat2d)                                         ## this is dirty, defined above but out of imov loop
     
    362374#                   if typefile in ['mesoideal']:    what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
    363375
    364                     if imov >= 0:
    365                         # Renew axis directives for movie frames which are not the first one.
    366                         zxmin, zxmax = xaxis ; zymin, zymax = yaxis
    367                         if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
    368                         if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
    369                         if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
    370                         if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
    371                         if invert_y:     lima,limb = mpl.pyplot.ylim() ; mpl.pyplot.ylim(limb,lima)
    372                         if ylog:         mpl.pyplot.semilogy()
    373                    
    374                     if which == "regular":
     376                    if which == "unidim":
     377                        lbl = ""
     378                        if indexlat is not None:  lbl = lbl + " ix" + str(indexlat[0])
     379                        if indexlon is not None:  lbl = lbl + " iy" + str(indexlon[0])
     380                        if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0])
     381                        if indextime is not None: lbl = lbl + " it" + str(indextime[0])
     382                        if mrate is not None: x = y  ## because swapaxes...
     383                        if indexvert is not None or indextime is None:    plot(x,what_I_plot_frame,label=lbl)  ## regular plot
     384                        else:                                             plot(what_I_plot_frame,x,label=lbl)  ## vertical profile
     385                        if nplot > 1: legend(loc='best')
     386                        if indextime is None and axtime is not None:      xlabel(axtime.upper()) ## define the right label
     387                        if save == 'txt':  writeascii(np.transpose(what_I_plot),'profile'+str(nplot*1000+imov)+'.txt')
     388
     389                    elif which == "regular":
    375390                        if hole:         what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
    376391                        else:            what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
     
    384399                            if mapmode == 1:       m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
    385400                            elif mapmode == 0:     pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
     401
    386402                        if colorb != 'nobar':       
    387                             if (fileref is not None) and (index_f is numplot-1):   daformat = "%.3f"
     403                            if (fileref is not None) and (index_f == numplot-1):   daformat = "%.3f"
    388404                            elif mult != 1:                                        daformat = "%.1f"
    389405                            else:                                                  daformat = fmtvar(fvar.upper())
     
    413429                        elif mapmode == 1:  cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
    414430
    415 
    416                     if which == "regular":
     431                    if which in ["regular","unidim"]:
     432
     433                        if nplot > 1 and which == "unidim":
     434                           pass  ## because we superimpose nplot instances
     435                        else:
     436                           # Axis directives for movie frames [including the first one).
     437                           zxmin, zxmax = xaxis ; zymin, zymax = yaxis
     438                           if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
     439                           if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
     440                           if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
     441                           if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
     442                           if ylog:      mpl.pyplot.semilogy()
     443                           if invert_y:  ax = mpl.pyplot.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
     444
    417445                        if mrate is not None:
    418446                           ### THIS IS A MENCODER MOVIE
     
    437465                                 myfile.write("last_image = "+str(iend)+";"+ '\n')
    438466                                 myfile.close()
    439                         if var2:  which = "contour"
     467                        if var2 and which == "regular":  which = "contour"
    440468                        imov = imov+1
    441469                    elif which == "contour":
    442470                        which = "regular"
    443 
    444                ##### 1D field
    445                elif len(what_I_plot.shape) is 1:
    446                  lbl = ""
    447                  if indexlat is not None:  lbl = lbl + " ix" + str(indexlat)
    448                  if indexlon is not None:  lbl = lbl + " iy" + str(indexlon)
    449                  if indexvert is not None: lbl = lbl + " iz" + str(indexvert)
    450                  if indextime is not None: lbl = lbl + " it" + str(indextime)
    451 
    452                  if indexvert is not None or indextime is None:    plot(x,what_I_plot,label=lbl)  ## regular plot
    453                  else:                                             plot(what_I_plot,x,label=lbl)  ## vertical profile
    454                  if nplot > 1: legend(loc='best')
    455                  if indextime is None and axtime is not None:      xlabel(axtime.upper()) ## define the right label
    456                  if save == 'txt':  writeascii(np.transpose(what_I_plot),'profile'+str(nplot)+'.txt')
    457 
    458                #### Other cases: (maybe plot 3-D field one day ??)
    459                else:
    460                  print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
    461                  errormess("Are you sure you did not forget to prescribe a dimension ?")
    462        else:
    463                errormess("There is an error in reducing field !")
    464471
    465472       ### Next subplot
     
    469476       if mrate is not None: basename = "movie_" + basename
    470477       if typefile in ['mesoapi','meso']:
    471             if slon is not None: basename = basename + "_lon_" + str(int(lon2d[indices[1],indices[0]]))
    472             if slat is not None: basename = basename + "_lat_" + str(int(lat2d[indices[1],indices[0]]))
     478            if slon is not None: basename = basename + "_lon_" + str(int(round(lonp)))
     479            if slat is not None: basename = basename + "_lat_" + str(int(round(latp)))
    473480            plottitle = basename
    474481            ### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus)
  • trunk/UTIL/PYTHON/pp.py

    r468 r475  
    107107          #####
    108108          elif typefile == "gcm":
     109            inputvar = zevars
     110            if opt.var2 is not None : inputvar = np.append(inputvar,opt.var2)
    109111            interpolated_files=""
    110112            interpolated_files=call_zrecast(interp_mode=opt.itp,\
    111113                    input_name=zefiles,\
    112                     fields=zevars,\
     114                    fields=inputvar,\
    113115                    predefined=opt.intas)
    114116
  • trunk/UTIL/farm_tour

    r454 r475  
    1717machine='levan'
    1818#####################
    19 ssh $machine -t "ps --sort=-pcpu -e -o '%C %u %c %x %p' -G lmdjus r | grep -v ps | grep -v root >> ~/log1" > /dev/null 2> /dev/null
     19ssh $machine -t "ps --sort=-pcpu -e -o '%C %u %c %x %p' -G lmdjus r | grep -v ps | grep -v root | grep -v grep | grep -v bash >> ~/log1" > /dev/null 2> /dev/null
    2020echo "*********************************************************" >> ~/log
    2121echo " $machine ... "`head -n -1 ~/log1 | wc -l`" processeurs utilises sur 8"  >> ~/log
     
    2727machine='penn'
    2828#####################
    29 ssh $machine -t "ps --sort=-pcpu -e -o '%C %u %c %x %p' -G lmdjus r | grep -v ps | grep -v root >> ~/log2" > /dev/null 2> /dev/null
     29ssh $machine -t "ps --sort=-pcpu -e -o '%C %u %c %x %p' -G lmdjus r | grep -v ps | grep -v root | grep -v grep | grep -v bash >> ~/log2" > /dev/null 2> /dev/null
    3030echo "*********************************************************" >> ~/log
    3131echo " $machine ... "`head -n -1 ~/log2 | wc -l`" processeurs utilises sur 8"  >> ~/log
     
    3737machine='viccaro'
    3838#####################
    39 ssh $machine -t "ps --sort=-pcpu -e -o '%C %u %c %x %p' -G lmdjus r | grep -v ps | grep -v root >> ~/log3" > /dev/null 2> /dev/null
     39ssh $machine -t "ps --sort=-pcpu -e -o '%C %u %c %x %p' -G lmdjus r | grep -v ps | grep -v root | grep -v grep | grep -v bash >> ~/log3" > /dev/null 2> /dev/null
    4040echo "*********************************************************" >> ~/log
    4141echo " $machine ... "`head -n -1 ~/log3 | wc -l`" processeurs utilises sur 8"  >> ~/log
Note: See TracChangeset for help on using the changeset viewer.