Interpolated Plotting and Smoothing:
Example Worksheet
The topic of this worksheet is interpolation and smoothing of given 2dimensional and 3dimensional data. The organization is by dimension, task, and regularity of the data.
One common task is to generate only a plot of the data, as a curve or surface which passes through or approximates the data points.
Another task is to generate a procedure or piecewise spline expression which approximates the data and which can be queried for a value at any individual point nearby the original data points.
Another choice is whether to produce a curve or surface which passes through the given data points (and so interpolates them) or which simply approximates them (by smoothing). Is the data known to contain noise or measurement error? If the dependent data is expected to contain error or noise then it is reasonable to fit a smoothed surface or curve to the data, where the smoothed fit does not necessarily pass through or match the dependent data at the given data points. On the other hand, if the data is expected to be fully correct then it is reasonable to produce an interpolated curve or surface which passes through all the given data points.
A curve or surface may also be plotted directly from the data, or indirectly by first creating a procedure or (in the 2D case) a piecewise spline, and then subsequently supplying that to the plot or plot3d command. The constructed procedure would accept individual 1D or 2D independent data points and compute a dependent scalar value for that input.
Another relevant distinction for 3D data is whether the independent portion of the data lies on a regular grid in the 2D plane or whether it is comprised of points which are irregularly spaced in both directions.
This worksheet elaborates and compares these approaches for 2D and 3D plotting.
>

restart;
with(LinearAlgebra): with(CurveFitting): with(Interpolation): with(plots): with(Statistics):
randomize():


1dimensional case


The xdata points are taken uniformly in the range, as the integers from 1 to 15. But they could just as well be unequally spaced, as long as they are ordered.
Here, the data2D is a 15x2 Matrix with the first column being the ordered xdata points and the second column being their yvalues.
>

xdata := <($1..15)>:
ydata := RandomVector(15, generator = 1.0 .. 5.0):
data2D:=< xdata  ydata >:

Here is the plot command producing a linear interpolation. The goal is to produce a smooth curve instead.
>

display( plot(data2D, color=black),
pointplot(data2D, symbolsize=10, color=red) );

We produce an interpolating function using splines, and plot this function.
>

f := Interpolate(xdata, ydata):
display(plot(f, 1..15, color=black),
pointplot(data2D, symbolsize=10, color=red) );

The above method of evaluating the interpolating spline is simplest. It uses the Interpolate command from the Interpolation package. Behind the scenes, it uses the ArrayInterpolation command from the CurveFitting package, which is called in an efficient manner for all the interpolating points that the plotting procedure wants to evaluate. Using options, one can select between different types of interpolants.
>

g := Interpolate(xdata, ydata, method=cubic):
h := Interpolate(xdata, ydata, method=radialbasisfunction, thinplatespline):
display(plot([f, g, h], 1..15, color=[black, blue, green]),
pointplot(data2D, symbolsize=10, color=red) );

Another use of an interpolant generated by the Interpolate command is to do numeric integration. The following commands determine the area under the three curves in the previous plot.
>

int(f, 1 .. 15, numeric);
int(g, 1 .. 15, numeric);
int(h, 1 .. 15, numeric);

Another option is local optimization. We have chosen the data randomly so far, which is not very suitable for demonstrating this functionality. Let's use an example that is more suitable. We assume the following data come from an experiment where we want to find the (locally) maximal value. We see that the interpolated curve has a maximum close to the maximum of the underlying ground truth (which would be the maximum 1 at $xequals;\mathrm{\pi}$). This example uses the command Optimization[Maximize].
>

ground_truth := x > 1/(1 + (x  Pi)^2):
xsampled := [seq(i, i = 0. .. 5, 0.7)]:
ysampled := map(ground_truth, xsampled):
interpolated := Interpolate(xsampled, ysampled):
optimum := Optimization[Maximize](unapply(interpolated(x), x, numeric), 0 .. 5);
plots:display(plot([ground_truth, interpolated], 0 .. 5, linestyle=["dash", "solid"]),
pointplot(<xsampled, ysampled>, symbolsize=10, color=blue),
pointplot([optimum[2][1], optimum[1]], symbolsize=25, color=green, symbol=soliddiamond));

$\left[{0.9635817101770758}{\,}\left[\begin{array}{c}3.1369851865331473\end{array}\right]\right]$
 
 

Efficiency


