follow-up to bin averages and looping

[Chris Weaver <weaver@xxxxxxxxxxxxxxxxxxxxxxxxxxx>]

Hi,

Steve Hankin graciously provided me with some solutions to my general problem of 
creating averages of one variable in "bins" of another, so I thought I would 
post a summary here in case anyone was interested.

Let's say we want to bin sst by slp and then average the result over a given XY 
region (e.g., calculate mean sst in each 2 mb slp bin from 991 to 1041 mb).  The 
key point is we need to define an axis over which the binning will take 
place--in this case the axis units are "slp":

  use coads_climatology
  define axis/X=991:1041:2/units=slp Xslp

The limits of each slp bin (width) can then be defined in terms of this axis:

  let bin_lo = x[gx=Xslp] - xbox[gx=Xslp]/2
  let bin_hi = x[gx=Xslp] + xbox[gx=Xslp]/2

And we can compute the XY average of sst in a particular bin (say 1001 mb) as 
follows:

  let binned_sst = if slp gt bin_lo[x=1001] and slp le bin_hi[x=1001] then sst
  let binned_sst_ave = binned_sst[x=@ave,y=@ave]
  
To loop over Xslp and calculate the average in _all_ bins, we just need to put 
the above 4 command lines in a script (one_bin.jnl), with the following change 
to line #3:

  let binned_sst = if slp gt bin_lo[x=$1] and slp le bin_hi[x=$1] then sst
  
Now call the script in a repeat loop:

  repeat/i=1:26 go one_bin `x[gx=Xslp]`
  
The grave accents allow the specific value of x at each i of the bin axis (e.g., 
991, 993, etc.) to be passed as a scalar to the script.  For output purposes, 
you can have the script list binned_sst_ave to the screen or to a file.  For 
example:

list/ binned_sst_ave[x=130E:160E,y=5N:45N,l=1] or
list/nohead/file=my_output.out/append/format=(f8.2)   
                                    binned_sst_ave[x=130E:160E,y=5N:45N,l=1]
                                    
To save the output, including information about the bin axis, to a NetCDF file, 
requires the following steps:

  ! save result in a temporary file with a time axis but no "bin" axis
  save/clobber/file=temp.cdf/x=130E:160E/y=5N:45N binned_sst_ave

  ! read back the temporary variable and "stamp" it with the correct bin index
  cancel variables/all
  use temp.cdf
  let bin_stamp0 = 0*x[gx=Xslp]
  let result_sst = binned_sst_ave + bin_stamp0
  
  ! save in the final file
  set variable/title="My binned SST"/units="deg. C" result_sst
  save/append/file=my_output.cdf/ilimits=1:26 result_sst[x=$1]


I hope this will be useful to someone else like it was to me.                            
  
Regards,

Chris Weaver



Department of Environmental Sciences
Rutgers University/Cook College
14 College Farm Road
New Brunswick, NJ 08901-8551
USA
(ph) 1 732 932 7902  (int) weaver@gaia.rutgers.edu