1GRDMATH(1) Generic Mapping Tools GRDMATH(1)
2
6 grdmath - Reverse Polish Notation calculator for grid files
7
9 grdmath [ -F ] [ -Ixinc[unit][=|+][/yinc[unit][=|+]] ] [ -M ] [ -N ] [
10 -Rwest/east/south/north[r] ] [ -V ] [ -bi[s|S|d|D[ncol]|c[var1/...]] ]
11 [ -fcolinfo ] operand [ operand ] OPERATOR [ operand ] OPERATOR ... =
12 outgrdfile
13
15 grdmath will perform operations like add, subtract, multiply, and
16 divide on one or more grid files or constants using Reverse Polish
17 Notation (RPN) syntax (e.g., Hewlett-Packard calculator-style). Arbi‐
18 trarily complicated expressions may therefore be evaluated; the final
19 result is written to an output grid file. When two grids are on the
20 stack, each element in file A is modified by the corresponding element
21 in file B. However, some operators only require one operand (see
22 below). If no grid files are used in the expression then options -R,
23 -I must be set (and optionally -F). The expression = outgrdfile can
24 occur as many times as the depth of the stack allows.
25 operand
27 If operand can be opened as a file it will be read as a grid
28 file. If not a file, it is interpreted as a numerical constant
29 or a special symbol (see below).
30 outgrdfile
32 The name of a 2-D grid file that will hold the final result.
33 (See GRID FILE FORMATS below).
34 OPERATORS
36 Choose among the following 145 operators. "args" are the number
37 of input and output arguments.
38 Operator args Returns
40 ABS 1 1 abs (A).
42 ACOS 1 1 acos (A).
43 ACOSH 1 1 acosh (A).
44 ACOT 1 1 acot (A).
45 ACSC 1 1 acsc (A).
46 ADD 2 1 A + B.
47 AND 2 1 NaN if A and B == NaN, B if A == NaN, else A.
48 ASEC 1 1 asec (A).
49 ASIN 1 1 asin (A).
50 ASINH 1 1 asinh (A).
51 ATAN 1 1 atan (A).
52 ATAN2 2 1 atan2 (A, B).
53 ATANH 1 1 atanh (A).
54 BEI 1 1 bei (A).
55 BER 1 1 ber (A).
56 CAZ 2 1 Cartesian azimuth from grid nodes to stack x,y.
57 CBAZ 2 1 Cartesian backazimuth from grid nodes to stack
58 x,y.
59 CDIST 2 1 Cartesian distance between grid nodes and stack
60 x,y.
61 CEIL 1 1 ceil (A) (smallest integer >= A).
62 CHICRIT 2 1 Critical value for chi-squared-distribution, with
63 alpha = A and n = B.
64 CHIDIST 2 1 chi-squared-distribution P(chi2,n), with chi2 = A
65 and n = B.
66 CORRCOEFF 2 1 Correlation coefficient r(A, B).
67 COS 1 1 cos (A) (A in radians).
68 COSD 1 1 cos (A) (A in degrees).
69 COSH 1 1 cosh (A).
70 COT 1 1 cot (A) (A in radians).
71 COTD 1 1 cot (A) (A in degrees).
72 CPOISS 2 1 Cumulative Poisson distribution F(x,lambda), with
73 x = A and lambda = B.
74 CSC 1 1 csc (A) (A in radians).
75 CSCD 1 1 csc (A) (A in degrees).
76 CURV 1 1 Curvature of A (Laplacian).
77 D2DX2 1 1 d^2(A)/dx^2 2nd derivative.
78 D2DXY 1 1 d^2(A)/dxdy 2nd derivative.
79 D2DY2 1 1 d^2(A)/dy^2 2nd derivative.
80 D2R 1 1 Converts Degrees to Radians.
81 DDX 1 1 d(A)/dx Central 1st derivative.
82 DDY 1 1 d(A)/dy Central 1st derivative.
83 DILOG 1 1 dilog (A).
84 DIV 2 1 A / B.
85 DUP 1 2 Places duplicate of A on the stack.
86 EQ 2 1 1 if A == B, else 0.
87 ERF 1 1 Error function erf (A).
88 ERFC 1 1 Complementary Error function erfc (A).
89 ERFINV 1 1 Inverse error function of A.
90 EXCH 2 2 Exchanges A and B on the stack.
91 EXP 1 1 exp (A).
92 EXTREMA 1 1 Local Extrema: +2/-2 is max/min, +1/-1 is saddle
93 with max/min in x, 0 elsewhere.
94 FACT 1 1 A! (A factorial).
95 FCRIT 3 1 Critical value for F-distribution, with alpha =
96 A, n1 = B, and n2 = C.
97 FDIST 3 1 F-distribution Q(F,n1,n2), with F = A, n1 = B,
98 and n2 = C.
99 FLIPLR 1 1 Reverse order of values in each row.
100 FLIPUD 1 1 Reverse order of values in each column.
101 FLOOR 1 1 floor (A) (greatest integer <= A).
102 FMOD 2 1 A % B (remainder after truncated division).
103 GE 2 1 1 if A >= B, else 0.
104 GT 2 1 1 if A > B, else 0.
105 HYPOT 2 1 hypot (A, B) = sqrt (A*A + B*B).
106 I0 1 1 Modified Bessel function of A (1st kind, order
107 0).
108 I1 1 1 Modified Bessel function of A (1st kind, order
109 1).
110 IN 2 1 Modified Bessel function of A (1st kind, order
111 B).
112 INRANGE 3 1 1 if B <= A <= C, else 0.
113 INSIDE 1 1 1 when inside or on polygon(s) in A, else 0.
114 INV 1 1 1 / A.
115 ISNAN 1 1 1 if A == NaN, else 0.
116 J0 1 1 Bessel function of A (1st kind, order 0).
117 J1 1 1 Bessel function of A (1st kind, order 1).
118 JN 2 1 Bessel function of A (1st kind, order B).
119 K0 1 1 Modified Kelvin function of A (2nd kind, order
120 0).
121 K1 1 1 Modified Bessel function of A (2nd kind, order
122 1).
123 KEI 1 1 kei (A).
124 KER 1 1 ker (A).
125 KN 2 1 Modified Bessel function of A (2nd kind, order
126 B).
127 KURT 1 1 Kurtosis of A.
128 LDIST 1 1 Compute distance from lines in multi-segment
129 ASCII file A.
130 LE 2 1 1 if A <= B, else 0.
131 LMSSCL 1 1 LMS scale estimate (LMS STD) of A.
132 LOG 1 1 log (A) (natural log).
133 LOG10 1 1 log10 (A) (base 10).
134 LOG1P 1 1 log (1+A) (accurate for small A).
135 LOG2 1 1 log2 (A) (base 2).
136 LOWER 1 1 The lowest (minimum) value of A.
137 LRAND 2 1 Laplace random noise with mean A and std. devia‐
138 tion B.
139 LT 2 1 1 if A < B, else 0.
140 MAD 1 1 Median Absolute Deviation (L1 STD) of A.
141 MAX 2 1 Maximum of A and B.
142 MEAN 1 1 Mean value of A.
143 MED 1 1 Median value of A.
144 MIN 2 1 Minimum of A and B.
145 MOD 2 1 A mod B (remainder after floored division).
146 MODE 1 1 Mode value (Least Median of Squares) of A.
147 MUL 2 1 A * B.
148 NAN 2 1 NaN if A == B, else A.
149 NEG 1 1 -A.
150 NEQ 2 1 1 if A != B, else 0.
151 NOT 1 1 NaN if A == NaN, 1 if A == 0, else 0.
152 NRAND 2 1 Normal, random values with mean A and std. devia‐
153 tion B.
154 OR 2 1 NaN if A or B == NaN, else A.
155 PDIST 1 1 Compute distance from points in ASCII file A.
156 PLM 3 1 Associated Legendre polynomial P(A) degree B
157 order C.
158 PLMg 3 1 Normalized associated Legendre polynomial P(A)
159 degree B order C (geophysical convention).
160 POP 1 0 Delete top element from the stack.
161 POW 2 1 A ^ B.
162 PQUANT 2 1 The B'th Quantile (0-100%) of A.
163 PSI 1 1 Psi (or Digamma) of A.
164 PV 3 1 Legendre function Pv(A) of degree v = real(B) +
165 imag(C).
166 QV 3 1 Legendre function Qv(A) of degree v = real(B) +
167 imag(C).
168 R2 2 1 R2 = A^2 + B^2.
169 R2D 1 1 Convert Radians to Degrees.
170 RAND 2 1 Uniform random values between A and B.
171 RINT 1 1 rint (A) (nearest integer).
172 ROTX 2 1 Rotate A by the (constant) shift B in x-direc‐
173 tion.
174 ROTY 2 1 Rotate A by the (constant) shift B in y-direc‐
175 tion.
176 SAZ 2 1 Spherical azimuth from grid nodes to stack x,y.
177 SBAZ 2 1 Spherical backazimuth from grid nodes to stack
178 x,y.
179 SDIST 2 1 Spherical (Great circle) distance (in degrees)
180 between grid nodes and stack lon,lat (A, B).
181 SEC 1 1 sec (A) (A in radians).
182 SECD 1 1 sec (A) (A in degrees).
183 SIGN 1 1 sign (+1 or -1) of A.
184 SIN 1 1 sin (A) (A in radians).
185 SINC 1 1 sinc (A) (sin (pi*A)/(pi*A)).
186 SIND 1 1 sin (A) (A in degrees).
187 SINH 1 1 sinh (A).
188 SKEW 1 1 Skewness of A.
189 SQR 1 1 A^2.
190 SQRT 1 1 sqrt (A).
191 STD 1 1 Standard deviation of A.
192 STEP 1 1 Heaviside step function: H(A).
193 STEPX 1 1 Heaviside step function in x: H(x-A).
194 STEPY 1 1 Heaviside step function in y: H(y-A).
195 SUB 2 1 A - B.
196 TAN 1 1 tan (A) (A in radians).
197 TAND 1 1 tan (A) (A in degrees).
198 TANH 1 1 tanh (A).
199 TCRIT 2 1 Critical value for Student's t-distribution, with
200 alpha = A and n = B.
201 TDIST 2 1 Student's t-distribution A(t,n), with t = A, and
202 n = B.
203 TN 2 1 Chebyshev polynomial Tn(-1<t<+1,n), with t = A,
204 and n = B.
205 UPPER 1 1 The highest (maximum) value of A.
206 XOR 2 1 B if A == NaN, else A.
207 Y0 1 1 Bessel function of A (2nd kind, order 0).
208 Y1 1 1 Bessel function of A (2nd kind, order 1).
209 YLM 2 2 Re and Im orthonormalized spherical harmonics
210 degree A order B.
211 YLMg 2 2 Cos and Sin normalized spherical harmonics degree
212 A order B (geophysical convention).
213 YN 2 1 Bessel function of A (2nd kind, order B).
214 ZCRIT 1 1 Critical value for the normal-distribution, with
215 alpha = A.
216 ZDIST 1 1 Cumulative normal-distribution C(x), with x = A.
217 SYMBOLS
219 The following symbols have special meaning:
220 PI 3.1415926...
222 E 2.7182818...
223 EULER 0.5772156...
224 XMIN Minimum x value
225 XMAX Maximum x value
226 XINC x increment
227 NX The number of x nodes
228 YMIN Minimum y value
229 YMAX Maximum y value
230 YINC y increment
231 NY The number of y nodes
232 X Grid with x-coordinates
233 Y Grid with y-coordinates
234 Xn Grid with normalized [-1 to +1] x-coordinates
235 Yn Grid with normalized [-1 to +1] y-coordinates
236
238 -F Force pixel node registration [Default is gridline registra‐
239 tion]. (Node registrations are defined in GMT Cookbook Appendix
240 B on grid file formats.) Only used with -R -I.
241 -I x_inc [and optionally y_inc] is the grid spacing. Optionally,
243 append a suffix modifier. Geographical (degrees) coordinates:
244 Append m to indicate arc minutes or c to indicate arc seconds.
245 If one of the units e, k, i, or n is appended instead, the
246 increment is assumed to be given in meter, km, miles, or nauti‐
247 cal miles, respectively, and will be converted to the equivalent
248 degrees longitude at the middle latitude of the region (the con‐
249 version depends on ELLIPSOID). If /y_inc is given but set to 0
250 it will be reset equal to x_inc; otherwise it will be converted
251 to degrees latitude. All coordinates: If = is appended then the
252 corresponding max x (east) or y (north) may be slightly adjusted
253 to fit exactly the given increment [by default the increment may
254 be adjusted slightly to fit the given domain]. Finally, instead
255 of giving an increment you may specify the number of nodes
256 desired by appending + to the supplied integer argument; the
257 increment is then recalculated from the number of nodes and the
258 domain. The resulting increment value depends on whether you
259 have selected a gridline-registered or pixel-registered grid;
260 see Appendix B for details. Note: if -Rgrdfile is used then
261 grid spacing has already been initialized; use -I to override
262 the values.
263 -M By default any derivatives calculated are in z_units/ x(or
265 y)_units. However, the user may choose this option to convert
266 dx,dy in degrees of longitude,latitude into meters using a flat
267 Earth approximation, so that gradients are in z_units/meter.
268 -N Turn off strict domain match checking when multiple grids are
270 manipulated [Default will insist that each grid domain is within
271 1e-4 * grid_spacing of the domain of the first grid listed].
272 -R xmin, xmax, ymin, and ymax specify the Region of interest. For
274 geographic regions, these limits correspond to west, east,
275 south, and north and you may specify them in decimal degrees or
276 in [+-]dd:mm[:ss.xxx][W|E|S|N] format. Append r if lower left
277 and upper right map coordinates are given instead of w/e/s/n.
278 The two shorthands -Rg and -Rd stand for global domain (0/360
279 and -180/+180 in longitude respectively, with -90/+90 in lati‐
280 tude). Alternatively, specify the name of an existing grid file
281 and the -R settings (and grid spacing, if applicable) are copied
282 from the grid. For calendar time coordinates you may either
283 give (a) relative time (relative to the selected TIME_EPOCH and
284 in the selected TIME_UNIT; append t to -JX|x), or (b) absolute
285 time of the form [date]T[clock] (append T to -JX|x). At least
286 one of date and clock must be present; the T is always required.
287 The date string must be of the form [-]yyyy[-mm[-dd]] (Gregorian
288 calendar) or yyyy[-Www[-d]] (ISO week calendar), while the clock
289 string must be of the form hh:mm:ss[.xxx]. The use of delim‐
290 iters and their type and positions must be exactly as indicated
291 (however, input, output and plot formats are customizable; see
292 gmtdefaults).
293 -V Selects verbose mode, which will send progress reports to stderr
295 [Default runs "silently"].
296 -bi Selects binary input. Append s for single precision [Default is
298 d (double)]. Uppercase S or D will force byte-swapping.
299 Optionally, append ncol, the number of columns in your binary
300 input file if it exceeds the columns needed by the program. Or
301 append c if the input file is netCDF. Optionally, append
302 var1/var2/... to specify the variables to be read. The binary
303 input option only applies to the data files needed by operators
304 LDIST, PDIST, and INSIDE.
305 -f Special formatting of input and/or output columns (time or geo‐
307 graphical data). Specify i or o to make this apply only to
308 input or output [Default applies to both]. Give one or more
309 columns (or column ranges) separated by commas. Append T (abso‐
310 lute calendar time), t (relative time in chosen TIME_UNIT since
311 TIME_EPOCH), x (longitude), y (latitude), or f (floating point)
312 to each column or column range item. Shorthand -f[i|o]g means
313 -f[i|o]0x,1y (geographic coordinates).
314
316 (1) The operator SDIST calculates spherical distances between the (lon,
317 lat) point on the stack and all node positions in the grid. The grid
318 domain and the (lon, lat) point are expected to be in degrees. Simi‐
319 larly, the SAZ and SBAZ operators calculate spherical azimuth and back-
320 azimuths in degrees, respectively. A few operators (PDIST, LDIST, and
321 INSIDE) expects their argument to be a single file with points, lines,
322 or polygons, respectively. Be aware that LDIST in particular can be
323 slow for large grids and numerous line segments. Note: If the current
324 ELLIPSOID is not spherical then geodesics are used in the calculations.
325 (2) The operator PLM calculates the associated Legendre polynomial of
327 degree L and order M (0 <= M <= L), and its argument is the sine of the
328 latitude. PLM is not normalized and includes the Condon-Shortley phase
329 (-1)^M. PLMg is normalized in the way that is most commonly used in
330 geophysics. The C-S phase can be added by using -M as argument. PLM
331 will overflow at higher degrees, whereas PLMg is stable until ultra
332 high degrees (at least 3000).
333 (3) The operators YLM and YLMg calculate normalized spherical harmonics
335 for degree L and order M (0 <= M <= L) for all positions in the grid,
336 which is assumed to be in degrees. YLM and YLMg return two grids, the
337 real (cosine) and imaginary (sine) component of the complex spherical
338 harmonic. Use the POP operator (and EXCH) to get rid of one of them,
339 or save both by giving two consecutive = file.grd calls.
340 The orthonormalized complex harmonics YLM are most commonly used in
341 physics and seismology. The square of YLM integrates to 1 over a
342 sphere. In geophysics, YLMg is normalized to produce unit power when
343 averaging the cosine and sine terms (separately!) over a sphere (i.e.
344 their squares each integrate to 4 pi). The Condon-Shortley phase
345 (-1)^M is not included in YLM or YLMg, but it can be added by using -M
346 as argument.
347 (4) All the derivatives are based on central finite differences, with
349 natural boundary conditions.
350 (5) Files that have the same names as some operators, e.g., ADD, SIGN,
352 =, etc. should be identified by prepending the current directory (i.e.,
353 ./LOG).
354 (6) Piping of files is not allowed.
356 (7) The stack depth limit is hard-wired to 100.
358 (8) All functions expecting a positive radius (e.g., LOG, KEI, etc.)
360 are passed the absolute value of their argument.
361
363 Regardless of the precision of the input data, GMT programs that create
364 grid files will internally hold the grids in 4-byte floating point
365 arrays. This is done to conserve memory and furthermore most if not
366 all real data can be stored using 4-byte floating point values. Data
367 with higher precision (i.e., double precision values) will lose that
368 precision once GMT operates on the grid or writes out new grids. To
369 limit loss of precision when processing data you should always consider
370 normalizing the data prior to processing.
371
373 By default GMT writes out grid as single precision floats in a COARDS-
374 complaint netCDF file format. However, GMT is able to produce grid
375 files in many other commonly used grid file formats and also facili‐
376 tates so called "packing" of grids, writing out floating point data as
377 2- or 4-byte integers. To specify the precision, scale and offset, the
378 user should add the suffix =id[/scale/offset[/nan]], where id is a two-
379 letter identifier of the grid type and precision, and scale and offset
380 are optional scale factor and offset to be applied to all grid values,
381 and nan is the value used to indicate missing data. When reading
382 grids, the format is generally automatically recognized. If not, the
383 same suffix can be added to input grid file names. See grdreformat(1)
384 and Section 4.17 of the GMT Technical Reference and Cookbook for more
385 information.
386 When reading a netCDF file that contains multiple grids, GMT will read,
388 by default, the first 2-dimensional grid that can find in that file. To
389 coax GMT into reading another multi-dimensional variable in the grid
390 file, append ?varname to the file name, where varname is the name of
391 the variable. Note that you may need to escape the special meaning of ?
392 in your shell program by putting a backslash in front of it, or by
393 placing the filename and suffix between quotes or double quotes. The
394 ?varname suffix can also be used for output grids to specify a variable
395 name different from the default: "z". See grdreformat(1) and Section
396 4.18 of the GMT Technical Reference and Cookbook for more information,
397 particularly on how to read splices of 3-, 4-, or 5-dimensional grids.
398
400 When the output grid type is netCDF, the coordinates will be labeled
401 "longitude", "latitude", or "time" based on the attributes of the input
402 data or grid (if any) or on the -f or -R options. For example, both
403 -f0x -f1t and -R90w/90e/0t/3t will result in a longitude/time grid.
404 When the x, y, or z coordinate is time, it will be stored in the grid
405 as relative time since epoch as specified by TIME_UNIT and TIME_EPOCH
406 in the .gmtdefaults file or on the command line. In addition, the unit
407 attribute of the time variable will indicate both this unit and epoch.
408
410 To take log10 of the average of 2 files, use
411 grdmath file1.grd file2.grd ADD 0.5 MUL LOG10 = file3.grd
413 Given the file ages.grd, which holds seafloor ages in m.y., use the
415 relation depth(in m) = 2500 + 350 * sqrt (age) to estimate normal
416 seafloor depths:
417 grdmath ages.grd SQRT 350 MUL 2500 ADD = depths.grd
419 To find the angle a (in degrees) of the largest principal stress from
421 the stress tensor given by the three files s_xx.grd s_yy.grd, and
422 s_xy.grd from the relation tan (2*a) = 2 * s_xy / (s_xx - s_yy), use
423 grdmath 2 s_xy.grd MUL s_xx.grd s_yy.grd SUB DIV ATAN2 2 DIV = direc‐
425 tion.grd
426 To calculate the fully normalized spherical harmonic of degree 8 and
428 order 4 on a 1 by 1 degree world map, using the real amplitude 0.4 and
429 the imaginary amplitude 1.1:
430 grdmath -R0/360/-90/90 -I1 8 4 YML 1.1 MUL EXCH 0.4 MUL ADD = harm.grd
432 To extract the locations of local maxima that exceed 100 mGal in the
434 file faa.grd:
435 grdmath faa.grd DUP EXTREMA 2 EQ MUL DUP 100 GT MUL 0 NAN = z.grd
437 grd2xyz z.grd -S > max.xyz
438
440 Abramowitz, M., and I. A. Stegun, 1964, Handbook of Mathematical Func‐
441 tions, Applied Mathematics Series, vol. 55, Dover, New York.
442 Holmes, S. A., and W. E. Featherstone, 2002, A unified approach to the
443 Clenshaw summation and the recursive computation of very high degree
444 and order normalised associated Legendre functions. Journal of
445 Geodesy, 76, 279-299.
446 Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
447 1992, Numerical Recipes, 2nd edition, Cambridge Univ., New York.
448 Spanier, J., and K. B. Oldman, 1987, An Atlas of Functions, Hemisphere
449 Publishing Corp.
450
452 GMT(1), gmtmath(1), grd2xyz(1), grdedit(1), grdinfo(1), xyz2grd(1)
453GMT 4.5.6 10 Mar 2011 GRDMATH(1)