The object returned by the Interpolate command uses the command CurveFitting[ArrayInterpolation] behind the scenes, for the methods that are listed on the Interpolate help page as accepting 1dimensional data only. There is some overhead in setting up each individual call to that command; the Interpolate object minimizes this overhead, but there is significant overhead inherent in the model where there is a call to ArrayInterpolation for each value computed. Conversely, if you can arrange to compute many values with a single call, this is more efficient. We will illustrate this below, reusing the first example set of data, and using the command CodeTools:Usage to print information about memory and time resources used.
A relatively large number of 5000 interpolation points will be used in order to more clearly demonstrate some performance differences, and each example will be run 10 times. This large number of points produces some jaggedness in the plot, but we're not primarily concerned with the appearance of the resulting plot. We create this plot in three different ways:
1.

With a single direct call to the ArrayInterpolation command, on all 5000 points at once.

2.

Using the result from Interpolate, f, to compute all 5000 points at once; this uses a calling sequence we haven't seen yet in this document, where you submit a Vector of data points to f.

3.

Letting the plot command call f, with a separate call for each point.

>

interpolationpoints := Vector(5000, i > i*15/5000, datatype=float[8]):

>

CodeTools:Usage(plot(<interpolationpoints  ArrayInterpolation(data2D, interpolationpoints, method=spline)>), iterations=10);

memory used=3.26MiB, alloc change=32.00MiB, cpu time=28.80ms, real time=29.20ms, gc time=5.60ms
 
>

CodeTools:Usage(plot(<interpolationpoints  f(interpolationpoints)>), iterations=10);

memory used=3.23MiB, alloc change=0 bytes, cpu time=22.40ms, real time=22.70ms, gc time=0ns
 
>

CodeTools:Usage(plot(f, 1 .. 15, sample=convert(interpolationpoints, list), adaptive=false), iterations=10);

memory used=89.99MiB, alloc change=8.00MiB, cpu time=869.60ms, real time=811.40ms, gc time=118.80ms
 
We see that the first two options are about equally fast, and the third is 3040 times slower. So: it is very advantageous for efficiency if one can pick a large number of points to interpolate at, at once.
In order to minimize overhead, the object created by Interpolate does the following:
 Upon construction of the object, it preprocesses the input data, by converting it to a floatingpoint rtable and sorting it. This means that this happens only once, rather than with each call.
 The object preallocates a single element rtable with the correct data type for the result. This prevents any allocation from taking place for calls that seek to compute a single value.


Obtaining a closedform expression


The above methods can be used for any application where you need to evaluate the interpolant numerically; it cannot, however, produce a closedform (typically piecewise) expression for the interpolating function. Such closedform expressions can be obtained from some commands in the CurveFitting package.
>

Spline(xdata, ydata, v);

