1GMTMATH(1) Generic Mapping Tools GMTMATH(1)
2
6 gmtmath - Reverse Polish Notation calculator for data tables
7
9 gmtmath [ -At_f(t).d ] [ -Ccols ] [ -Fcols ] [ -H[i][nrec] ] [ -I ] [
10 -M[i|o][flag] ] [ -Nn_col/t_col ] [ -Q ] [ -S[f|l] ] [
11 -Tt_min/t_max/t_inc[+]|tfile ] [ -V ] [
12 -b[i|o][s|S|d|D[ncol]|c[var1/...]] ] [ -f[i|o]colinfo ] operand [ oper‐
13 and ] OPERATOR [ operand ] OPERATOR ... = [ outfile ]
14
16 gmtmath will perform operations like add, subtract, multiply, and
17 divide on one or more table data files or constants using Reverse Pol‐
18 ish Notation (RPN) syntax (e.g., Hewlett-Packard calculator-style).
19 Arbitrarily complicated expressions may therefore be evaluated; the
20 final result is written to an output file [or standard output]. When
21 two data tables are on the stack, each element in file A is modified by
22 the corresponding element in file B. However, some operators only
23 require one operand (see below). If no data tables are used in the
24 expression then options -T, -N can be set (and optionally -b to indi‐
25 cate the data domain). If STDIN is given, <stdin> will be read and
26 placed on the stack as if a file with that content had been given on
27 the command line. By default, all columns except the "time" column are
28 operated on, but this can be changed (see -C).
29 operand
31 If operand can be opened as a file it will be read as an ASCII
32 (or binary, see -bi) table data file. If not a file, it is
33 interpreted as a numerical constant or a special symbol (see
34 below). The special argument STDIN means that stdin will be
35 read and placed on the stack; STDIN can appear more than once if
36 necessary.
37 outfile
39 The name of a table data file that will hold the final result.
40 If not given then the output is sent to stdout.
41 OPERATORS
43 Choose among the following 127 operators. "args" are the number
44 of input and output arguments.
45 Operator args Returns
47 ABS 1 1 abs (A).
49 ACOS 1 1 acos (A).
50 ACOSH 1 1 acosh (A).
51 ACOT 1 1 acot (A).
52 ACSC 1 1 acsc (A).
53 ADD 2 1 A + B.
54 AND 2 1 NaN if A and B == NaN, B if A == NaN, else A.
55 ASEC 1 1 asec (A).
56 ASIN 1 1 asin (A).
57 ASINH 1 1 asinh (A).
58 ATAN 1 1 atan (A).
59 ATAN2 2 1 atan2 (A, B).
60 ATANH 1 1 atanh (A).
61 BEI 1 1 bei (A).
62 BER 1 1 ber (A).
63 CEIL 1 1 ceil (A) (smallest integer >= A).
64 CHICRIT 2 1 Critical value for chi-squared-distribution, with
65 alpha = A and n = B.
66 CHIDIST 2 1 chi-squared-distribution P(chi2,n), with chi2 = A
67 and n = B.
68 COL 1 1 Places column A on the stack.
69 CORRCOEFF 2 1 Correlation coefficient r(A, B).
70 COS 1 1 cos (A) (A in radians).
71 COSD 1 1 cos (A) (A in degrees).
72 COSH 1 1 cosh (A).
73 COT 1 1 cot (A) (A in radians).
74 COTD 1 1 cot (A) (A in degrees).
75 CPOISS 2 1 Cumulative Poisson distribution F(x,lambda), with
76 x = A and lambda = B.
77 CSC 1 1 csc (A) (A in radians).
78 CSCD 1 1 csc (A) (A in degrees).
79 D2DT2 1 1 d^2(A)/dt^2 2nd derivative.
80 D2R 1 1 Converts Degrees to Radians.
81 DDT 1 1 d(A)/dt 1st derivative.
82 DILOG 1 1 dilog (A).
83 DIV 2 1 A / B.
84 DUP 1 2 Places duplicate of A on the stack.
85 EQ 2 1 1 if A == B, else 0.
86 ERF 1 1 Error function erf (A).
87 ERFC 1 1 Complementary Error function erfc (A).
88 ERFINV 1 1 Inverse error function of A.
89 EXCH 2 2 Exchanges A and B on the stack.
90 EXP 1 1 exp (A).
91 FACT 1 1 A! (A factorial).
92 FCRIT 3 1 Critical value for F-distribution, with alpha =
93 A, n1 = B, and n2 = C.
94 FDIST 3 1 F-distribution Q(F,n1,n2), with F = A, n1 = B,
95 and n2 = C.
96 FLIPUD 1 1 Reverse order of each column.
97 FLOOR 1 1 floor (A) (greatest integer <= A).
98 FMOD 2 1 A % B (remainder).
99 GE 2 1 1 if A >= B, else 0.
100 GT 2 1 1 if A > B, else 0.
101 HYPOT 2 1 hypot (A, B) = sqrt (A*A + B*B).
102 I0 1 1 Modified Bessel function of A (1st kind, order
103 0).
104 I1 1 1 Modified Bessel function of A (1st kind, order
105 1).
106 IN 2 1 Modified Bessel function of A (1st kind, order
107 B).
108 INT 1 1 Numerically integrate A.
109 INV 1 1 1 / A.
110 ISNAN 1 1 1 if A == NaN, else 0.
111 J0 1 1 Bessel function of A (1st kind, order 0).
112 J1 1 1 Bessel function of A (1st kind, order 1).
113 JN 2 1 Bessel function of A (1st kind, order B).
114 K0 1 1 Modified Kelvin function of A (2nd kind, order
115 0).
116 K1 1 1 Modified Bessel function of A (2nd kind, order
117 1).
118 KEI 1 1 kei (A).
119 KER 1 1 ker (A).
120 KN 2 1 Modified Bessel function of A (2nd kind, order
121 B).
122 KURT 1 1 Kurtosis of A.
123 LE 2 1 1 if A <= B, else 0.
124 LMSSCL 1 1 LMS scale estimate (LMS STD) of A.
125 LOG 1 1 log (A) (natural log).
126 LOG10 1 1 log10 (A) (base 10).
127 LOG1P 1 1 log (1+A) (accurate for small A).
128 LOG2 1 1 log2 (A) (base 2).
129 LOWER 1 1 The lowest (minimum) value of A.
130 LRAND 2 1 Laplace random noise with mean A and std. devia‐
131 tion B.
132 LSQFIT 1 0 Let current table be [A | b]; return least
133 squares solution x = A \ b.
134 LT 2 1 1 if A < B, else 0.
135 MAD 1 1 Median Absolute Deviation (L1 STD) of A.
136 MAX 2 1 Maximum of A and B.
137 MEAN 1 1 Mean value of A.
138 MED 1 1 Median value of A.
139 MIN 2 1 Minimum of A and B.
140 MODE 1 1 Mode value (Least Median of Squares) of A.
141 MUL 2 1 A * B.
142 NAN 2 1 NaN if A == B, else A.
143 NEG 1 1 -A.
144 NEQ 2 1 1 if A != B, else 0.
145 NRAND 2 1 Normal, random values with mean A and std. devia‐
146 tion B.
147 OR 2 1 NaN if A or B == NaN, else A.
148 PLM 3 1 Associated Legendre polynomial P(A) degree B
149 order C.
150 PLMg 3 1 Normalized associated Legendre polynomial P(A)
151 degree B order C (geophysical convention).
152 POP 1 0 Delete top element from the stack.
153 POW 2 1 A ^ B.
154 PQUANT 2 1 The B'th Quantile (0-100%) of A.
155 PSI 1 1 Psi (or Digamma) of A.
156 PV 3 1 Legendre function Pv(A) of degree v = real(B) +
157 imag(C).
158 QV 3 1 Legendre function Qv(A) of degree v = real(B) +
159 imag(C).
160 R2 2 1 R2 = A^2 + B^2.
161 R2D 1 1 Convert Radians to Degrees.
162 RAND 2 1 Uniform random values between A and B.
163 RINT 1 1 rint (A) (nearest integer).
164 ROOTS 2 1 Treats col A as f(t) = 0 and returns its roots.
165 ROTT 2 1 Rotate A by the (constant) shift B in the t-
166 direction.
167 SEC 1 1 sec (A) (A in radians).
168 SECD 1 1 sec (A) (A in degrees).
169 SIGN 1 1 sign (+1 or -1) of A.
170 SIN 1 1 sin (A) (A in radians).
171 SINC 1 1 sinc (A) (sin (pi*A)/(pi*A)).
172 SIND 1 1 sin (A) (A in degrees).
173 SINH 1 1 sinh (A).
174 SKEW 1 1 Skewness of A.
175 SQRT 1 1 sqrt (A).
176 STD 1 1 Standard deviation of A.
177 STEP 1 1 Heaviside step function H(A).
178 STEPT 1 1 Heaviside step function H(t-A).
179 SUB 2 1 A - B.
180 SUM 1 1 Cumulative sum of A.
181 TAN 1 1 tan (A) (A in radians).
182 TAND 1 1 tan (A) (A in degrees).
183 TANH 1 1 tanh (A).
184 TCRIT 2 1 Critical value for Student's t-distribution, with
185 alpha = A and n = B.
186 TDIST 2 1 Student's t-distribution A(t,n), with t = A, and
187 n = B.
188 TN 2 1 Chebyshev polynomial Tn(-1<A<+1) of degree B.
189 UPPER 1 1 The highest (maximum) value of A.
190 XOR 2 1 B if A == NaN, else A.
191 Y0 1 1 Bessel function of A (2nd kind, order 0).
192 Y1 1 1 Bessel function of A (2nd kind, order 1).
193 YN 2 1 Bessel function of A (2nd kind, order B).
194 ZCRIT 1 1 Critical value for the normal-distribution, with
195 alpha = A.
196 ZDIST 1 1 Cumulative normal-distribution C(x), with x = A.
197 SYMBOLS
199 The following symbols have special meaning:
200 PI 3.1415926...
202 E 2.7182818...
203 EULER 0.5772156...
204 TMIN Minimum t value
205 TMAX Maximum t value
206 TINC t increment
207 N The number of records
208 T Table with t-coordinates
209
211 -A Requires -N and will partially initialize a table with values
212 from the given file containing t and f(t) only. The t is placed
213 in column t_col while f(t) goes into column n_col - 1 (see -N).
214 -C Select the columns that will be operated on until next occur‐
216 rence of -C. List columns separated by commas; ranges like
217 1,3-5,7 are allowed. -C (no arguments) resets the default
218 action of using all columns except time column (see -N). -Ca
219 selects all columns, including time column, while -Cr reverses
220 (toggles) the current choices.
221 -F Give a comma-separated list of desired columns or ranges that
223 should be part of the output (0 is first column) [Default out‐
224 puts all columns].
225 -H Input file(s) has Header record(s). Number of header records
227 can be changed by editing your .gmtdefaults4 file. If used, GMT
228 default is 1 header record. Use -Hi if only input data should
229 have header records [Default will write out header records if
230 the input data have them]. Blank lines and lines starting with #
231 are always skipped.
232 -I Reverses the output row sequence from ascending time to descend‐
234 ing [ascending].
235 -M Multiple segment file(s). Segments are separated by a special
237 record. For ASCII files the first character must be flag
238 [Default is '>']. For binary files all fields must be NaN and
239 -b must set the number of output columns explicitly. By default
240 the -M setting applies to both input and output. Use -Mi and
241 -Mo to give separate settings.
242 -N Select the number of columns and the column number that contains
244 the "time" variable. Columns are numbered starting at 0 [2/0].
245 -Q Quick mode for scalar calculation. Shorthand for -Ca -N1/0
247 -T0/0/1.
248 -S Only report the first or last row of the results [Default is all
250 rows]. This is useful if you have computed a statistic (say the
251 MODE) and only want to report a single number instead of numer‐
252 ous records with identical values. Append l to get the last row
253 and f to get the first row only [Default].
254 -T Required when no input files are given. Sets the t-coordinates
256 of the first and last point and the equidistant sampling inter‐
257 val for the "time" column (see -N). Append + if you are speci‐
258 fying the number of equidistant points instead. If there is no
259 time column (only data columns), give -T with no arguments; this
260 also implies -Ca. Alternatively, give the name of a file whose
261 first column contains the desired t-coordinates which may be
262 irregular.
263 -V Selects verbose mode, which will send progress reports to stderr
265 [Default runs "silently"].
266 -bi Selects binary input. Append s for single precision [Default is
268 d (double)]. Uppercase S or D will force byte-swapping.
269 Optionally, append ncol, the number of columns in your binary
270 input file if it exceeds the columns needed by the program. Or
271 append c if the input file is netCDF. Optionally, append
272 var1/var2/... to specify the variables to be read.
273 -bo Selects binary output. Append s for single precision [Default
275 is d (double)]. Uppercase S or D will force byte-swapping.
276 Optionally, append ncol, the number of desired columns in your
277 binary output file. [Default is same as input, but see -F]
278
280 The ASCII output formats of numerical data are controlled by parameters
281 in your .gmtdefaults4 file. Longitude and latitude are formatted
282 according to OUTPUT_DEGREE_FORMAT, whereas other values are formatted
283 according to D_FORMAT. Be aware that the format in effect can lead to
284 loss of precision in the output, which can lead to various problems
285 downstream. If you find the output is not written with enough preci‐
286 sion, consider switching to binary output (-bo if available) or specify
287 more decimals using the D_FORMAT setting.
288
290 (1) The operators PLM and PLMg calculate the associated Legendre poly‐
291 nomial of degree L and order M in x which must satisfy -1 <= x <= +1
292 and 0 <= M <= L. x, L, and M are the three arguments preceeding the
293 operator. PLM is not normalized and includes the Condon-Shortley phase
294 (-1)^M. PLMg is normalized in the way that is most commonly used in
295 geophysics. The C-S phase can be added by using -M as argument. PLM
296 will overflow at higher degrees, whereas PLMg is stable until ultra
297 high degrees (at least 3000).
298 (2) Files that have the same names as some operators, e.g., ADD, SIGN,
300 =, etc. should be identified by prepending the current directory (i.e.,
301 ./LOG).
302 (3) The stack depth limit is hard-wired to 100.
304 (4) All functions expecting a positive radius (e.g., LOG, KEI, etc.)
306 are passed the absolute value of their argument.
307 (5) The DDT and D2DT2 functions only work on regularly spaced data.
309 (6) All derivatives are based on central finite differences, with natu‐
311 ral boundary conditions.
312 (7) ROOTS must be the last operator on the stack, only followed by =.
314
316 To take the square root of the content of the second data column being
317 piped through gmtmath by process1 and pipe it through a 3rd process,
318 use
319 process1 | gmtmath STDIN SQRT = | process3
321 To take log10 of the average of 2 data files, use
323 gmtmath file1.d file2.d ADD 0.5 MUL LOG10 = file3.d
325 Given the file samples.d, which holds seafloor ages in m.y. and
327 seafloor depth in m, use the relation depth(in m) = 2500 + 350 * sqrt
328 (age) to print the depth anomalies:
329 gmtmath samples.d T SQRT 350 MUL 2500 ADD SUB = | lpr
331 To take the average of columns 1 and 4-6 in the three data sets
333 sizes.1, sizes.2, and sizes.3, use
334 gmtmath -C1,4-6 sizes.1 sizes.2 ADD sizes.3 ADD 3 DIV = ave.d
336 To take the 1-column data set ages.d and calculate the modal value and
338 assign it to a variable, try
339 set mode_age = `gmtmath -S -T ages.d MODE =`
341 To evaluate the dilog(x) function for coordinates given in the file
343 t.d:
344 gmtmath -Tt.d T DILOG = dilog.d
346 To use gmtmath as a RPN Hewlett-Packard calculator on scalars (i.e., no
348 input files) and calculate arbitrary expressions, use the -Q option.
349 As an example, we will calculate the value of Kei (((1 + 1.75)/2.2) +
350 cos (60)) and store the result in the shell variable z:
351 set z = `gmtmath -Q 1 1.75 ADD 2.2 DIV 60 COSD ADD KEI =`
353 To use gmtmath as a general least squares equation solver, imagine that
355 the current table is the augmented matrix [ A | b ] and you want the
356 least squares solution x to the matrix equation A * x = b. The opera‐
357 tor LSQFIT does this; it is your job to populate the matrix correctly
358 first. The -A option will facilitate this. Suppose you have a 2-col‐
359 umn file ty.d with t and b(t) and you would like to fit a the model
360 y(t) = a + b*t + c*H(t-t0), where H is the Heaviside step function for
361 a given t0 = 1.55. Then, you need a 4-column augmented table loaded
362 with t in column 0 and your observed y(t) in column 3. The calculation
363 becomes
364 gmtmath -N4/1 -Aty.d -C0 1 ADD -C2 1.55 STEPT ADD -Ca LSQFIT = solu‐
366 tion.d
367 Note we use the -C option to select which columns we are working on,
369 then make active all the columns we need (here all of them, with -Ca)
370 before calling LSQFIT. The second and fourth columns (col numbers 1
371 and 3) are preloaded with t and y(t), respectively, the other columns
372 are zero. If you already have a precalculated table with the augmented
373 matrix [ A | b ] in a file (say lsqsys.d), the least squares solution
374 is simply
375 gmtmath -T lsqsys.d LSQFIT = solution.d
377
379 Abramowitz, M., and I. A. Stegun, 1964, Handbook of Mathematical Func‐
380 tions, Applied Mathematics Series, vol. 55, Dover, New York.
381 Holmes, S. A., and W. E. Featherstone, 2002, A unified approach to the
382 Clenshaw summation and the recursive computation of very high degree
383 and order normalised associated Legendre functions. Journal of
384 Geodesy, 76, 279-299.
385 Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
386 1992, Numerical Recipes, 2nd edition, Cambridge Univ., New York.
387 Spanier, J., and K. B. Oldman, 1987, An Atlas of Functions, Hemisphere
388 Publishing Corp.
389
391 GMT(1), grdmath(1)
392GMT 4.3.1 15 May 2008 GMTMATH(1)