Next: , Previous: Irregular grids, Up: ncap2 netCDF Arithmetic Processor


4.1.18 Bilinear interpolation

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]. 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 ncap2 interpolation function is

     var_out=bilinear_interp(grid_in,grid_out,grid_out_x,grid_out_y,grid_in_x,grid_in_y)

where

grid_in
Input function data. Usually a two dimensional variable. It must be of size grid_in_x.size()*grid_in_y.size()
grid_out
This variable is the shape of var_out. Usually a two dimensional variable. It must be of size grid_out_x.size()*grid_out_y.size()
grid_out_x
X output values
grid_out_y
Y output values
grid_in_x
X input values values. Must be montonic (increasing or decreasing).
grid_in_y
Y input values values. Must be montonic (increasing or decreasing).
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 };
     
     // Coordinate variables
     x_in[$X]={0.0,1.0,2.0,3.01};
     y_in[$Y]={1.0,2.0,3.0,4.0,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 ;

It is possible 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 endpoints 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, though typically these are in the range (-90,90). This ncap2 script illustrates both wrapping and extrapolation of end points:

     defdim("lat_in",6);
     defdim("lon_in",5);
     
     // Coordinate input 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);
     
     // Coordinate variables
     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 ;