$\left\{\begin{array}{cc}{6.40737377947708}{}{3.29095019314989}{}{v}{+}{1.38269059814989}{}{\left({v}{}{1}\right)}^{{3}}& {v}{<}{2}\\ {}{0.506079211218933}{+}{0.857121601299780}{}{v}{+}{4.14807179444967}{}{\left({v}{}{2}\right)}^{{2}}{}{2.86512326974945}{}{\left({v}{}{2}\right)}^{{3}}& {v}{<}{3}\\ {1.67454797453381}{+}{0.557895380950769}{}{v}{}{4.44729801479868}{}{\left({v}{}{3}\right)}^{{2}}{+}{1.92150000984791}{}{\left({v}{}{3}\right)}^{{3}}& {v}{<}{4}\\ {11.6691339690577}{}{2.57220061910286}{}{v}{+}{1.31720201474505}{}{\left({v}{}{4}\right)}^{{2}}{+}{0.988945535357804}{}{\left({v}{}{4}\right)}^{{3}}& {v}{<}{5}\\ {}{14.0309216587970}{+}{3.02904001646066}{}{v}{+}{4.28403862081847}{}{\left({v}{}{5}\right)}^{{2}}{}{3.81310060227913}{}{\left({v}{}{5}\right)}^{{3}}& {v}{<}{6}\\ {3.66736375129320}{+}{0.157815451260210}{}{v}{}{7.15526318601892}{}{\left({v}{}{6}\right)}^{{2}}{+}{4.56321490775871}{}{\left({v}{}{6}\right)}^{{3}}& {v}{<}{7}\\ {5.42148701435834}{}{0.463066197501500}{}{v}{+}{6.53438153725721}{}{\left({v}{}{7}\right)}^{{2}}{}{3.98741633475571}{}{\left({v}{}{7}\right)}^{{3}}& {v}{<}{8}\\ {}{0.883660344840631}{+}{0.643447872745792}{}{v}{}{5.42786746700992}{}{\left({v}{}{8}\right)}^{{2}}{+}{2.95974476326412}{}{\left({v}{}{8}\right)}^{{3}}& {v}{<}{9}\\ {14.4367227496527}{}{1.33305277148167}{}{v}{+}{3.45136682278246}{}{\left({v}{}{9}\right)}^{{2}}{}{1.79965285430079}{}{\left({v}{}{9}\right)}^{{3}}& {v}{<}{10}\\ {1.05068589165074}{+}{0.170722311180876}{}{v}{}{1.94759174011991}{}{\left({v}{}{10}\right)}^{{2}}{+}{1.21867068593904}{}{\left({v}{}{10}\right)}^{{3}}& {v}{<}{11}\\ {2.95265048389397}{}{0.0684491112418367}{}{v}{+}{1.70842031769720}{}{\left({v}{}{11}\right)}^{{2}}{}{0.717063285455365}{}{\left({v}{}{11}\right)}^{{3}}& {v}{<}{12}\\ {}{11.2438018327704}{+}{1.19720166778647}{}{v}{}{0.442769538668893}{}{\left({v}{}{12}\right)}^{{2}}{+}{0.313852219882421}{}{\left({v}{}{12}\right)}^{{3}}& {v}{<}{13}\\ {}{12.1009477207848}{+}{1.25321925009595}{}{v}{+}{0.498787120978371}{}{\left({v}{}{13}\right)}^{{2}}{}{1.00162767607432}{}{\left({v}{}{13}\right)}^{{3}}& {v}{<}{14}\\ {15.4985347316823}{}{0.754089536170272}{}{v}{}{2.50609590724459}{}{\left({v}{}{14}\right)}^{{2}}{+}{0.835365302414864}{}{\left({v}{}{14}\right)}^{{3}}& {\mathrm{otherwise}}\end{array}\right.$
 (1.2.1) 



2 and higherdimensional case


Below we'll distinguish between data interpolation, which matches the given data points exactly, and data smoothing, which approximates noisy data.

Smoothing


In this section, we are assuming that we have a collection of 2 or higherdimensional independent data for which the 1dimensional dependent data has some error component or is known to be noisy. The presence of noise implies that whatever surface this data represents does not necessarily pass through the given dependent data points.
In this case, a surface is approximated by numerically smoothing the data using the lowess algorithm. For any given point to be plotted a window of close enough, surrounding data points will be used to compute a local, weighted, low order fit.

Scatterplot3D


The command ScatterPlot3D from the Statistics package provides a smoothed plot of the 2D noisy data. The independent data does not need to be regularly spaced, and is supplied as an nby3 Array or Matrix. Each of the n rows represents an individual point. The columns are the x, y, and zvalues.
Here is an example. First, we'll look down upon the projection of the data points onto the xy plane, and thus visualize the layout and spacing of the 2D dependent data values.
>

X := Sample(Uniform(50,50),175):

>

Y := Sample(Uniform(50,50),175):

>

Zerror := Sample(Normal(0,100),175):

>

Z := Array(1..175,(i)>(sin(Y[i]/20)*(X[i]6)^2+(Y[i]7)^2+Zerror[i])):

>

XYZ := Matrix([[X],[Y],[Z]],datatype=float[8])^%T;

>

ScatterPlot3D(XYZ, axes=box, orientation=[20,0,0]);

>

ScatterPlot3D(XYZ, lowess, grid=[25,25], axes=box, orientation=[20,70,0]);

Here is another example, which reads in the nby3 data from a file.
>

M:=ImportMatrix(cat(kernelopts(mapledir),"/data/plotting/irregulardata3D.csv"),
source=csv,datatype=float[8]):

>

Ppt := pointplot3d(M,axes=box,style=patchnogrid,symbolsize=10,color=red,
labels=[x,y,z]):

>

display(Ppt,orientation=[90,0,0]);

The default behavior for the ScatterPlot3D command is to use an order 2 quadratic form, and to use a window of width 1/3 of the original domain in each of the x and ydirections.
On the right below is the smoothed surface, while on the left is a slightly transparent view of the same surface overlaid with a pointplot of the original data points.
>

Ploess := ScatterPlot3D( M, lowess ):

>

display(Array([
display([Ppt,Ploess], transparency=0.2),
display(Ploess,style=patch,color=RGB(0.0,0.4,0.6))
]),view=700..860,axes=box,orientation=[54,78,15],
transparency=0.0);



Lowess


If you need to do more than just plot the surface, you can use the CurveFitting[Lowess] command. It returns a procedure that runs the lowess algorithm. This procedure can deal with arbitrarydimensional input. Of course, if the input points are 2dimensional, one can use this procedure for a plot, as well.
${841.914000008457}$
 (2.1.2.1) 
>

optimum := Optimization:Maximize(p(x, y), x=6000 .. 6000, y=6000 .. 6000);

${\mathrm{optimum}}{\u2254}\left[{849.801754007939735}{\,}\left[{x}{=}{1865.27968390193}{\,}{y}{=}{\mathrm{1327.90514681541}}\right]\right]$
 (2.1.2.2) 
>

plots:display(
plot3d(p, 6000 .. 6000, 6000 .. 6000),
pointplot3d([op(eval([x, y], optimum[2])), optimum[1]], symbolsize=25,
color=red, orientation=[54,78,15]));

Below, we plot the (numerically estimated) derivative of the smoothed surface in the x and ydirection, in red and blue, respectively.
>

plot3d([fdiff(p, [1], [x, y]), fdiff(p, [2], [x, y])], x = 6000 .. 6000, y = 6000 .. 6000,
color=[red, blue]);




Interpolation


In this section we are assuming that we have a collection of 2 or higherdimensional independent data for which the 1dimensional dependent data is accepted as correct. The goal is to produce a surface which must pass through all of the data points. That is, at every xy point of the independent data the height of the plotted surface must match the corresponding value of the dependent (z) data.
In this case, a surface is computed numerically by interpolating the data.
This falls into one of two distinct situations. The first case has all the independent xy data forming a regular grid, and the second case consists of irregularly spaced independent data points in the xy plane. In both cases, the surface may be plotted using the surfdata command.

Data points in a grid


In the case of a uniform 2dimensional grid of data a surface plot can be generated using the surfdata command.
The data is comprised of a grid of x and ypoints. For the examples of this section, the data points are taken uniformly in each direction. The xvalues are taken as the integers from 1 to 7, and the yvalues are taken as the integers from 1 to 9. (Nonuniform grid data can be processed in the same way, but for nongrid data, see the later sections.)

View of data points in the xy plane


Here is a view of the input data points, which in this example is a full grid, uniformly spaced in both the x and y directions.
>

xpts:=<($1..7)>:
ypts:=<($1..9)>:

>

display( seq(plottools:line([xpts[i],1],[xpts[i],9]), i=1..7),
seq(plottools:line([1,ypts[j]],[7,ypts[j]]), j=1..9),
pointplot([seq(seq([xpts[i],ypts[j]],j=1..9),i=1..7)], color=red) );


The given data values corresponding to each data point are stored in data which is a 7x9 Array. A random collection of values is generated for this example.
>

data3D:=Array( LinearAlgebra:RandomMatrix(7, 9, generator=0.2 .. 0.8), datatype=float[8]):

Computations for such data can be done using some commands in the Interpolation package: the commands listed on the help page for SplineInterpolation. These commands call CurveFitting[ArrayInterpolation] behind the scenes. They accept arbitrarydimensional grid data; in order to be able to show plots, we use 2dimensional input points.
>

f := SplineInterpolation([xpts, ypts], data3D);

${\mathrm{Interpolation}}{:}{\mathrm{SplineInterpolation}}{}\left(\left[\left[\begin{array}{c}1.0\\ 2.0\\ 3.0\\ 4.0\\ 5.0\\ 6.0\\ 7.0\end{array}\right]{\,}\left[\begin{array}{c}1.0\\ 2.0\\ 3.0\\ 4.0\\ 5.0\\ 6.0\\ 7.0\\ 8.0\\ 9.0\end{array}\right]\right]{\,}\left[\begin{array}{ccccccccc}0.4394237129718905& 0.2334864716054194& 0.3530279668184577& 0.5811382828416922& 0.329132906753384& 0.4257115852204389& 0.7407807077973092& 0.6045309713542488& 0.5830824563358177\\ 0.6537581716442946& 0.3283500920129456& 0.5803027267434254& 0.49909705039389723& 0.5344330142446085& 0.36088329441371786& 0.4782985287500852& 0.739674339722364& 0.4606682853258782\\ 0.44754345383522176& 0.37997366187623427& 0.6282167186778296& 0.4567297883426742& 0.39488827113152647& 0.38565812448548487& 0.40260018867213376& 0.43445161412115535& 0.493520873105826\\ 0.26813123992450677& 0.7684096138696237& 0.2725921057892862& 0.5356245165007512& 0.2854508694100452& 0.3507484652973021& 0.5119334605079597& 0.2318917785390107& 0.7145829159779697\\ 0.2622267959819816& 0.740189323829678& 0.4680381928798836& 0.2872738022939716& 0.6146244293662912& 0.5322672878666747& 0.7086220044857923& 0.22497458657142758& 0.4467450022571856\\ 0.6487037917023593& 0.3293460551570467& 0.28331834999781264& 0.42989787998858314& 0.5682950163278608& 0.6582734443248712& 0.24373262362755188& 0.4052123669097832& 0.511058941604045\\ 0.3682588366844777& 0.5216256708793106& 0.6899911099144254& 0.4393607107008185& 0.26521738844863446& 0.5812471055229301& 0.39596458852053373& 0.6142710156337192& 0.775137710060507\end{array}\right]{\,}{\mathrm{verify}}{\=}{\mathrm{false}}\right)$
 (2.2.1.1) 
>

g := LinearInterpolation([xpts, ypts], data3D);

${\mathrm{Interpolation}}{:}{\mathrm{LinearInterpolation}}{}\left(\left[\left[\begin{array}{c}1.0\\ 2.0\\ 3.0\\ 4.0\\ 5.0\\ 6.0\\ 7.0\end{array}\right]{\,}\left[\begin{array}{c}1.0\\ 2.0\\ 3.0\\ 4.0\\ 5.0\\ 6.0\\ 7.0\\ 8.0\\ 9.0\end{array}\right]\right]{\,}\left[\begin{array}{ccccccccc}0.4394237129718905& 0.2334864716054194& 0.3530279668184577& 0.5811382828416922& 0.329132906753384& 0.4257115852204389& 0.7407807077973092& 0.6045309713542488& 0.5830824563358177\\ 0.6537581716442946& 0.3283500920129456& 0.5803027267434254& 0.49909705039389723& 0.5344330142446085& 0.36088329441371786& 0.4782985287500852& 0.739674339722364& 0.4606682853258782\\ 0.44754345383522176& 0.37997366187623427& 0.6282167186778296& 0.4567297883426742& 0.39488827113152647& 0.38565812448548487& 0.40260018867213376& 0.43445161412115535& 0.493520873105826\\ 0.26813123992450677& 0.7684096138696237& 0.2725921057892862& 0.5356245165007512& 0.2854508694100452& 0.3507484652973021& 0.5119334605079597& 0.2318917785390107& 0.7145829159779697\\ 0.2622267959819816& 0.740189323829678& 0.4680381928798836& 0.2872738022939716& 0.6146244293662912& 0.5322672878666747& 0.7086220044857923& 0.22497458657142758& 0.4467450022571856\\ 0.6487037917023593& 0.3293460551570467& 0.28331834999781264& 0.42989787998858314& 0.5682950163278608& 0.6582734443248712& 0.24373262362755188& 0.4052123669097832& 0.511058941604045\\ 0.3682588366844777& 0.5216256708793106& 0.6899911099144254& 0.4393607107008185& 0.26521738844863446& 0.5812471055229301& 0.39596458852053373& 0.6142710156337192& 0.775137710060507\end{array}\right]{\,}{\mathrm{verify}}{\=}{\mathrm{false}}\right)$
 (2.2.1.2) 
>

plot3d([f(x,y), g(x,y) + 0.2], x=1..7, y=1..9, color=[red, blue]);

>

Optimization:Minimize(f(x, y), x = 1..7, y=1..9);

$\left[{0.236597474643406969}{\,}\left[{x}{=}{3.83202160187721}{\,}{y}{=}{5.29959188777011}\right]\right]$
 (2.2.1.3) 
>

fsolve({(a,b)>ab, (a,b)>f(a,b)0.5}, [1..7, 1..9]);

$\left[{2.4966859680690522}{\,}{2.4966859680690522}\right]$
 (2.2.1.4) 
${\mathrm{7.10542735760100}}{}{{10}}^{{\mathrm{15}}}$
 (2.2.1.5) 
The plot above uses the normal plot3d command, which does not know to select its sample points at the exact locations where our data points are, because it does not receive the actual data, only the interpolating object. The surfdata command does know to do this, because it is given the actual data.
The first 3D plot below is the piecewiseplanar surface produced by default. The goal is to produce a smoother surface instead, and this will be accomplished by using the gridsize and interpolation options of the surfdata command.
The plotted surface is being overlaid here with both the 3D pointplot as well as the (linear) patchsurface, so that the computed surface may be visually demonstrated to be authentic with respect to the original discrete data.
>

ptsplot:=PLOT3D( GRID(1..7, 1..9, data3D, STYLE(POINT) ) ):
gridplot:=PLOT3D( GRID( 1..7, 1..9, data3D, STYLE(WIREFRAME)) ):

The following plot3d options control the general look of these plots, and are reused in this subsection. They are assigned to a single, reusable name so that the differences between methods is more clearly illustrated.
>

lookandfeel := axes=boxed, style=surface, glossiness=1.0, lightmodel=light1,
color=RGB(204/255, 84/255, 0/255), view=[1..7, 1..9, 0..1],
orientation = [60, 51, 6]:

>

display(
ptsplot,gridplot,
surfdata(data3D, 1..7, 1..9, lookandfeel)
);

On the left below is the default surface produced by the surfdata command, without any additional interpolation.
On the right below is the surface produced by the surfdata command using the new interpolating optional keyword parameters gridsize and interpolation.
The gridsize parameter is an equation of the form gridsize=list, where the the list contains a pair of positive integers denoting the size of the uniform grid of points at which interpolation of the original data will be performed.
The interpolation parameter is an equation of the form interpolation=list, where the list contains options accepted by the CurveFitting:ArrayInterpolation command.
>

display( Array([
surfdata( data3D, 1..7, 1..9, lookandfeel ),
surfdata( data3D, 1..7, 1..9,
gridsize=[28,36], interpolation=[method=spline],
lookandfeel)
]) );

>

display( Array([
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9, lookandfeel ) ),
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9,
gridsize=[28,36], interpolation=[method=spline],
lookandfeel) )
]) );

