1GREENSPLINE(1) Generic Mapping Tools GREENSPLINE(1)
2
6 greenspline - Interpolate 1-D, 2-D, 3-D Cartesian or spherical surface
7 data using Green's function splines.
8
10 greenspline [ datafile(s) ] [ -A[1|2|3|4|5,]gradfile ] [ -Ccut[/file] ]
11 [ -Dmode ] [ -F ] [ -Ggrdfile ] [ -H[i][nrec] ] [ -Ixinc[yinc[zinc]] ]
12 [ -L ] [ -Nnodefile ] [ -Qaz|x/y/z ] [ -Rxmin/xmax[/ymin/ymax[/zminz‐
13 max]] ] [ -Sc|t|g|p|q[pars] ] [ -Tmaskgrid ] [ -V ] [ -:[i|o] ] [
14 -bi[s|S|d|D[ncol]|c[var1/...]] ] [ -bo[s|S|d|D[ncol]|c[var1/...]] ]
15
17 greenspline uses the Green's function G(x; x') for the chosen spline
18 and geometry to interpolate data at regular [or arbitrary] output loca‐
19 tions. Mathematically, the solution is composed as w(x) = sum {c(i)
20 G(x; x(i))}, for i = 1, n, the number of data points {x(i), w(i)}.
21 Once the n coefficients c(i) have been found then the sum can be evalu‐
22 ated at any output point x. Choose between ten minimum curvature, reg‐
23 ularized, or continuous curvature splines in tension for either 1-D,
24 2-D, or 3-D Cartesian coordinates or spherical surface coordinates.
25 After first removing a linear or planar trend (Cartesian geometries) or
26 mean value (spherical surface) and normalizing these residuals, the
27 least-squares matrix solution for the spline coefficients c(i) is found
28 by solving the n by n linear system w(j) = sum-over-i {c(i) G(x(j);
29 x(i))}, for j = 1, n; this solution yields an exact interpolation of
30 the supplied data points. Alternatively, you may choose to perform a
31 singular value decomposition (SVD) and eliminate the contribution from
32 the smallest eigenvalues; this approach yields an approximate solution.
33 Trends and scales are restored when evaluating the output.
34
36 datafile(s)
37 The name of one or more ASCII [or binary, see -bi] files holding
38 the x, w data points. If no file is given then we read standard
39 input instead.
40 -A The solution will partly be constrained by surface gradients v =
42 v*n, where v is the gradient magnitude and n its unit vector
43 direction. The gradient direction may be specified either by
44 Cartesian components (either unit vector n and magnitude v sepa‐
45 rately or gradient components v directly) or angles w.r.t. the
46 coordinate axes. Specify one of five input formats: 0: For 1-D
47 data there is no direction, just gradient magnitude (slope) so
48 the input format is x, gradient. Options 1-2 are for 2-D data
49 sets: 1: records contain x, y, azimuth, gradient (azimuth in
50 degrees is measured clockwise from the vertical (north)
51 [Default]). 2: records contain x, y, gradient, azimuth (azimuth
52 in degrees is measured clockwise from the vertical (north)).
53 Options 3-5 are for either 2-D or 3-D data: 3: records contain
54 x, direction(s), v (direction(s) in degrees are measured
55 counter-clockwise from the horizontal (and for 3-D the vertical
56 axis). 4: records contain x, v. 5: records contain x, n, v.
57 Append name of ASCII file with the surface gradients (following
58 a comma if a format is specified).
59 -C Find an approximate surface fit: Solve the linear system for the
61 spline coefficients by SVD and eliminate the contribution from
62 all eigenvalues whose ratio to the largest eigenvalue is less
63 than cut [Default uses Gauss-Jordan elimination to solve the
64 linear system and fit the data exactly]. Optionally, append
65 /file to save the eigenvalue ratios to the specified file for
66 further analysis. Finally, if a negative cut is given then
67 /file is required and execution will stop after saving the ei‐
68 genvalues, i.e., no surface output is produced.
69 -D Sets the distance flag that determines how we calculate dis‐
71 tances between data points. Select mode 0 for Cartesian 1-D
72 spline interpolation: -D0 means (x) in user units, Cartesian
73 distances, Select mode 1-3 for Cartesian 2-D surface spline
74 interpolation: -D1 means (x,y) in user units, Cartesian dis‐
75 tances, -D2 for (x,y) in degrees, flat Earth distances, and -D3
76 for (x,y) in degrees, spherical distances in km. Then, if
77 ELLIPSOID is spherical, we compute great circle arcs, otherwise
78 geodesics. Option mode = 4 applies to spherical surface spline
79 interpolation only: -D4 for (x,y) in degrees, use cosine of
80 great circle (or geodesic) arcs. Select mode 5 for Cartesian
81 3-D surface spline interpolation: -D5 means (x,y,z) in user
82 units, Cartesian distances.
83 -F Force pixel registration. [Default is gridline registration].
85 -G Name of resulting output file. (1) If options -R, -I, and pos‐
87 sibly -F are set we produce an equidistant output table. This
88 will be written to stdout unless -G is specified. Note: for 2-D
89 grids the -G option is required. (2) If option -T is selected
90 then -G is required and the output file is a 2-D binary grid
91 file. Applies to 2-D interpolation only. (3) If -N is selected
92 then the output is an ASCII (or binary; see -bo) table; if -G is
93 not given then this table is written to standard output.
94 Ignored if -C or -C0 is given.
95 -H Input file(s) has header record(s). If used, the default number
97 of header records is N_HEADER_RECS. Use -Hi if only input data
98 should have header records [Default will write out header
99 records if the input data have them]. Blank lines and lines
100 starting with # are always skipped.
101 -I Specify equidistant sampling intervals, on for each dimension,
103 separated by slashes.
104 -L Do not remove a linear (1-D) or planer (2-D) trend when -D
106 selects mode 0-3 [For those Cartesian cases a least-squares line
107 or plane is modeled and removed, then restored after fitting a
108 spline to the residuals]. However, in mixed cases with both
109 data values and gradients, or for spherical surface data, only
110 the mean data value is removed (and later and restored).
111 -N ASCII file with coordinates of desired output locations x in the
113 first column(s). The resulting w values are appended to each
114 record and written to the file given in -G [or stdout if not
115 specified]; see -bo for binary output instead. This option
116 eliminates the need to specify options -R, -I, and -F.
117 -Q Rather than evaluate the surface, take the directional deriva‐
119 tive in the az azimuth and return the magnitude of this deriva‐
120 tive instead. For 3-D interpolation, specify the three compo‐
121 nents of the desired vector direction (the vector will be nor‐
122 malized before use).
123 -R Specify the domain for an equidistant lattice where output pre‐
125 dictions are required. Requires -I and optionally -F.
126 1-D: Give xmin/xmax, the minimum and maximum x coordinates.
127 2-D: Give xmin/xmax/ymin/ymax, the minimum and maximum x and y
128 coordinates. These may be Cartesian or geographical. If geo‐
129 graphical, then west, east, south, and north specify the Region
130 of interest, and you may specify them in decimal degrees or in
131 [+-]dd:mm[:ss.xxx][W|E|S|N] format. The two shorthands -Rg and
132 -Rd stand for global domain (0/360 and -180/+180 in longitude
133 respectively, with -90/+90 in latitude).
134 3-D: Give xmin/xmax/ymin/ymax/zmin/zmax, the minimum and maxi‐
135 mum x, y and z coordinates. See the 2-D section if your hori‐
136 zontal coordinates are geographical; note the shorthands -Rg and
137 -Rd cannot be used if a 3-D domain is specified.
138 -S Select one of five different splines. The first two are used
140 for 1-D, 2-D, or 3-D Cartesian splines (see -D for discussion).
141 Note that all tension values are expected to be normalized ten‐
142 sion in the range 0 < t < 1: (c) Minimum curvature spline
143 [Sandwell, 1987], (t) Continuous curvature spline in tension
144 [Wessel and Bercovici, 1998]; append tension[/scale] with ten‐
145 sion in the 0-1 range and optionally supply a length scale
146 [Default is the average grid spacing]. The next is a 2-D or 3-D
147 spline: (r) Regularized spline in tension [Mitasova and Mitas,
148 1993]; again, append tension and optional scale. The last two
149 are spherical surface splines and both imply -D4 -fg: (p) Mini‐
150 mum curvature spline [Parker, 1994], (q) Continuous curvature
151 spline in tension [Wessel and Becker, 2008]; append tension.
152 The G(x; x') for the last method is slower to compute; by speci‐
153 fying OPT(SQ) you can speed up calculations by first pre-calcu‐
154 lating G(x; x') for a dense set of x values (e.g., 100,001 nodes
155 between -1 to +1) and store them in look-up tables. Optionally
156 append /N (an odd integer) to specify how many points in the
157 spline to set [100001]
158 -T For 2-D interpolation only. Only evaluate the solution at the
160 nodes in the maskgrid that are not equal to NaN. This option
161 eliminates the need to specify options -R, -I, and -F.
162 -V Selects verbose mode, which will send progress reports to stderr
164 [Default runs "silently"].
165 -bi Selects binary input. Append s for single precision [Default is
167 d (double)]. Uppercase S or D will force byte-swapping.
168 Optionally, append ncol, the number of columns in your binary
169 input file if it exceeds the columns needed by the program. Or
170 append c if the input file is netCDF. Optionally, append
171 var1/var2/... to specify the variables to be read. [Default is
172 2-4 input columns (x,w); the number depends on the chosen dimen‐
173 sion].
174 -bo Selects binary output. Append s for single precision [Default
176 is d (double)]. Uppercase S or D will force byte-swapping.
177 Optionally, append ncol, the number of desired columns in your
178 binary output file.
179
181 To resample the x,y Gaussian random data created by gmtmath and stored
182 in 1D.txt, requesting output every 0.1 step from 0 to 10, and using a
183 minimum cubic spline, try
184 gmtmath -T0/10/1 0 1 NRAND = 1D.txt
186 psxy -R0/10/-5/5 -JX6i/3i -B2f1/1 -Sc0.1 -Gblack 1D.txt -K > 1D.ps
187 greenspline 1D.txt -R0/10 -I0.1 -Sc -V | psxy -R -J -O -Wthin >> 1D.ps
188 To apply a spline in tension instead, using a tension of 0.7, try
190 psxy -R0/10/-5/5 -JX6i/3i -B2f1/1 -Sc0.1 -Gblack 1D.txt -K > 1Dt.ps
192 greenspline 1D.txt -R0/10 -I0.1 -St0.7 -V | psxy -R -J -O -Wthin >>
193 1Dt.ps
194
196 To make a uniform grid using the minimum curvature spline for the same
197 Cartesian data set from Davis (1986) that is used in the GMT Cookbook
198 example 16, try
199 greenspline table_5.11 -R0/6.5/-0.2/6.5 -I0.1 -Sc -V -D1 -GS1987.grd
201 psxy -R0/6.5/-0.2/6.5 -JX6i -B2f1 -Sc0.1 -Gblack table_5.11 -K > 2D.ps
202 grdcontour -JX6i -B2f1 -O -C25 -A50 S1987.grd >> 2D.ps
203 To use Cartesian splines in tension but only evaluate the solution
205 where the input mask grid is not NaN, try
206 greenspline table_5.11 -Tmask.grd -St0.5 -V -D1 -GWB1998.grd
208 To use Cartesian generalized splines in tension and return the magni‐
210 tude of the surface slope in the NW direction, try
211 greenspline table_5.11 -R0/6.5/-0.2/6.5 -I0.1 -Sr0.95 -V -D1 -Q-45
213 -Gslopes.grd Finally, to use Cartesian minimum curvature splines in
214 recovering a surface where the input data is a single surface value
215 (pt.d) and the remaining constraints specify only the surface slope and
216 direction (slopes.d), use
217 greenspline pt.d -R-3.2/3.2/-3.2/3.2 -I0.1 -Sc -V -D1 -A1,slopes.d
219 -Gslopes.grd
220
222 To create a uniform 3-D Cartesian grid table based on the data in ta‐
223 ble_5.23 in Davis (1986) that contains x,y,z locations and a measure of
224 uranium oxide concentrations (in percent), try
225 greenspline table_5.23 -R5/40/-5/10/5/16 -I0.25 -Sr0.85 -V -D5
227 -G3D_UO2.txt
228
230 To recreate Parker's [1994] example on a global 1x1 degree grid, assum‐
231 ing the data are in file mag_obs_1990.d, try
232 greenspline -V -Rg -fg -Sp -D3 -I1 -GP1994.grd mag_obs_1990.d
234 To do the same problem but applying tension and use pre-calculated
236 Green functions, use
237 greenspline -V -Rg -fg -SQ0.85 -D3 -I1 -GWB2008.grd mag_obs_1990.d
239
241 (1) For the Cartesian cases we use the free-space Green functions,
242 hence no boundary conditions are applied at the edges of the specified
243 domain. For most applications this is fine as the region typically is
244 arbitrarily set to reflect the extent of your data. However, if your
245 application requires particular boundary conditions then you may con‐
246 sider using surface instead.
247 (2) In all cases, the solution is obtained by inverting a n x n double
248 precision matrix for the Green function coefficients, where n is the
249 number of data constraints. Hence, your computer's memory may place
250 restrictions on how large data sets you can process with greenspline.
251 Pre-processing your data with blockmean, blockmedian, or blockmode is
252 recommended to avoid aliasing and may also control the size of n. For
253 information, if n = 1024 then only 8 Mb memory is needed, but for n =
254 10240 we need 800 Mb. Note that greenspline is fully 64-bit compliant
255 if compiled as such.
256 (3) The inversion for coefficients can become numerically unstable when
257 data neighbors are very close compared to the overall span of the data.
258 You can remedy this by pre-processing the data, e.g., by averaging
259 closely spaced neighbors. Alternatively, you can improve stability by
260 using the SVD solution and discard information associated with the
261 smallest eigenvalues (see -C).
262
264 Tension is generally used to suppress spurious oscillations caused by
265 the minimum curvature requirement, in particular when rapid gradient
266 changes are present in the data. The proper amount of tension can only
267 be determined by experimentation. Generally, very smooth data (such as
268 potential fields) do not require much, if any tension, while rougher
269 data (such as topography) will typically interpolate better with moder‐
270 ate tension. Make sure you try a range of values before choosing your
271 final result. Note: the regularized spline in tension is only stable
272 for a finite range of scale values; you must experiment to find the
273 valid range and a useful setting. For more information on tension see
274 the references below.
275
277 Davis, J. C., 1986, Statistics and Data Analysis in Geology, 2nd Edi‐
278 tion, 646 pp., Wiley, New York,
279 Mitasova, H., and L. Mitas, 1993, Interpolation by regularized spline
280 with tension: I. Theory and implementation, Math. Geol., 25, 641-655.
281 Parker, R. L., 1994, Geophysical Inverse Theory, 386 pp., Princeton
282 Univ. Press, Princeton, N.J.
283 Sandwell, D. T., 1987, Biharmonic spline interpolation of Geos-3 and
284 Seasat altimeter data, Geophys. Res. Lett., 14, 139-142.
285 Wessel, P., and D. Bercovici, 1998, Interpolation with splines in ten‐
286 sion: a Green's function approach, Math. Geol., 30, 77-93.
287 Wessel, P., and J. M. Becker, 2008, Interpolation using a generalized
288 Green's function for a spherical surface spline in tension, Geophys. J.
289 Int, 174, 21-28.
290 Wessel, P., 2009, A general-purpose Green's function interpolator, Com‐
291 puters & Geosciences, 35, 1247-1254, doi:10.1016/j.cageo.2008.08.012.
292
294 GMT(1), gmtmath(1), nearneighbor(1), psxy(1), surface(1), triangu‐
295 late(1), xyz2grd(1)
296GMT 4.5.6 10 Mar 2011 GREENSPLINE(1)