1GREENSPLINE(1) GMT GREENSPLINE(1)
2
6 greenspline - Interpolate using Green's functions for splines in 1-3
7 dimensions
8
10 greenspline [ table ] [ -Agradfile+f1|2|3|4|5 ] [
11 -C[n|r|v]value[+ffile] ] [ -Dmode ] [ -E[misfitfile] ] [ -Ggrdfile ]
12 [ -Ixinc[/yinc[/zinc]] ] [ -L ] [ -Nnodefile ] [ -Qaz|x/y/z ] [
13 -Rwest/east/south/north[/zmin/zmax][+r] ] [ -Sc|t|l|r|p|q[pars] ] [
14 -Tmaskgrid ] [ -V[level] ] [ -W[w]] [ -bbinary ] [ -dnodata ] [
15 -eregexp ] [ -fflags ] [ -hheaders ] [ -oflags ] [ -x[[-]n] ] [ -:[i|o]
16 ]
17 Note: No space is allowed between the option flag and the associated
19 arguments.
20
22 greenspline uses the Green's function G(x; x') for the chosen spline
23 and geometry to interpolate data at regular [or arbitrary] output loca‐
24 tions. Mathematically, the solution is composed as w(x) = sum {c(i)
25 G(x'; x(i))}, for i = 1, n, the number of data points {x(i), w(i)}.
26 Once the n coefficients c(i) have been found the sum can be evaluated
27 at any output point x. Choose between minimum curvature, regularized,
28 or continuous curvature splines in tension for either 1-D, 2-D, or 3-D
29 Cartesian coordinates or spherical surface coordinates. After first
30 removing a linear or planar trend (Cartesian geometries) or mean value
31 (spherical surface) and normalizing these residuals, the least-squares
32 matrix solution for the spline coefficients c(i) is found by solving
33 the n by n linear system w(j) = sum-over-i {c(i) G(x(j); x(i))}, for j
34 = 1, n; this solution yields an exact interpolation of the supplied
35 data points. Alternatively, you may choose to perform a singular value
36 decomposition (SVD) and eliminate the contribution from the smallest
37 eigenvalues; this approach yields an approximate solution. Trends and
38 scales are restored when evaluating the output.
39
41 None.
42
44 table The name of one or more ASCII [or binary, see -bi] files holding
45 the x, w data points. If no file is given then we read standard
46 input instead.
47 -Agradfile+f1|2|3|4|5
49 The solution will partly be constrained by surface gradients v =
50 v*n, where v is the gradient magnitude and n its unit vector
51 direction. The gradient direction may be specified either by
52 Cartesian components (either unit vector n and magnitude v sepa‐
53 rately or gradient components v directly) or angles w.r.t. the
54 coordinate axes. Append name of ASCII file with the surface gra‐
55 dients. Use +f to select one of five input formats: 0: For 1-D
56 data there is no direction, just gradient magnitude (slope) so
57 the input format is x, gradient. Options 1-2 are for 2-D data
58 sets: 1: records contain x, y, azimuth, gradient (azimuth in
59 degrees is measured clockwise from the vertical (north)
60 [Default]). 2: records contain x, y, gradient, azimuth (azimuth
61 in degrees is measured clockwise from the vertical (north)).
62 Options 3-5 are for either 2-D or 3-D data: 3: records contain
63 x, direction(s), v (direction(s) in degrees are measured
64 counter-clockwise from the horizontal (and for 3-D the vertical
65 axis). 4: records contain x, v. 5: records contain x, n, v.
66 -C[n|r|v]value[+ffile]
68 Find an approximate surface fit: Solve the linear system for the
69 spline coefficients by SVD and eliminate the contribution from
70 all eigenvalues whose ratio to the largest eigenvalue is less
71 than value [Default uses Gauss-Jordan elimination to solve the
72 linear system and fit the data exactly]. Optionally, append
73 +ffile to save the eigenvalue ratios to the specified file for
74 further analysis. Finally, if a negative value is given then
75 +ffile is required and execution will stop after saving the ei‐
76 genvalues, i.e., no surface output is produced. Specify -Cv to
77 use the largest eigenvalues needed to explain approximately
78 value % of the data variance. Specify -Cr to use the largest
79 eigenvalues needed to leave approximately value as the model
80 misfit. If value is not given then -W is required and we com‐
81 pute value as the rms of the data uncertainties. Alternatively,
82 use -Cn to select the value largest eigenvalues. If a file is
83 given with -Cv then we save the eigenvalues instead of the
84 ratios.
85 -Dmode Sets the distance flag that determines how we calculate dis‐
87 tances between data points. Select mode 0 for Cartesian 1-D
88 spline interpolation: -D0 means (x) in user units, Cartesian
89 distances, Select mode 1-3 for Cartesian 2-D surface spline
90 interpolation: -D1 means (x,y) in user units, Cartesian dis‐
91 tances, -D2 for (x,y) in degrees, Flat Earth distances, and -D3
92 for (x,y) in degrees, Spherical distances in km. Then, if
93 PROJ_ELLIPSOID is spherical, we compute great circle arcs, oth‐
94 erwise geodesics. Option mode = 4 applies to spherical surface
95 spline interpolation only: -D4 for (x,y) in degrees, use cosine
96 of great circle (or geodesic) arcs. Select mode 5 for Cartesian
97 3-D surface spline interpolation: -D5 means (x,y,z) in user
98 units, Cartesian distances.
99 E[misfitfile]
101 Evaluate the spline exactly at the input data locations and report
102 statistics of the misfit (mean, standard deviation, and rms).
103 Optionally, append a filename and we will write the data table, aug‐
104 mented by two extra columns holding the spline estimate and the mis‐
105 fit.
106 -Ggrdfile
108 Name of resulting output file. (1) If options -R, -I, and possi‐
109 bly -r are set we produce an equidistant output table. This will
110 be written to stdout unless -G is specified. Note: for 2-D grids
111 the -G option is required. (2) If option -T is selected then -G
112 is required and the output file is a 2-D binary grid file.
113 Applies to 2-D interpolation only. (3) If -N is selected then
114 the output is an ASCII (or binary; see -bo) table; if -G is not
115 given then this table is written to standard output. Ignored if
116 -C or -C0 is given.
117 -Ixinc[/yinc[/zinc]]
119 Specify equidistant sampling intervals, on for each dimension,
120 separated by slashes.
121 -L Do not remove a linear (1-D) or planer (2-D) trend when -D
123 selects mode 0-3 [For those Cartesian cases a least-squares line
124 or plane is modeled and removed, then restored after fitting a
125 spline to the residuals]. However, in mixed cases with both data
126 values and gradients, or for spherical surface data, only the
127 mean data value is removed (and later and restored).
128 -Nnodefile
130 ASCII file with coordinates of desired output locations x in the
131 first column(s). The resulting w values are appended to each
132 record and written to the file given in -G [or stdout if not
133 specified]; see -bo for binary output instead. This option elim‐
134 inates the need to specify options -R, -I, and -r.
135 -Qaz|x/y/z
137 Rather than evaluate the surface, take the directional deriva‐
138 tive in the az azimuth and return the magnitude of this deriva‐
139 tive instead. For 3-D interpolation, specify the three compo‐
140 nents of the desired vector direction (the vector will be nor‐
141 malized before use).
142 -Rxmin/xmax[/ymin/ymax[/zmin/zmax]]
144 Specify the domain for an equidistant lattice where output pre‐
145 dictions are required. Requires -I and optionally -r.
146 1-D: Give xmin/xmax, the minimum and maximum x coordinates.
148 2-D: Give xmin/xmax/ymin/ymax, the minimum and maximum x and y
150 coordinates. These may be Cartesian or geographical. If geo‐
151 graphical, then west, east, south, and north specify the Region
152 of interest, and you may specify them in decimal degrees or in
153 [±]dd:mm[:ss.xxx][W|E|S|N] format. The two shorthands -Rg and
154 -Rd stand for global domain (0/360 and -180/+180 in longitude
155 respectively, with -90/+90 in latitude).
156 3-D: Give xmin/xmax/ymin/ymax/zmin/zmax, the minimum and maximum
158 x, y and z coordinates. See the 2-D section if your horizontal
159 coordinates are geographical; note the shorthands -Rg and -Rd
160 cannot be used if a 3-D domain is specified.
161 -Sc|t|l|r|p|q[pars]
163 Select one of six different splines. The first two are used for
164 1-D, 2-D, or 3-D Cartesian splines (see -D for discussion). Note
165 that all tension values are expected to be normalized tension in
166 the range 0 < t < 1: (c) Minimum curvature spline [Sandwell,
167 1987], (t) Continuous curvature spline in tension [Wessel and
168 Bercovici, 1998]; append tension[/scale] with tension in the 0-1
169 range and optionally supply a length scale [Default is the aver‐
170 age grid spacing]. The next is a 1-D or 2-D spline: (l) Linear
171 (1-D) or Bilinear (2-D) spline; these produce output that do not
172 exceed the range of the given data. The next is a 2-D or 3-D
173 spline: (r) Regularized spline in tension [Mitasova and Mitas,
174 1993]; again, append tension and optional scale. The last two
175 are spherical surface splines and both imply -D4: (p) Minimum
176 curvature spline [Parker, 1994], (q) Continuous curvature spline
177 in tension [Wessel and Becker, 2008]; append tension. The G(x';
178 x') for the last method is slower to compute (a series solution)
179 so we pre-calculate values and use cubic spline interpolation
180 lookup instead. Optionally append +nN (an odd integer) to
181 change how many points to use in the spline setup [10001]. The
182 finite Legendre sum has a truncation error [1e-6]; you can lower
183 that by appending +elimit at the expense of longer run-time.
184 -Tmaskgrid
186 For 2-D interpolation only. Only evaluate the solution at the
187 nodes in the maskgrid that are not equal to NaN. This option
188 eliminates the need to specify options -R, -I, and -r.
189 -V[level] (more ...)
191 Select verbosity level [c].
192 -W[w] Data one-sigma uncertainties are provided in the last column.
194 We then compute weights that are inversely proportional to the
195 uncertainties. Append w if weights are given instead of uncer‐
196 tainties. This results in a weighted least squares fit. Note
197 that this only has an effect if -C is used. [Default uses no
198 weights or uncertainties].
199 -bi[ncols][t] (more ...)
201 Select native binary input. [Default is 2-4 input columns (x,w);
202 the number depends on the chosen dimension].
203 -bo[ncols][type] (more ...)
205 Select native binary output.
206 -d[i|o]nodata (more ...)
208 Replace input columns that equal nodata with NaN and do the
209 reverse on output.
210 -e[~]"pattern" | -e[~]/regexp/[i] (more ...)
212 Only accept data records that match the given pattern.
213 -f[i|o]colinfo (more ...)
215 Specify data types of input and/or output columns.
216 -h[i|o][n][+c][+d][+rremark][+rtitle] (more ...)
218 Skip or produce header record(s).
219 -icols[+l][+sscale][+ooffset][,...] (more ...)
221 Select input columns and transformations (0 is first column).
222 -ocols[,...] (more ...)
224 Select output columns (0 is first column).
225 -r (more ...)
227 Set pixel node registration [gridline].
228 -x[[-]n] (more ...)
230 Limit number of cores used in multi-threaded algorithms (OpenMP
231 required).
232 -^ or just -
234 Print a short message about the syntax of the command, then
235 exits (NOTE: on Windows just use -).
236 -+ or just +
238 Print an extensive usage (help) message, including the explana‐
239 tion of any module-specific option (but not the GMT common
240 options), then exits.
241 -? or no arguments
243 Print a complete usage (help) message, including the explanation
244 of all options, then exits.
245
247 To resample the x,y Gaussian random data created by gmtmath and stored
248 in 1D.txt, requesting output every 0.1 step from 0 to 10, and using a
249 minimum cubic spline, try
250 gmt math -T0/10/1 0 1 NRAND = 1D.txt
252 gmt psxy -R0/10/-5/5 -JX6i/3i -B2f1/1 -Sc0.1 -Gblack 1D.txt -K > 1D.ps
253 gmt greenspline 1D.txt -R0/10 -I0.1 -Sc -V | psxy -R -J -O -Wthin >> 1D.ps
254 To apply a spline in tension instead, using a tension of 0.7, try
256 gmt psxy -R0/10/-5/5 -JX6i/3i -B2f1/1 -Sc0.1 -Gblack 1D.txt -K > 1Dt.ps
258 gmt greenspline 1D.txt -R0/10 -I0.1 -St0.7 -V | psxy -R -J -O -Wthin >> 1Dt.ps
259
261 To make a uniform grid using the minimum curvature spline for the same
262 Cartesian data set from Davis (1986) that is used in the GMT Technical
263 Reference and Cookbook example 16, try
264 gmt greenspline table_5.11 -R0/6.5/-0.2/6.5 -I0.1 -Sc -V -D1 -GS1987.nc
266 gmt psxy -R0/6.5/-0.2/6.5 -JX6i -B2f1 -Sc0.1 -Gblack table_5.11 -K > 2D.ps
267 gmt grdcontour -JX6i -B2f1 -O -C25 -A50 S1987.nc >> 2D.ps
268 To use Cartesian splines in tension but only evaluate the solution
270 where the input mask grid is not NaN, try
271 gmt greenspline table_5.11 -Tmask.nc -St0.5 -V -D1 -GWB1998.nc
273 To use Cartesian generalized splines in tension and return the magni‐
275 tude of the surface slope in the NW direction, try
276 gmt greenspline table_5.11 -R0/6.5/-0.2/6.5 -I0.1 -Sr0.95 -V -D1 -Q-45 -Gslopes.nc
278 Finally, to use Cartesian minimum curvature splines in recovering a
280 surface where the input data is a single surface value (pt.txt) and the
281 remaining constraints specify only the surface slope and direction
282 (slopes.txt), use
283 gmt greenspline pt.txt -R-3.2/3.2/-3.2/3.2 -I0.1 -Sc -V -D1 -Aslopes.txt+f1 -Gslopes.nc
285
287 To create a uniform 3-D Cartesian grid table based on the data in ta‐
288 ble_5.23 in Davis (1986) that contains x,y,z locations and a measure of
289 uranium oxide concentrations (in percent), try
290 gmt greenspline table_5.23 -R5/40/-5/10/5/16 -I0.25 -Sr0.85 -V -D5 -G3D_UO2.txt
292
294 To recreate Parker's [1994] example on a global 1x1 degree grid, assum‐
295 ing the data are in file mag_obs_1990.txt, try
296 greenspline -V -Rg -Sp -D3 -I1 -GP1994.nc mag_obs_1990.txt
298 To do the same problem but applying tension of 0.85, use
300 greenspline -V -Rg -Sq0.85 -D3 -I1 -GWB2008.nc mag_obs_1990.txt
302
304 1. For the Cartesian cases we use the free-space Green functions, hence
305 no boundary conditions are applied at the edges of the specified
306 domain. For most applications this is fine as the region typically
307 is arbitrarily set to reflect the extent of your data. However, if
308 your application requires particular boundary conditions then you
309 may consider using surface instead.
310 2. In all cases, the solution is obtained by inverting a n x n double
312 precision matrix for the Green function coefficients, where n is the
313 number of data constraints. Hence, your computer's memory may place
314 restrictions on how large data sets you can process with green‐
315 spline. Pre-processing your data with doc:blockmean, doc:blockme‐
316 dian, or doc:blockmode is recommended to avoid aliasing and may also
317 control the size of n. For information, if n = 1024 then only 8 Mb
318 memory is needed, but for n = 10240 we need 800 Mb. Note that green‐
319 spline is fully 64-bit compliant if compiled as such. For spherical
320 data you may consider decimating using doc:gmtspatial nearest neigh‐
321 bor reduction.
322 3. The inversion for coefficients can become numerically unstable when
324 data neighbors are very close compared to the overall span of the
325 data. You can remedy this by pre-processing the data, e.g., by
326 averaging closely spaced neighbors. Alternatively, you can improve
327 stability by using the SVD solution and discard information associ‐
328 ated with the smallest eigenvalues (see -C).
329 4. The series solution implemented for -Sq was developed by Robert L.
331 Parker, Scripps Institution of Oceanography, which we gratefully
332 acknowledge.
333 5. If you need to fit a certain 1-D spline through your data points you
335 may wish to consider sample1d instead. It will offer traditional
336 splines with standard boundary conditions (such as the natural cubic
337 spline, which sets the curvatures at the ends to zero). In con‐
338 trast, greenspline's 1-D spline, as is explained in note 1, does not
339 specify boundary conditions at the end of the data domain.
340
342 Tension is generally used to suppress spurious oscillations caused by
343 the minimum curvature requirement, in particular when rapid gradient
344 changes are present in the data. The proper amount of tension can only
345 be determined by experimentation. Generally, very smooth data (such as
346 potential fields) do not require much, if any tension, while rougher
347 data (such as topography) will typically interpolate better with moder‐
348 ate tension. Make sure you try a range of values before choosing your
349 final result. Note: the regularized spline in tension is only stable
350 for a finite range of scale values; you must experiment to find the
351 valid range and a useful setting. For more information on tension see
352 the references below.
353
355 Davis, J. C., 1986, Statistics and Data Analysis in Geology, 2nd Edi‐
356 tion, 646 pp., Wiley, New York,
357 Mitasova, H., and L. Mitas, 1993, Interpolation by regularized spline
359 with tension: I. Theory and implementation, Math. Geol., 25, 641-655.
360 Parker, R. L., 1994, Geophysical Inverse Theory, 386 pp., Princeton
362 Univ. Press, Princeton, N.J.
363 Sandwell, D. T., 1987, Biharmonic spline interpolation of Geos-3 and
365 Seasat altimeter data, Geophys. Res. Lett., 14, 139-142.
366 Wessel, P., and D. Bercovici, 1998, Interpolation with splines in ten‐
368 sion: a Green's function approach, Math. Geol., 30, 77-93.
369 Wessel, P., and J. M. Becker, 2008, Interpolation using a generalized
371 Green's function for a spherical surface spline in tension, Geophys. J.
372 Int, 174, 21-28.
373 Wessel, P., 2009, A general-purpose Green's function interpolator, Com‐
375 puters & Geosciences, 35, 1247-1254, doi:10.1016/j.cageo.2008.08.012.
376
378 gmt, gmtmath, nearneighbor, psxy, sample1d, sphtriangulate, surface,
379 triangulate, xyz2grd
380
382 2019, P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe
3835.4.5 Feb 24, 2019 GREENSPLINE(1)