Note that the 3D plot renderer does its own small amount smoothing of the surface. Hence, even when using the purely linear method of the computational interpolation scheme, the plot on the right below shows a modest level of surface smoothing.
>

display( Array([
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9, lookandfeel ) ),
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9,
gridsize=[25,33], interpolation=[method=linear],
lookandfeel) )
]) );

The refined grid on which values are interpolated does not need to contain a whole integer multiple (in either direction) of the number points of the original grid of data points.
>

display( Array([
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9, lookandfeel ) ),
display( ptsplot, gridplot,
surfdata( data3D, 1..7, 1..9,
gridsize=[33,33], interpolation=[method=cubic],
lookandfeel) )
]) );



Nonuniform data grid


This is the case of irregularly spaced data. In the case of twodimensional independent data and one dependent dimension, this is typically supplied as an nby3 Array or Matrix. Each row represents a point, which the columns interpreted as x, y, and zvalue.
Here below is an example. First, we'll look down upon the projection of the data points onto the xy plane and thus visualize the layout and spacing of the 2D dependent data values.
This is the same data that was used for the ScatterPlot3D example in the Smoothing section. But in contrast to the smoothing example we are now supposing that the data points are to be interpolated, which is to say that the surface must pass directly through the data points.
>

M:=ImportMatrix(cat(kernelopts(mapledir),"/data/plotting/irregulardata3D.csv"),
source=csv,datatype=float[8]):
points := M[.., 1..2]:
values := M[.., 3]:

