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 142 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 1st derivative.
82 DDY 1 1 d(A)/dy 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).
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 MODE 1 1 Mode value (Least Median of Squares) of A.
146 MUL 2 1 A * B.
147 NAN 2 1 NaN if A == B, else A.
148 NEG 1 1 -A.
149 NEQ 2 1 1 if A != B, else 0.
150 NRAND 2 1 Normal, random values with mean A and std. devia‐
151 tion B.
152 OR 2 1 NaN if A or B == NaN, else A.
153 PDIST 1 1 Compute distance from points in ASCII file A.
154 PLM 3 1 Associated Legendre polynomial P(A) degree B
155 order C.
156 PLMg 3 1 Normalized associated Legendre polynomial P(A)
157 degree B order C (geophysical convention).
158 POP 1 0 Delete top element from the stack.
159 POW 2 1 A ^ B.
160 PQUANT 2 1 The B'th Quantile (0-100%) of A.
161 PSI 1 1 Psi (or Digamma) of A.
162 PV 3 1 Legendre function Pv(A) of degree v = real(B) +
163 imag(C).
164 QV 3 1 Legendre function Qv(A) of degree v = real(B) +
165 imag(C).
166 R2 2 1 R2 = A^2 + B^2.
167 R2D 1 1 Convert Radians to Degrees.
168 RAND 2 1 Uniform random values between A and B.
169 RINT 1 1 rint (A) (nearest integer).
170 ROTX 2 1 Rotate A by the (constant) shift B in x-direc‐
171 tion.
172 ROTY 2 1 Rotate A by the (constant) shift B in y-direc‐
173 tion.
174 SAZ 2 1 Spherical azimuth from grid nodes to stack x,y.
175 SBAZ 2 1 Spherical backazimuth from grid nodes to stack
176 x,y.
177 SDIST 2 1 Spherical (Great circle) distance (in degrees)
178 between grid nodes and stack lon,lat (A, B).
179 SEC 1 1 sec (A) (A in radians).
180 SECD 1 1 sec (A) (A in degrees).
181 SIGN 1 1 sign (+1 or -1) of A.
182 SIN 1 1 sin (A) (A in radians).
183 SINC 1 1 sinc (A) (sin (pi*A)/(pi*A)).
184 SIND 1 1 sin (A) (A in degrees).
185 SINH 1 1 sinh (A).
186 SKEW 1 1 Skewness of A.
187 SQRT 1 1 sqrt (A).
188 STD 1 1 Standard deviation of A.
189 STEP 1 1 Heaviside step function: H(A).
190 STEPX 1 1 Heaviside step function in x: H(x-A).
191 STEPY 1 1 Heaviside step function in y: H(y-A).
192 SUB 2 1 A - B.
193 TAN 1 1 tan (A) (A in radians).
194 TAND 1 1 tan (A) (A in degrees).
195 TANH 1 1 tanh (A).
196 TCRIT 2 1 Critical value for Student's t-distribution, with
197 alpha = A and n = B.
198 TDIST 2 1 Student's t-distribution A(t,n), with t = A, and
199 n = B.
200 TN 2 1 Chebyshev polynomial Tn(-1<t<+1,n), with t = A,
201 and n = B.
202 UPPER 1 1 The highest (maximum) value of A.
203 XOR 2 1 B if A == NaN, else A.
204 Y0 1 1 Bessel function of A (2nd kind, order 0).
205 Y1 1 1 Bessel function of A (2nd kind, order 1).
206 YLM 2 2 Re and Im orthonormalized spherical harmonics
207 degree A order B.
208 YLMg 2 2 Cos and Sin normalized spherical harmonics degree
209 A order B (geophysical convention).
210 YN 2 1 Bessel function of A (2nd kind, order B).
211 ZCRIT 1 1 Critical value for the normal-distribution, with
212 alpha = A.
213 ZDIST 1 1 Cumulative normal-distribution C(x), with x = A.
214 SYMBOLS
216 The following symbols have special meaning:
217 PI 3.1415926...
219 E 2.7182818...
220 EULER 0.5772156...
221 XMIN Minimum x value
222 XMAX Maximum x value
223 XINC x increment
224 NX The number of x nodes
225 YMIN Minimum y value
226 YMAX Maximum y value
227 YINC y increment
228 NY The number of y nodes
229 X Grid with x-coordinates
230 Y Grid with y-coordinates
231 Xn Grid with normalized [-1 to +1] x-coordinates
232 Yn Grid with normalized [-1 to +1] y-coordinates
233
235 -F Force pixel node registration [Default is gridline registra‐
236 tion]. (Node registrations are defined in GMT Cookbook Appendix
237 B on grid file formats.) Only used with -R, -I.
238 -I x_inc [and optionally y_inc] is the grid spacing. Optionally,
240 append a suffix modifier. Geographical (degrees) coordinates:
241 Append m to indicate arc minutes or c to indicate arc seconds.
242 If one of the units e, k, i, or n is appended instead, the
243 increment is assumed to be given in meter, km, miles, or nauti‐
244 cal miles, respectively, and will be converted to the equivalent
245 degrees longitude at the middle latitude of the region (the con‐
246 version depends on ELLIPSOID). If /y_inc is given but set to 0
247 it will be reset equal to x_inc; otherwise it will be converted
248 to degrees latitude. All coordinates: If = is appended then the
249 corresponding max x (east) or y (north) may be slightly adjusted
250 to fit exactly the given increment [by default the increment may
251 be adjusted slightly to fit the given domain]. Finally, instead
252 of giving an increment you may specify the number of nodes
253 desired by appending + to the supplied integer argument; the
254 increment is then recalculated from the number of nodes and the
255 domain. The resulting increment value depends on whether you
256 have selected a gridline-registered or pixel-registered grid;
257 see Appendix B for details.
258 -M By default any derivatives calculated are in z_units/ x(or
260 y)_units. However, the user may choose this option to convert
261 dx,dy in degrees of longitude,latitude into meters using a falt
262 Earth approximation, so that gradients are in z_units/meter.
263 -N Turn off strict domain match checking when multiple grids are
265 manipulated [Default will insist that each grid domain is within
266 1e-4 * grid_spacing of the domain of the first grid listed].
267 -R xmin, xmax, ymin, and ymax specify the Region of interest. For
269 geographic regions, these limits correspond to west, east,
270 south, and north and you may specify them in decimal degrees or
271 in [+-]dd:mm[:ss.xxx][W|E|S|N] format. Append r if lower left
272 and upper right map coordinates are given instead of w/e/s/n.
273 The two shorthands -Rg and -Rd stand for global domain (0/360
274 and -180/+180 in longitude respectively, with -90/+90 in lati‐
275 tude). For calendar time coordinates you may either give (a)
276 relative time (relative to the selected TIME_EPOCH and in the
277 selected TIME_UNIT; append t to -JX|x), or (b) absolute time of
278 the form [date]T[clock] (append T to -JX|x). At least one of
279 date and clock must be present; the T is always required. The
280 date string must be of the form [-]yyyy[-mm[-dd]] (Gregorian
281 calendar) or yyyy[-Www[-d]] (ISO week calendar), while the clock
282 string must be of the form hh:mm:ss[.xxx]. The use of delim‐
283 iters and their type and positions must be exactly as indicated
284 (however, input, output and plot formats are customizable; see
285 gmtdefaults).
286 -V Selects verbose mode, which will send progress reports to stderr
288 [Default runs "silently"].
289 -f Special formatting of input and/or output columns (time or geo‐
291 graphical data). Specify i or o to make this apply only to
292 input or output [Default applies to both]. Give one or more
293 columns (or column ranges) separated by commas. Append T (abso‐
294 lute calendar time), t (relative time in chosen TIME_UNIT since
295 TIME_EPOCH), x (longitude), y (latitude), or f (floating point)
296 to each column or column range item. Shorthand -f[i|o]g means
297 -f[i|o]0x,1y (geographic coordinates).
298
300 (1) The operator SDIST calculates spherical distances between the (lon,
301 lat) point on the stack and all node positions in the grid. The grid
302 domain and the (lon, lat) point are expected to be in degrees. Simi‐
303 larly, the SAZ and SBAZ operators calculate spherical azimuth and back‐
304 azimuths in degrees, respectively. Note: If the current ELLIPSOID is
305 not Sphere then geodesics are used in the calculations.
306 (2) The operator PLM calculates the associated Legendre polynomial of
308 degree L and order M (0 <= M <= L), and its argument is the sine of the
309 latitude. PLM is not normalized and includes the Condon-Shortley phase
310 (-1)^M. PLMg is normalized in the way that is most commonly used in
311 geophysics. The C-S phase can be added by using -M as argument. PLM
312 will overflow at higher degrees, whereas PLMg is stable until ultra
313 high degrees (at least 3000).
314 (3) The operators YLM and YLMg calculate normalized spherical harmonics
316 for degree L and order M (0 <= M <= L) for all positions in the grid,
317 which is assumed to be in degrees. YLM and YLMg return two grids, the
318 real (cosine) and imaginary (sine) component of the complex spherical
319 harmonic. Use the POP operator (and EXCH) to get rid of one of them,
320 or save both by giving two consequtive = file.grd calls.
321 The orthonormalized complex harmonics YLM are most commonly used in
322 physics and seismology. The square of YLM integrates to 1 over a
323 sphere. In geophysics, YLMg is normalized to produce unit power when
324 averaging the cosine and sine terms (separately!) over a sphere (i.e.
325 their squares each integrate to 4 pi). The Condon-Shortley phase
326 (-1)^M is not included in YLM or YLMg, but it can be added by using -M
327 as argument.
328 (4) All the derivatives are based on central finite differences, with
330 natural boundary conditions.
331 (5) Files that have the same names as some operators, e.g., ADD, SIGN,
333 =, etc. should be identified by prepending the current directory (i.e.,
334 ./LOG).
335 (6) Piping of files is not allowed.
337 (7) The stack depth limit is hard-wired to 100.
339 (8) All functions expecting a positive radius (e.g., LOG, KEI, etc.)
341 are passed the absolute value of their argument.
342
344 Regardless of the precision of the input data, GMT programs that create
345 gridded files will internally hold the grids in 4-byte floating point
346 arrays. This is done to conserve memory and futhermore most if not all
347 real data can be stored using 4-byte floating point values. Data with
348 higher precision (i.e., double precision values) will lose that preci‐
349 sion once GMT operates on the grid or writes out new grids. To limit
350 loss of precision when processing data you should always consider nor‐
351 malizing the data prior to processing.
352
354 By default GMT writes out grid as single precision floats in a COARDS-
355 complaint netCDF file format. However, GMT is able to produce grid
356 files in many other commonly used grid file formats and also facili‐
357 tates so called "packing" of grids, writing out floating point data as
358 2- or 4-byte integers. To specify the precision, scale and offset, the
359 user should add the suffix =id[/scale/offset[/nan]], where id is a two-
360 letter identifier of the grid type and precision, and scale and offset
361 are optional scale factor and offset to be applied to all grid values,
362 and nan is the value used to indicate missing data. When reading
363 grids, the format is generally automatically recognized. If not, the
364 same suffix can be added to input grid file names. See grdreformat(1)
365 and Section 4.17 of the GMT Technical Reference and Cookbook for more
366 information.
367 When reading a netCDF file that contains multiple grids, GMT will read,
369 by default, the first 2-dimensional grid that can find in that file. To
370 coax GMT into reading another multi-dimensional variable in the grid
371 file, append ?varname to the file name, where varname is the name of
372 the variable. Note that you may need to escape the special meaning of ?
373 in your shell program by putting a backslash in front of it, or by
374 placing the filename and suffix between quotes or double quotes. The
375 ?varname suffix can also be used for output grids to specify a variable
376 name different from the default: "z". See grdreformat(1) and Section
377 4.18 of the GMT Technical Reference and Cookbook for more information,
378 particularly on how to read splices of 3-, 4-, or 5-dimensional grids.
379
381 When the output grid type is netCDF, the coordinates will be labeled
382 "longitude", "latitude", or "time" based on the attributes of the input
383 data or grid (if any) or on the -f or -R options. For example, both
384 -f0x -f1t and -R90w/90e/0t/3t will result in a longitude/time grid.
385 When the x, y, or z coordinate is time, it will be stored in the grid
386 as relative time since epoch as specified by TIME_UNIT and TIME_EPOCH
387 in the .gmtdefaults file or on the command line. In addition, the unit
388 attribute of the time variable will indicate both this unit and epoch.
389
391 To take log10 of the average of 2 files, use
392 grdmath file1.grd file2.grd ADD 0.5 MUL LOG10 = file3.grd
394 Given the file ages.grd, which holds seafloor ages in m.y., use the
396 relation depth(in m) = 2500 + 350 * sqrt (age) to estimate normal
397 seafloor depths:
398 grdmath ages.grd SQRT 350 MUL 2500 ADD = depths.grd
400 To find the angle a (in degrees) of the largest principal stress from
402 the stress tensor given by the three files s_xx.grd s_yy.grd, and
403 s_xy.grd from the relation tan (2*a) = 2 * s_xy / (s_xx - s_yy), use
404 grdmath 2 s_xy.grd MUL s_xx.grd s_yy.grd SUB DIV ATAN2 2 DIV = direc‐
406 tion.grd
407 To calculate the fully normalized spherical harmonic of degree 8 and
409 order 4 on a 1 by 1 degree world map, using the real amplitude 0.4 and
410 the imaginary amplitude 1.1:
411 grdmath -R0/360/-90/90 -I1 8 4 YML 1.1 MUL EXCH 0.4 MUL ADD = harm.grd
413 To extract the locations of local maxima that exceed 100 mGal in the
415 file faa.grd:
416 grdmath faa.grd DUP EXTREMA 2 EQ MUL DUP 100 GT MUL 0 NAN = z.grd
418 grd2xyz z.grd -S > max.xyz
419
421 Abramowitz, M., and I. A. Stegun, 1964, Handbook of Mathematical Func‐
422 tions, Applied Mathematics Series, vol. 55, Dover, New York.
423 Holmes, S. A., and W. E. Featherstone, 2002, A unified approach to the
424 Clenshaw summation and the recursive computation of very high degree
425 and order normalised associated Legendre functions. Journal of
426 Geodesy, 76, 279-299.
427 Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
428 1992, Numerical Recipes, 2nd edition, Cambridge Univ., New York.
429 Spanier, J., and K. B. Oldman, 1987, An Atlas of Functions, Hemisphere
430 Publishing Corp.
431
433 GMT(1), gmtmath(1), grd2xyz(1), grdedit(1), grdinfo(1), xyz2grd(1)
434GMT 4.3.1 15 May 2008 GRDMATH(1)