2.8.14 interp2

Menu Information

Interpolate Z from XY

Brief Information

Perform interpolation on XYZ data at a given set of XY values

Additional Information

Minimum Origin Version Required: 9.1 SR0

Command Line Usage

1.interp2 ixy:=[XYZRandomGaus]"XYZ Random Gaussian"!(D,E) ir:=[XYZRandomGaus]"XYZ Random Gaussian"!(A,B,C) method:=rk pts:=3 smooth:=4;

X-Function Execution Options

Please refer to the page for additional option switches when accessing the x-function from script


XY Values to Interpolate ixy



Specify the input XY range for the given XY values.
Input ir



Specify the input XYZ range for the group of XYZ data.
Method method



Specify the interpolation method to be used. 8 methods are available.

Option list:

  • nr:Nearest
  • rk:Random Kriging
  • rc:Renka Renka Cline
  • rs:Random Shepard
  • tps:Random TPS
  • sp:Spline
  • ta:Triangle
  • wa:Weight Average

More details can be found in the Algorithm section below.

Search Radius radius



Specify a search radius. Only available when the method is either Random Kriging or Weight Average.
Minimum Points pts



Specify a minimum number of points. Only available when the method is Random Kriging.
Smoothing Factor smooth



Specify a Kriging or Spline smoothing factor. Only available when the method is either Random Kriging or Spline.
Input Weight weight



Specify a vector as the weight. Only available when the method is Spline.
Weight Function Locality Factor nw



Specify a weight function locality factor. Only available when the method is Random Shepard.
Quadratic Interpolant Locality Factor nq



Specify a quadratic interpolant locality factor. Only available when the method is Random Shepard.
Result of Interpolation oz



Specify where to output the interpolation result.


This X-Function is used to interpolate a group of existing xyz data to find z value for given (x,y), using one of eight (8) methods: Nearest, Random Kriging, Random Renka Cline, Random Shepard, Random TPS, Spline, Triangle and Weight Average.

In order to interpolate z from xy, using the interp2 X-Function:

  1. In Script Window, run
    interp2 -d;
    to open the interp2 dialog box.
  2. Pick data (x,y) from worksheet as XY Values To Interpolate, and a range of xyz data as Input for interpolation.
  3. Select a Method from the drop-down menu.
  4. Assign destination columns as Result of Interpolation.


Interp2 provides estimation of a variable at an unmeasured location from observed values at surrounding locations. There are 8 methods used in X-Function Interp2, the correspondence for each method and algorithm are listed in table below.

Methods Algorithm
Nearest Nearest-neighbor interpolation
Random Kriging Ordinary Kriging interpolation
Random Renka Cline Triangle-based Renka Cline interpolation
Random Shepard Modification of Shepard’s method interpolation
Random TPS Random TPS interpolation
Spline Bicubic B-spline interpolation
Triangle Triangle-based Renka Cline interpolation
Weight Average Inverse Distance Weighting interpolation

Nearest neighbor interpolation

In the nearest neighbor algorithm, it finds the nearest point of the given to-be-interpolated X and Y values, and then uses this found point's Z value as the interpolated result. This procedure does not include other neighboring points, that will yield a piece-wise-constant interpolation. This algorithm is so simple that its resulting surface is not smooth.

(Random)Ordinary Kriging interpolation

Origin has called this function ocmath_2d_kriging_scat to perform ordinary Kriging interpolation. The basic procedure of ordinary Kriging is described as following.

Assume x_1, x_2, ..., x_n are a series of observed random points and z(x_1), z(x_2), ..., z(x_n) are the measured values of the series points respectively. At a given point, x_0, which is not measured yet, according to ordinary Kriging, the estimated value, z^*(x_0), for this unmeasured point can be derived by a weighted linear combination of the measured points, which can be written as:

z^*(x_0) = \sum_{i=1}^nw_iz(x_i)

where w_i is the weight for its corresponding point.

To get the estimated value, z^*(x_0), for the unmeasured point x_0, it needs to figure out the weight for each measured point, which has contribution to the calculation.

First of all, we need to assumed that the random function for the points is stationary, which makes sure the unbiasedness of result. That is to say, the expected value of the following expression needs to be zero.