>

Ppt := pointplot3d(M,axes=box,style=patchnogrid,symbolsize=15,color=red,
labels=[x,y,z]):

>

display(Ppt,orientation=[90,0,0]);

There are many different approaches to working with such data; we will examine the five that are implemented in Maple. The first two methods, natural neighbor interpolation and linear triangular interpolation, determine a socalled Delaunay triangulation of the convex hull of the input points, and use this triangulation to determine which points contribute to the interpolated value. For visualizing such a Delaunay triangulation, one can use the surfdata command.
The option `source=irregular` is supplied to denote that this data is not to be interpreted as forming a regularly spaced grid with only three values along one of the two independent data dimensions.
>

sd := surfdata(M,source=irregular,axes=box, thickness=0):
plots:display(Ppt, sd);

If you want to do computations with the interpolant, rather than just generating a plot, you can use the linear triangular method of interpolation. The interpolant is exactly the surface you see in the plot. You can obtain this interpolant from the Interpolation[Interpolate] command by giving it the option $\mathrm{method}equals;\mathrm{lineartriangular}$, or equivalently by calling the LinearTriangularInterpolation command in the same package. The interpolant is undefined at points outside the convex hull of the input points.
>

f_lt := Interpolate(points, values, method = lineartriangular);

${\mathrm{f\_lt}}{\u2254}\left({\mathrm{Linear\; Triangular\; interpolation\; object\; with\; 176\; sample\; points}}\right)$
 (2.2.2.1) 
