[gradsusr] bug on boundary limits of ave() function

[Henrique Barbosa]

Dear all,

I am using '-b' parameter in ave() function to calculate the average
at the exact boundary limits. However I am getting strange results
when both start and end boundaries lie inside the same bin. In that
case, ave() is returning the value at the centre of the gridbox... I
did some tests and got to the conclusion that the algorithm is
assuming a constant value in the whole grid-box. Wouldn't it be better
to do a linear interpolation instead?

In fact, if boundaries limits are inside different griboxes both
approaches will give exactly the same result... But if boundaries lie
inside the same gridbox the linear interpolation will give a more
realistic result. More over, if ave() uses a linear interpolation it
would be possible to draw a field at any world coordinate. For
instance:

'set y 1'
'd ave(prec, lat=0, lat=0)'
'd ave(prec, lat=0.1, lat=0.1)'

would draw two different curves, the first interpolated to
lat=0degrees and the second to lat=0.1degrees. Note that this would be
different from:

'd prec(lat=0)'
'd prec(lat=0.1)'

because in this second case grads tries to convert world to grid
coordinates and gives an error message if the requested latitudes
convert to non-integer  grid coordinates. Below is the output of 'q
config' and a script that can reproduce this error using the
SampleDataset distributed with OpenGrads.

[]'s
Henrique

-------------------------------------------------------------------------------------------
reinit
rc=gsfallow('on')

* open sample data file
'open /usr/local/grads-2.0.a7.oga.3/Contents/Resources/SampleDatasets/model.ctl'

* number of bins
'q dims'
line=sublin(result,2); nx=subwrd(line,13); say 'nx='nx' points'
line=sublin(result,3); ny=subwrd(line,13); say 'ny='ny' points'

* grid spacing
'set y 1 2'
l1=subwrd(result,4)
l2=subwrd(result,5)
dy=l2-l1; say 'dy='dy' degrees'

* fix latitude
'set y 1'
'set lon 150 250'

* 2m-Temperature averaged between different latitudes

ll=5.0
while(ll>0)
  say 'Click to draw average lat=[0, 'll']... check graph window to
see if it is different from previous result...'
  'q pos'
  'd ave(ts, lat=0, lat='ll',-b)'
  ll=ll-0.2
endwhile

*end
--------------------------------------------------------------------------------------------
ga->q config
Config: v2.0.a7.oga.3 little-endian readline printim grib2 netcdf
hdf4-sds hdf5 opendap-grids,stn athena geotiff
Grid Analysis and Display System (GrADS) Version 2.0.a7.oga.3
Copyright (c) 1988-2009 by Brian Doty and the
Institute for Global Environment and Society (IGES)
This program is distributed WITHOUT ANY WARRANTY
See file COPYRIGHT for more information.

Built Thu Oct 29 17:57:08 EDT 2009 for x86_64-unknown-linux-gnu

This version of GrADS has been configured with the following options:
  o Built on a LITTLE ENDIAN machine
  o Command line editing ENABLED
      http://tiswww.case.edu/php/chet/readline/rltop.html
  o printim command for image output ENABLED
      http://www.zlib.net
      http://www.libpng.org/pub/png/libpng.html
      http://www.libgd.org/Main_Page
  o GRIB2 interface ENABLED
      http://www.ijg.org
      http://www.ece.uvic.ca/~mdadams/jasper
      http://www.nco.ncep.noaa.gov/pmb/codes/GRIB2
      g2clib-1.0.5
  o NetCDF interface ENABLED
      http://www.opendap.org
      libnc-dap 4.0.1-beta3-snapshot2009021712 of Mar  3 2009 14:13:33 $
  o HDF interface ENABLED
      http://hdfgroup.org
      HDF 4.2r3
      HDF5 1.8.2
  o Athena Widget GUI ENABLED
  o OPeNDAP gridded data interface ENABLED
      http://www.opendap.org
      libdap 3.7.10
  o OPeNDAP station data interface ENABLED
      http://iges.org/grads/gadoc/supplibs.html
      libgadap 2.0.oga.1
  o GeoTIFF and KML output ENABLED
      http://www.libtiff.org
      http://geotiff.osgeo.org

For additional information please consult http://iges.org/grads