pscontour (for contouring) and triangulate
(for gridding) use the simplest method of interpolating
data: a Delaunay triangulation (see Example 12) which
forms 
as a union of planar triangular facets.
One advantage of this method is that it will not extrapolate
beyond the convex hull of the input (x, y)
data. Another is that it will not estimate a z value
above or below the local bounds on any triangle.
A disadvantage is that the surface is not
differentiable, but has sharp kinks at triangle edges and
thus also along contours. This may not look physically
reasonable, but it can be filtered later (last panel below).
surface can be used to generate a higher-order
(smooth and differentiable) interpolation of onto
a grid, after which the grid may be illustrated (grdcontour,
grdimage, grdview). surface will interpolate
to all (x, y) points in a rectangular region, and thus
will extrapolate beyond the convex hull of the data. However,
this can be masked out in various ways (see Example 15).
A more serious objection is that surface may estimate
z values outside the local range of the data (note area
near x = 0.8, y = 5.3). This commonly happens when
the default tension value of zero is used to create a ``minimum
curvature'' (most smooth) interpolant. surface can be
used with non-zero tension to partially overcome this problem.
The limiting value 
should approximate the triangulation,
while a value between 0 and 1 may yield a good compromise between
the above two cases. A value of 0.5 is shown here
(Figure 7.16). A side
effect of the tension is that it tends to make the contours turn
near the edges of the domain so that they approach the edge from
a perpendicular direction. A solution is to use surface
in a larger area and then use grdcut to cut out the desired
smaller area. Another way to achieve a compromise is to
interpolate the data to a grid and then filter the grid using
grdfft or grdfilter. The latter can handle grids
containing ``NaN'' values and it can do median and mode filters
as well as convolutions. Shown here is triangulate followed
by grdfilter. Note that the filter has done some
extrapolation beyond the convex hull of the original x, y
values. The ``best'' smooth approximation of depends
on the errors in the data and the physical laws obeyed by z.
GMT cannot always do the ``best'' thing but it offers great
flexibility through its combinations of tools. We illustrate all
four solutions using a cpt file that contains color fills,
patterns, and a ``skip slice'' request for 
.
#!/bin/csh # GMT EXAMPLE 16 # # $Id: job16.csh,v 1.7 2004/04/13 21:32:27 pwessel Exp $ # # Purpose: Illustrates interpolation methods using same data as Example 12. # GMT progs: gmtset, grdview, grdfilter, pscontour, psscale, pstext, surface, triangulate # Unix progs: echo, rm # # First make a cpt file as in example 12: # #set z = `minmax table_5.11 -C -I25` #makecpt -Crainbow -T$z[5]/$z[6]/25 >! ex16.cpt #Hand edit to add patterns and skips # # Now illustrate various means of contouring, using triangulate and surface. # gmtset ANNOT_FONT_SIZE_PRIMARY 9 # pscontour -R0/6.5/-0.2/6.5 -Jx0.45i -P -K -Y5.5i -Ba2f1WSne table_5.11 -Cex16.cpt -I \ >! example_16.ps echo "3.25 7 18 0 4 CB pscontour (triangulate)" | pstext -R -J -O -K -N >> example_16.ps # surface table_5.11 -R -I0.1 -Graws0.grd grdview raws0.grd -R -J -Ba2f1WSne -Cex16.cpt -Qs -O -K -X3.5i >> example_16.ps echo "3.25 7 18 0 4 CB surface (tension = 0)" | pstext -R -J -O -K -N >> example_16.ps # surface table_5.11 -R -I0.1 -Graws5.grd -T0.5 grdview raws5.grd -R -J -Ba2f1WSne -Cex16.cpt -Qs -O -K -Y-3.75i -X-3.5i >> example_16.ps echo "3.25 7 18 0 4 CB surface (tension = 0.5)" | pstext -R -J -O -K -N >> example_16.ps # triangulate table_5.11 -Grawt.grd -R -I0.1 > /dev/null grdfilter rawt.grd -Gfiltered.grd -D0 -Fc1 grdview filtered.grd -R -J -Ba2f1WSne -Cex16.cpt -Qs -O -K -X3.5i >> example_16.ps echo "3.25 7 18 0 4 CB triangulate @~\256@~ grdfilter" | pstext -R -J -O -K -N >> example_16.ps echo "3.2125 7.5 32 0 4 CB Gridding of Data" | pstext -R0/10/0/10 -Jx1i -O -K -N -X-3.5i \ >> example_16.ps psscale -D3.25i/0.35i/5i/0.25ih -Cex16.cpt -O -U"Example 16 in Cookbook" -Y-0.75i >> example_16.ps # \rm -f *.grd .gmt*
