The lterp function performs intepolation between two grids. It is a built-in function as of version 2.0 -- in earlier versions it was implemented as a user defined function. Additional capabilities were added with version 2.0.2. The syntax is:
lterp (source, dest) (Version 2.0.1 or earlier)
lterp (source, dest <,method> <,pct>) (Version 2.0.2 or later)
Where
source a grid expression that contains the data values to be interpolated
dest a grid expression that describes the grid that the data are to be interpolated to.
The data values in dest are ignored, only the grid information is used.
method an optional keyword that describes the interpolation method to be used. Options are:
bilin bilinear interpolation -- this is the default, the behavior when no method is specified, and is identical to the pre-2.0.2 behavior
bessel an enhancement to bilinear interpolation that uses 3rd order bessel interpolation
aave area average with latitude weighting -- may be better than bilin when interpolating from higher to lower grid resolution
amean area average without latitude weighting -- same as aave, but when weighting the grid box area by latitude is not desired
pct the minimum percentage of area of each destination grid box that must contain non-missing input data to be considered valid
(pct is optional and only relevant when used with aave and amean and must be between 0 and 100; default is 50%)
The lterp function works only with gridded data. The returned result is a grid expression on the same grid as dest.
The source and dest expressions must have the same varying dimensions, which may be X, Y, or T. Interpolation is not performed in the Z or E dimension.
The source and dest expressions may vary in 1 or 2 dimensions, unless you are using the aave or amean methods, in which case the grids must be 2-D with X and Y as the varying dimensiions.
If the domain of
source is larger than the domain of dest, the returned result will have an expanded grid to cover the requested domain.
For interpolation in the time dimension, you may interpolate (A) between monthly and yearly time axes, or (B) between minute, hourly, and daily time axes.
For the bilin method, each of the grid points in the dest grid contains the average of the nearest four grid points in the source grid, weighted by their distance from the destination grid point. If any of the four surrounding input grid points contain missing data, the interpolated value will be flagged as missing.
(Version 2.0.2 or later) The bessel method is adapted from the regrid user defined function. It uses a bessel interpolation routine developed at Fleet Numerical
Meteorology and Oceanography Center (courtesy of D. Hensen).
The bilin method of interpolation is performed at all points to
produce a "first guess." Improvements to the "first
guess" are made using the higher-order terms in the bessel
interpolator. Bessel interpolation uses not just the nearst 4 grid boxes in the source grid, but also the secondary ring of grid points that form a 4x4 matrix around the destination grid point. If any of the grid boxes in the outer ring of the 4x4 matrix
contain missing data, the output grid is assigned the bilin value. The bessel method may produce a closer fit to the source data
in regions where there are large gradients (e.g., around low pressure centers).
(Version 2.0.2 or later) For the aave and amean methods, a spatial average using source data is calculated in the area outlined by each grid box in the dest grid.
Note that in GrADS the boundaries of all grid boxes are the midpoints between their centers. If a grid box in the source grid partially overlaps the area of a dest grid box, its contribution to the spatial average is weighted by the fraction of area within the dest grid box domain. When the source grid contains missing data, the pct argument may be used to specify a threshold for the percentage of the destination grid box area that must be covered with non-missing data in order to avoid being flagged as missing. The default pct value is 50%. The aave and amean methods may be used when interpolating from a higher to lower resolution grid in order to preserve statistical properties of the grid such as a gloabal mean -- e.g. the result of the expression "aave(source,global)" will be the same as "aave(lterp(source,dest,aave),global)" as long as source does not contain any missing data. The only difference between aave and amean methods is that with aave the grid box area calculations are weighted by the difference in the sine of the latitude at the northern and southern boundaries of the grid box. Use amean when the Y axis in the source grid does not represent Latitude. The aave and amean methods are based on the "box average" feature of the regrid user defined function.
open fine.ctl
open coarse.ctl
d lterp(fine,coarse.2,aave)
d lterp(fine,coarse.2,aave,100)
d lterp(coarse.2,fine,bessel)
open hourly_station_data.ctl
open 3hourly_grid_data.ctl
set x 1
set y 1
set time 00z1jan 00z1feb
d lterp(s2g1d(var.1(stid=kbwi)),var.2(lon=-77,lat=39))
'open obs.ctl'
'q file'
line5 = sublin(result,5)
sx = subwrd(line5,3)
'set x 1 'sx
'open model.ctl'
'set gxout fwrite'
'set fwrite model2obs.dat'
'd lterp(var.2,var.1)'
'disable fwrite' 'open coarse.ctl' 'set x 1 72' 'open fine.ctl' 'define newvar = lterp(ts.2(t=1),ps.1(t=1),aave)' 'set sdfwrite newvar.nc' 'sdfwrite newvar'
'open fine.ctl' 'open coarse.ctl' 'set dfile 2' 'set x 1 72' 'define newvar = lterp(ts.1(t=1),ps.2(t=1),aave)' 'set sdfwrite newvar.nc' 'sdfwrite newvar'