${790.016559118906}$
 (2.2.2.2) 
Natural neighbor interpolation is the default method employed by Interpolation[Interpolate] for twodimensional input points (so the problem is threedimensional if you include the values). You can get the same object by calling the NaturalNeighborInterpolation command in the same package. This interpolation method is based on the Delaunay triangulation shown above, but leads to a smoother interpolant than the linear triangular method. As a consequence of its dependence on the Delaunay triangulation, the interpolated value is undefined if the sample point is not within the convex hull of the input points. This, combined with the sampling grid of the $\mathrm{plot3d}$ command, causes the jagged edges in the plot below.
>

f_nn := Interpolate(points, values);

${\mathrm{f\_nn}}{\u2254}\left({\mathrm{Natural\; Neighbor\; interpolation\; object\; with\; 176\; sample\; points}}\right)$
 (2.2.2.3) 
>

f_nn(4000, 2000);
plots:display(Ppt, plot3d(f_nn, 6100 .. 6000, 6600 .. 6600, grid=[149,149]));

The previous two implementation methods are implemented in Maple only for twodimensional input points. The radial basis function, inverse distance weighted, and Kriging interpolation methods are implemented for any dimension  but here we will show them with twodimensional input points.
The radial basis function interpolant is a linear combination of radial basis functions centered at each input point, with the coefficients chosen so that the result is an interpolant. You can use this interpolation method by calling the Interpolation[Interpolate] command with the option $\mathrm{method}equals;\mathrm{radialbasisfunction}$, or equivalently by calling the RadialBasisFunctionInterpolation command in the same package. Several different radial basis functions are available; they are listed on the RadialBasisFunctionInterpolation help page and they can be selected by passing the name of the radial basis function to either of these commands (in some cases with an optional numeric shape parameter).
>