E[z^*(x_0)-z(x_0)] = E[\sum_{i=1}^nw_iz(x_i)-z(x_0)] = \sum_{i=1}^nw_iE[z(x_i)]-E[z(x_0)] =0

Obviously, E[z(x_i)] = E[z(x_0)] = m here, where m is the mean, so we get \sum_{i=1}^nw_i = 1.

Secondly, it should satisfy the optimal variance for the error. That is to minimize the variance, which shows as the expression below.

\sigma_0^2 = Var[z^*(x_0)-z(x_0)] 
= E[((z^*(x_0)-z(x_0))-E[z^*(x_0)-z(x_0)])^2] 
= E[(z^*(x_0)-z(x_0))^2] 
= E[(\sum_{i=1}^nw_iz(x_i)-z(x_0))^2] 
= E[(z(x_0))^2] - 2E[z(x_0)\sum_{i=1}^nw_iz(x_i)] + E[(\sum_{i=1}^nw_iz(x_i))^2]
= C(x_0, x_0) - 2\sum_{i=1}^nw_iC(x_i, x_0) + \sum_{i=1}^n\sum_{j=1}^nw_iw_jC(x_i, x_j)
= minimum

To minimize the above expression, the Lagrange parameter \mu is introduced, and let:

F = \sigma_0^2 - 2\mu (\sum_{i=1}^nw_i - 1)

Calculate the partial first derivatives of the above equation with respect to w_i and \mu and get:

