As of version 4.0.0 NCO has internal routines to perform bilinear interpolation on gridded data sets.
In mathematics, bilinear interpolation is an extension of linear interpolation for interpolating functions of two variables on a regular grid. The idea is to perform linear interpolation first in one direction, and then again in the other direction.
Suppose we have an irregular grid of data temperature[lat,lon]
, with co-ordinate vars lat[lat], lon[lon]
. And we wish to find the temperature at an arbitary point [X,Y] within the grid. If we can locate lat_min,lat_max and lon_min,lon_max such that lat_min <= X <= lat_max
and lon_min <=Y <=lon_max
then we can interpolate in two dimensions to find the temperature at [X,Y].
the general form of the ncap interpolation function is as follows:-
var_out=bilinear_interp(grid_in, grid_out, grid_out_x, grid_out_y, grid_in_x, grid_in_y)
grid_in
grid_in_x.size()*grid_in_y.size()
grid_out
grid_out_x.size()*grid_out_y.size()
grid_out_x
grid_out_y
grid_in_x
grid_in_y
Prior to calculations all arguments are converted to type NC_DOUBLE
. After calculations var_out
is converted to the input type of grid_in
.
Suppose the first part of an ncap2 script is:
/****************************************/ defdim("X",4); defdim("Y",5); //Temperature T_in[$X,$Y]= {100, 200, 300, 400, 500, 101, 202, 303, 404, 505, 102, 204, 306, 408, 510, 103, 206, 309, 412, 515.0 }; //Co-ordinate Vars x_in[$X]={ 0.0,1.0,2.0,3.01 }; y_in[$Y]={ 1.0,2.0,3,4,5 }; /***************************************/
Now we interpolate with the following variables:
/***************************************/ defdim("Xn",3); defdim("Yn",4); T_out[$Xn,$Yn]=0.0; x_out[$Xn]={0.0,0.02,3.01 }; y_out[$Yn]={1.1,2.0,3,4 }; var_out=bilinear_interp(T_in,T_out,x_out,y_out,x_in,y_in); print(var_out); // 110, 200, 300, 400, // 110.022, 200.04, 300.06, 400.08, // 113.3, 206, 309, 412 ; /***************************************/
Its possible to use the call to interpolate a single point:
/***************************************/ var_out=bilinear_interp(T_in,0.0,3.0,4.99,x_in,y_in); print(var_out); // 513.920594059406 /***************************************/
Wrapping and Extrapolation
The function bilinear_interp_wrap()
takes the same arguments as bilinear_interp()
but performs wrapping (Y) and extrapolation (X) for points off the edge of the grid. If the given range of longitude is say (25-335) and we have a point at 20 degrees- then the end points of the range are used for the interpolation. This is what wrapping means.
For wrapping to occur Y must be longitude and must be in the range (0,360) or (-180,180). There are no restrictions on the longitude (X) values , but typically these are in the range (-90,90).
The follwing ncap script illustrates both wrapping and extrapolation of end points.
/****************************************/ defdim("lat_in",6); defdim("lon_in",5); // co-ordinate in vars lat_in[$lat_in]={ -80,-40,0,30,60.0,85.0 }; lon_in[$lon_in]={ 30, 110, 190, 270, 350.0 }; T_in[$lat_in,$lon_in]= { 10,40,50,30,15, 12,43,52,31,16, 14,46,54,32,17, 16,49,56,33,18, 18,52,58,34,19, 20,55,60,35,20.0 }; defdim("lat_out",4); defdim("lon_out",3); // co-ordinate vars lat_out[$lat_out]={ -90, 0, 70, 88.0 }; lon_out[$lon_out]={ 0, 190, 355.0 }; T_out[$lat_out,$lon_out]=0.0; T_out=bilinear_interp_wrap(T_in,T_out,lat_out,lon_out,lat_in,lon_in); print(T_out); // 13.4375, 49.5, 14.09375, // 16.25, 54, 16.625, // 19.25, 58.8, 19.325, // 20.15, 60.24, 20.135 ; /****************************************/