f_rb := Interpolate(points, values, method = radialbasisfunction);

${\mathrm{f\_rb}}{\u2254}\left(\begin{array}{c}{\mathrm{Ra\ⅆial\; Basis\; Function\; int\ⅇrpolation\; ob\u0237\ⅇct\; with\; 176\; sampl\ⅇ\; points}}\\ {\mathrm{Ra\ⅆial\; Basis\; Function:\; multiqua\ⅆric}}\end{array}\right)$
 (2.2.2.4) 
>

f_rb_g := Interpolate(points, values, method = radialbasisfunction, inversequadratic, .0008);

${\mathrm{f\_rb\_g}}{\u2254}\left(\begin{array}{c}{\mathrm{Ra\ⅆial\; Basis\; Function\; int\ⅇrpolation\; ob\u0237\ⅇct\; with\; 176\; sampl\ⅇ\; points}}\\ {\mathrm{Ra\ⅆial\; Basis\; Function:\; inv\ⅇrs\ⅇqua\ⅆratic}}\end{array}\right)$
 (2.2.2.5) 
>

f_rb(4000, 2000);
f_rb_g(4000, 2000);
plots:display(Array(1..2, map(f >
plots:display(Ppt,
plot3d(f, 6100 .. 6000, 6600 .. 6600),
view = [DEFAULT, DEFAULT, 700..900]), [f_rb, f_rb_g])));

The inverse distance weighted interpolant is obtained by weighting the values associated with all points with the distance to the query point. An extra parameter can be supplied, which is the distance beyond which input points should be disregarded completely. You can use this interpolation method by calling the Interpolation[Interpolate] command with the option $\mathrm{method}equals;\mathrm{inversedistanceweighted}$, or equivalently by calling the InverseDistanceWeightedInterpolation command in the same package
>

f_id := Interpolate(points, values, method=inversedistanceweighted);

${\mathrm{f\_id}}{\u2254}\left(\begin{array}{c}{\mathrm{Inv\ⅇrs\ⅇ\; Distanc\ⅇ\; W\ⅇight\ⅇ\ⅆ\; int\ⅇrpolation\; ob\u0237\ⅇct\; with\; 176\; sampl\ⅇ\; points}}\\ {\mathrm{Ra\ⅆius\; of\; influ\ⅇnc\ⅇ:\; infinity}}\end{array}\right)$
 (2.2.2.6) 
