Some variations on Matlab gridding

Informal notes by S. Pierce

These examples attempt to reconstruct the function z = x * exp ( -x^2 - y^2 ) using random observation locations. The test data are available in a mat file here: testdat.mat.

Native matlab:
matlab griddata with default method (not generally recommended)
>> zi = griddata(x,y,z,xi,yi);
matlab griddata using option 'v4', a minimum curvature spline (Sandwell, 1987)
>> zi = griddata(x,y,z,xi,yi,'v4');

Some good freely available alternatives:
spline2d regspline2d
spline2d with tension=0.5
>> zi = spline2d(xi,yi,x,y,z,0.5);
regspline2d with tension=0.5
>> zi = regspline2d(xi,yi,x,y,z,0.5);
spline2d with tension=0.2
>> zi = spline2d(xi,yi,x,y,z,0.2);
regspline2d with tension=0.2
>> zi = regspline2d(xi,yi,x,y,z,0.2);
spline2d with tension=0.01
>> zi = spline2d(xi,yi,x,y,z,0.01);
regspline2d with tension=0.01
>> zi = regspline2d(xi,yi,x,y,z,0.01);
spline2d uses the method of Wessel & Bercovici (1998),
and is available here: spline2d.m
(modified version of Wessel TSPLINE original).
regspline2d uses the method of Mitasova & Mitas (1993),
and is available here: regspline2d.m.

To experiment on your own, and create similar figures (eg.):

Notes

The minimum curvature method (griddata 'v4') can produce spurious oscillations with some data, and it has no adjustment parameter. In this example, it is able to reproduce the original function with rmsfit = 1.04e-4.

The "tension" methods can effectively suppress spurious oscillations. The appropriate amount of tension is best determined with some experimentation. These routines use a normalized tension (0 <= t < 1).

The spline2d (Wessel & Bercovici, 1998), a favorite method of GMT package users, outperforms the griddata 'v4' in this case if the tension parameter is <=0.01 or so, with rmsfit = 0.70e-4.

The regspline2d (Mitasova & Mitas, 1993) can also produce excellent results, and is a method used within the GRASS GIS package. For tension parameters in the range 0.14-0.3, the rmsfit is better in this example than the other methods, achieving a minimum rmsfit=0.13e-4. But a cautionary note: this method is only stable for a certain range of tension values (different for each case). Outside of this, it becomes clearly unstable (see regspline2d example above with tension=0.01).

Another great option, when some smoothing on chosen length scales is desirable:
barnes
barnes with length scales xr=1.0 and yr=1.0
>> zi = barnes(x,y,z,xi,yi,1.0,1.0);
barnes with length scales xr=0.5 and yr=0.5
>> zi = barnes(x,y,z,xi,yi,0.5,0.5);
barnes uses Barnes objective analysis (Barnes, 1994), and is available here: barnes.m. See
>> help barnes
for more information.

References

Barnes, S. L. (1994) Applications of the Barnes objective analysis scheme. Part I: effects of undersampling, wave position, and station randomness. J. of Atmos. and Oceanic Tech., 11, 1433-1448.

Mitasova, H., and L. Mitas (1993) Interpolation by regularized spline with tension: I. Theory and implementation, Math. Geol., 25, 641-655.

Sandwell, D. T. (1987) Biharmonic spline interpolation of Geos-3 and Seasat altimeter data, Geophys. Res. Lett., 14, 139-142.

Wessel, P., and D. Bercovici (1998) Interpolation with splines in tension: a Green's function approach, Math. Geol., 30, 77-93.