\left\{\begin{matrix}\frac{\partial F}{\partial w_i} = 2\sum_{j=1}^nw_jC(x_i, x_j) - 2C(x_i, x_0) - 2\mu = 0
\\ \frac{\partial F}{\partial \mu} = -2(\sum_{i=1}^nw_i-1) = 0

Then we get the simplified one:

\left\{\begin{matrix}\sum_{j=1}^nw_jC(x_i, x_j) - \mu = C(x_i, x_0)
\\ \sum_{i=1}^nw_i = 1

and \sigma_0^2 = C(x_0, x_0) - \sum_{i=1}^nw_iC(x_i, x_0) + \mu


C = \begin{bmatrix}C(x_1, x_1) & C(x_1, x_2) & \cdots & C(x_1, x_n) & 1
\\ C(x_2, x_1) & C(x_2, x_2) & \cdots & C(x_2, x_n) & 1
\\ \vdots & \vdots & \: & \vdots & \vdots
\\ C(x_n, x_1) & C(x_n, x_2) & \cdots & C(x_n, x_n) & 1
\\ 1 & 1 & \cdots & 1 & 0
W = \begin{bmatrix}w_1 \\ w_2 \\ \vdots \\ w_n \\ -\mu\end{bmatrix},
D = \begin{bmatrix}C(x_1, x_0) \\ C(x_2, x_0) \\ \vdots \\ C(x_n, x_0) \\ 1\end{bmatrix}

Then the partial first derivatives can be expanded and rewritten as matrix form, C\cdot W = D

To solve for the weights, multiply both sides of the matrix form by C^{-1}, the inverse matrix of the C matrix, so to get W = C^{-1}\cdot D.

The final step is to decide the covariance. For all the above calculations are based on the assumption that the random variables in random function model all have the same mean and variance. Such assumption leads us to the relationship between the model variogram and the model covariance as following.

\gamma (x_i, x_j) = \frac{1}{2}E[(z(x_i)-z(x_j))^2] = \frac{1}{2}E[(z(x_i))^2] + \frac{1}{2}E[(z(x_j))^2] - E[z(x_i) \cdot z(x_j)] = E[z^2]-E[z(x_i) \cdot z(x_j)] = (E[z^2] - m^2)-(E[z(x_i) \cdot z(x_j)] - m^2) = \sigma^2-C(x_i, x_j)

Put this relationship equation to the partial first derivatives, we can get another form (the matrix form is similar):

\left\{\begin{matrix}\sum_{j=1}^nw_j\gamma (x_i, x_j) + \mu =\gamma (x_i, x_0)
\\ \sum_{i=1}^nw_i = 1

In Origin, the power model is used for estimating the varigram, which is:

\gamma (h) = \left\{\begin{matrix}0 \;\;\;\;\;\;\;\;\;\;\;\;\;\; |h| = 0
\\ a + b|h|^\lambda \;\;\;\; |h| > 0

where a = 0, b = 1, and 0 < \lambda < 2. h can refer to the distance between x_i and x_j, so \gamma (h) = \gamma (x_i, x_j)

In the dialog of this X-Function, when this method is used, \lambda corresponds to the Smoothing Factor, and Search Radius (for the search radius) and Minimum Points (minimum points included for interpolation) correspond to the radius and noctMin parameters in function ocmath_2d_kriging_scat respectively.

Random Renka Cline

The algorithm of the Random Renka Cline is based on nag_2d_scat_interpolant (e01sac) function provided by NAG. However, according to the NAG Library: Advice on Replacement Calls for Withdrawn/Superseded Functions, the nag_2d_scat_interpolant (e01sac) is replaced by nag_2d_shep_interp (e01sgc) or nag_2d_triang_interp (e01sjc). Thus, the Random Renka Cline method has the algorithm similar to Triangle method which based on nag_2d_triang_interp (e01sjc), while the latter is recommended due to the smaller maximum interpolation errors could be obtained.

Random Shepard

Please refer to the (Random) Shepard for algorithm introduction.

Random TPS

Please refer to the Random TPS for algorithm introduction.


This method based on the nag_2d_spline_fit_scat (e02ddc) which computes a bicubic spline approximation to a set of scattered data.


This method involves firstly creating Delaunay triangles by connecting all data points. The triangles are as nearly equiangular as possible and no two triangles should intersect. Refer to this document for the triangulation theory.

After the triangulation been created, the triangle-based interpolation method computes a value at a point based only on data values and first partial derivatives at the three vertices of the triangle containing the point. It calls nag_2d_triang_interp (e01sjc), which uses a method due to Renka and Cline.

Weighting average

The weighted average method is based on ocmath_gridding_weighted_average. The unknown value P(x_p, y_p, z_p) can be calculated from the surrounded true values Q_i within the specified search radius r. The interpolate value can be expressed by:


d_i is the distance between P and Q_i. If there is no value within the search radius, Zp will be missing value. Increasing the search radius means that each point is more interrelated to neighboring points, so to produce a smoother surface that may lose fine details.


Here we offer the comparison of the eight methods. We apply the eight methods to six functions from F1 to F6, which corresponds to exponential peak, cliff surface, saddle surface, gentle peak, steep peak, and sphere surface, respectively. The method is derived from Ref. 2. The interpolation for 1089 points is based on 100 nodes uniform distribution. Mean error and Maximum error are average and max of the difference between real values and interpolated values for 1089 points, respectively. The bar chart for mean error and maximum error are shown below:

Analysis 3.png

Function expression:

F1: 0.75e^{-\frac{(9x-2)^2+(9y-2)^2}{4}}+0.75e^{\frac{-(9x+1)^2}{49}-\frac{9y+1}{10}}+0.5e^{-\frac{(9x-7)^2+(9y-3)^2}{4}}-0.2e^{-(9x-4)^2-(9y-7)^2}

F2: \frac{tanh(9y-9x)+1}{9}

F3: \frac{1.25+cos(5.4y)}{6+6(3x-1)^2}

F4: e^{\frac{-(81/16)((x-0.5)^2+(y-0.5)^2)}{3}}

F5: e^{\frac{-(81/4)((x-0.5)^2+(y-0.5)^2)}{3}}

F6: \frac{\sqrt{64-81((x-0.5)^2+(y-0.5)^2}}{9}-0.5


  1. Isaaks and Srivastava. An Introduction to Applied Geostatistics, Chapter 12. Oxford University Press (1989)
  2. R. J. Renka and A. K. Cline. A Triangle-Based C^1 Interpolation Method, Rocky Mountain Journal of Mathematics. 14, 223-237 (1984)
  3. P. Dierckx. An Algorithm for Surface-Fitting with Spline Functions, IMA Journal of Numerical Analysis. 1, 267-283 (1981)

Related X-Functions