>

f_id_2000 := Interpolate(points, values, method=inversedistanceweighted, 2000);

${\mathrm{f\_id\_2000}}{\u2254}\left(\begin{array}{c}{\mathrm{Inv\ⅇrs\ⅇ\; Distanc\ⅇ\; W\ⅇight\ⅇ\ⅆ\; int\ⅇrpolation\; ob\u0237\ⅇct\; with\; 176\; sampl\ⅇ\; points}}\\ {\mathrm{Ra\ⅆius\; of\; influ\ⅇnc\ⅇ:\; 2000.}}\end{array}\right)$
 (2.2.2.7) 
>

f_id(4000, 2000);
f_id_2000(4000, 2000);
plots:display(Array(1..2, map(f >
plots:display(Ppt,
plot3d(f, 6100 .. 6000, 6600 .. 6600, grid=[149, 149])), [f_id, f_id_2000])));

Finally, Kriging is an interpolation method based on a statistical model of how distance between two points influences the correlation between the values at those two points. This model is called a variogram. Several parameterized variograms are implemented in Maple, as listed on the Interpolation[Kriging][SetVariogram] page. Variogram parameters can be inferred automatically, but the best results are obtained when the variogram is chosen manually. Sometimes theoretical considerations about the problem, for example from literature research, can give reasons to choose a particular variogram; at other times, you can estimate a variogram from the input points and values. That is what we do in this case. We start by examining the empirical variogram, that is, the graph of distances between pairs of points and absolute differences in value between them. This can be done with the Interpolation[Kriging][DisplayVariogram] command. We can call this command once we have created a Kriging object, using the command Interpolation[Interpolate] with the option $\mathrm{method}equals;\mathrm{kriging}$, or the equivalent command Interpolation[Kriging].
>

f_kr := Interpolate(points, values, method=kriging);

${\mathrm{f\_kr}}{\u2254}\left(\begin{array}{c}{\mathrm{Kriging\; int\ⅇrpolation\; ob\u0237\ⅇct\; with\; 176\; sampl\ⅇ\; points}}\\ {\mathrm{Variogram:\; Sph\ⅇrical(81.,1444.,5634.358496)}}\end{array}\right)$
 (2.2.2.8) 
>

DisplayVariogram(f_kr, output=groupedpoints);

This shows that the values are generally closer when points are closer, but it is hard to see what the relationship might look like. We tweak the display a bit by grouping distancedifference pairs together and averaging them. That has happened already for the graph above, but we need to do so more aggressively.
>

binned := DisplayVariogram(f_kr, output=bins, bounds=25):
binmeans := map(convert, map(Statistics:Mean, binned), list):
binmeans_t := ListTools:Transpose(binmeans):

>

pointplot(binmeans, view=[0 .. max(binmeans_t[1]), 0 .. max(binmeans_t[2])]);

This suggests that the minimum expected difference is about 100, the maximum is about 1600, and that maximum is reached when the distance between the points is around 8000. This would be reasonably well represented by the variogram $\mathrm{Circular}\left(100\,1600comma;8000\right)$.
>

SetVariogram(f_kr, Circular(100, 1600, 8000));

$\left(\begin{array}{c}{\mathrm{Kriging\; int\ⅇrpolation\; ob\u0237\ⅇct\; with\; 176\; sampl\ⅇ\; points}}\\ {\mathrm{Variogram:\; Circular(100,1600,8000)}}\end{array}\right)$
 (2.2.2.9) 
Now that we have set the variogram, we can evaluate the Kriging interpolator at any point, and create a plot of the interpolator.
>

f_kr(4000, 2000);
plots:display(Ppt,
plot3d(f_kr, 6100 .. 6000, 6600 .. 6600, grid=[99,99]));

The statistics behind the method give us another output: the expected variance associated with a predicted value, which gives a measure for our confidence in that value. We can ask for that for one point, or show it as the surface color in a 3D plot.
>

f_kr(4000, 2000, output=variance);

${227.482481165330682}$
 (2.2.2.10) 
>

plots:display(Ppt,
plot3d(f_kr(x, y), x=6100 .. 6000, y=6600 .. 6600,
color = [f_kr(x, y, output=variance)/800, 0.8, 0.8, colortype=HSV]));

We see that areas close to many points give higher confidence than areas that are further away from most points.



Return to Index for Example Worksheets
