grdmath will perform operations like add, subtract, multiply, and
hundreds of other operands on one or more grid files or constants using
Reverse Polish Notation (RPN)
syntax. Arbitrarily complicated expressions may therefore be evaluated; the
final result is written to an output grid file. Grid operations are
element-by-element, not matrix manipulations. Some operators only
require one operand (see below). If no grid files are used in the
expression then options -R, -I must be set (and optionally
-rreg). The expression =outgrid can occur as many times as
the depth of the stack allows in order to save intermediate results.
Complicated or frequently occurring expressions may be coded as a macro
for future use or stored and recalled via named memory locations.
Hewlett-Packard made lots of calculators (left) using Reverse Polish Notation,
which is a post-fix system for mathematical notion originally developed by
Jan_Łukasiewicz (right).
Here, operands are entered first followed by an operator, e.g., “3 5 +” instead
of “3 + 5 =” (photo courtesy of John W. Robbins).
If operand can be opened as a file it will be read as a grid file.
If not a file, it is interpreted as a numerical constant or a
special symbol (see below).
outgrid
The name of a 2-D grid file that will hold the final result. (See
Grid File Formats).
Features with an area smaller than min_area in km2 or of
hierarchical level that is lower than min_level or higher than
max_level will not be plotted [Default is 0/0/4 (all features)].
Level 2 (lakes) contains regular lakes and wide river bodies which
we normally include as lakes. Several modifiers provide further control:
+a - Control special aspects of the Antarctica coastline. Append one
of g|s|s|S:
g - Selects the Antarctica ice grounding line as the coastline.
i - Selects the ice shelf boundary as the coastline for Antarctica [Default].
s - Skip all GSHHG features below 60S (For users who wish to utilize their
own Antarctica (with islands) coastline).
S - Like s but skip instead all features north of 60S.
+l - Only regular lakes and exclude river-lakes.
+p - Append percent to exclude polygons whose percentage area of the
corresponding full-resolution feature is less than percent.
+r - Only select river-lakes and exclude regular lakes.
See GSHHG Information below for more details. (-A is only relevant to the LDISTG operator)
-C[cpt]
Retain the grid’s default CPT (if it has one), or alternatively replace it with a
new default cpt [Default removes any default CPT from the output grid].
-Dresolution[+f]
Selects the resolution of the data set to use with the operator LDISTG
((f)ull, (h)igh, (i)ntermediate, (l)ow, and (c)rude). The
resolution drops off by 80% between data sets [Default is l].
Append +f to automatically select a lower resolution should the one
requested not be available [abort if not found].
-Ix_inc[+e|n][/y_inc[+e|n]]
Set the grid spacing as x_inc [and optionally y_inc].
Geographical (degrees) coordinates: Optionally, append an increment unit.
Choose among:
d - Indicate arc degrees
m - Indicate arc minutes
s - Indicate arc seconds
If one of e (meter), f (foot), k (km), M (mile), n (nautical
mile) or u (US survey foot), the the increment will be
converted to the equivalent degrees longitude at the middle latitude
of the region (the conversion depends on PROJ_ELLIPSOID). If
y_inc is not given or given but set to 0 it will be reset equal to x_inc;
otherwise it will be converted to degrees latitude.
All coordinates: The following modifiers are supported:
+e - Slightly adjust the max x (east) or y (north) to fit
exactly the given increment if needed [Default is to slightly adjust
the increment to fit the given domain].
+n - Define the number of nodes rather than the increment, in
which case the increment is recalculated from the number of nodes,
the registration (see GMT File Formats), and the domain.
Note: If -Rgrdfile is used then the grid spacing and
the registration have already been initialized; use -I and -R
to override these values.
-M
By default any derivatives calculated are in z_units/x(or
y)_units. However, the user may choose this option to convert dx,dy
in degrees of longitude,latitude into meters using a flat Earth
approximation, so that gradients are in z_units/meter.
-N
Turn off strict domain match checking when multiple grids are
manipulated [Default will insist that each grid domain is within
\(10^{-4}\) times the grid spacing of the domain of the first grid listed].
Reduce (i.e., collapse) the entire stack to a single grid by applying the
next operator to all co-registered nodes across the entire stack. You
must specify -Safter listing all of your grids. Note: You can only
follow -S with a reducing operator, i.e., from the list ADD, AND, MAD,
LMSSCL, MAX, MEAN, MEDIAN, MIN, MODE, MUL, RMS,
STD, SUB, VAR or XOR.
Select native binary format for primary table input. The binary input option
only applies to the data files needed by operators LDIST,
PDIST, and INSIDE.
For Cartesian grids the operators MEAN, MEDIAN, MODE,
LMSSCL, MAD, PQUANT, RMS, STD, and VAR return the
expected value from the given matrix. However, for geographic grids
we perform a spherically weighted calculation where each node value
is weighted by the geographic area represented by that node.
The operator SDIST calculates spherical distances in km between the
(lon, lat) point on the stack and all node positions in the grid. The
grid domain and the (lon, lat) point are expected to be in degrees.
Similarly, the SAZ and SBAZ operators calculate spherical
azimuth and back-azimuths in degrees, respectively. The operators
LDIST and PDIST compute spherical distances in km if -fg is
set or implied, else they return Cartesian distances. Note: If the current
PROJ_ELLIPSOID is ellipsoidal then
geodesics are used in calculations of distances, which can be slow.
You can trade speed with accuracy by changing the algorithm used to
compute the geodesic (see PROJ_GEODESIC).
The operator LDISTG is a version of LDIST that operates on the
GSHHG data. Instead of reading an ASCII file, it directly accesses one of
the GSHHG data sets as determined by the -D and -A options.
The operator POINT reads a ASCII table, computes the mean x and mean
y values and places these on the stack. If geographic data then we use
the mean 3-D vector to determine the mean location.
The operator PLM calculates the associated Legendre polynomial
of degree L and order M (\(0 \leq M \leq L)\), and its argument is the sine of
the latitude. PLM is not normalized and includes the Condon-Shortley
phase \((-1)^M\). PLMg is normalized in the way that is most commonly
used in geophysics. The Condon-Shortley phase can be added by using -M as argument.
PLM will overflow at higher degrees, whereas PLMg is stable
until ultra high degrees (at least 3000).
The operators YLM and YLMg calculate normalized spherical
harmonics for degree L and order M (\(0 \leq M \leq L)\) for all positions in
the grid, which is assumed to be in degrees. YLM and YLMg return
two grids, the real (cosine) and imaginary (sine) component of the
complex spherical harmonic. Use the POP operator (and EXCH) to
get rid of one of them, or save both by giving two consecutive = file.nc calls.
The orthonormalized complex harmonics YLM are most commonly used in
physics and seismology. The square of YLM integrates to 1 over a
sphere. In geophysics, YLMg is normalized to produce unit power when
averaging the cosine and sine terms (separately!) over a sphere (i.e.,
their squares each integrate to \(4 \pi\)). The Condon-Shortley phase \((-1)^M\)
is not included in YLM or YLMg, but it can be added by using -M
as argument.
All the derivatives are based on central finite differences, with
natural boundary conditions, and are Cartesian derivatives.
Files that have the same names as some operators, e.g., ADD,
SIGN, =, etc. should be identified by prepending the current
directory (i.e., ./LOG).
Piping of files is not allowed.
The stack depth limit is hard-wired to 100.
All functions expecting a positive radius (e.g., LOG, KEI,
etc.) are passed the absolute value of their argument. (9) The bitwise
operators (BITAND, BITLEFT, BITNOT, BITOR, BITRIGHT,
BITTEST, and BITXOR) convert a grid’s single precision values to
unsigned 32-bit ints to perform the bitwise operations. Consequently,
the largest whole integer value that can be stored in a float grid is
2^24 or 16,777,216. Any higher result will be masked to fit in the lower
24 bits. Thus, bit operations are effectively limited to 24 bit. All
bitwise operators return NaN if given NaN arguments or bit-settings <= 0.
When OpenMP support is compiled in, a few operators will take advantage
of the ability to spread the load onto several cores. At present, the
list of such operators is: LDIST, LDIST2, PDIST, PDIST2,
SAZ, SBAZ, SDIST, YLM, and grd_YLMg.
Operators DEG2KM and KM2DEG are only exact when a spherical Earth
is selected with PROJ_ELLIPSOID.
Operator DOT normalizes 2-D vectors before the dot-product takes place.
For 3-D vector they are all unit vectors to begin with.
The color-triplet conversion functions (RGB2HSV, etc.) includes not
only r,g,b and h,s,v triplet conversions, but also l,a,b (CIE L a b ) and
sRGB (x,y,z) conversions between all four color spaces.
The DAYNIGHT operator returns a grid with ones on the side facing the given
sun location at (A,B). If the transition width (C) is zero then we get
either 1 or 0, but if C is nonzero then we approximate the step function
using an atan-approximation instead. Thus, the values are never exactly
0 or 1, but close, and the smaller C the closer we get.
The VPDF operator expects angles in degrees.
The CUMSUM operator normally resets the accumulated sums at the end of a
row or column. Use ±3 or ±4 to have the accumulated sums continue with the
start of the next row or column.
Regardless of the precision of the input data, GMT programs that create
grid files will internally hold the grids in 4-byte floating point
arrays. This is done to conserve memory and furthermore most if not all
real data can be stored using 4-byte floating point values. Data with
higher precision (i.e., double precision values) will lose that
precision once GMT operates on the grid or writes out new grids. To
limit loss of precision when processing data you should always consider
normalizing the data prior to processing.
When the output grid type is netCDF, the coordinates will be labeled
“longitude”, “latitude”, or “time” based on the attributes of the input
data or grid (if any) or on the -f or -R options. For example,
both -f0x-f1t and -R90w/90e/0t/3t will result in a
longitude/time grid. When the x, y, or z coordinate is time, it will be
stored in the grid as relative time since epoch as specified by
TIME_UNIT and TIME_EPOCH in the
gmt.conf file or on the
command line. In addition, the unit attribute of the time variable
will indicate both this unit and epoch.
You may store intermediate calculations to a named variable that you may
recall and place on the stack at a later time. This is useful if you
need access to a computed quantity many times in your expression as it
will shorten the overall expression and improve readability. To save a
result you use the special operator STO@label, where label
is the name you choose to give the quantity. To recall the stored result
to the stack at a later time, use [RCL]@label, i.e., RCL
is optional. To clear memory you may use CLR@label. Note that
both STO and CLR leave the stack unchanged.
The coastline database is GSHHG (formerly GSHHS) which is compiled from
three sources: World Vector Shorelines (WVS, not including Antarctica),
CIA World Data Bank II (WDBII), and Atlas of the Cryosphere (AC, for Antarctica only).
Apart from Antarctica, all level-1 polygons (ocean-land boundary) are derived from
the more accurate WVS while all higher level polygons (level 2-4,
representing land/lake, lake/island-in-lake, and
island-in-lake/lake-in-island-in-lake boundaries) are taken from WDBII.
The Antarctica coastlines come in two flavors: ice-front or grounding line,
selectable via the -A option.
Much processing has taken place to convert WVS, WDBII, and AC data into
usable form for GMT: assembling closed polygons from line segments,
checking for duplicates, and correcting for crossings between polygons.
The area of each polygon has been determined so that the user may choose
not to draw features smaller than a minimum area (see -A); one may
also limit the highest hierarchical level of polygons to be included (4
is the maximum). The 4 lower-resolution databases were derived from the
full resolution database using the Douglas-Peucker line-simplification
algorithm. The classification of rivers and borders follow that of the
WDBII. See The Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) for further details.
To determine if a point is inside, outside, or exactly on the boundary of a polygon
we need to balance the complexity (and execution time) of the algorithm with the
type of data and shape of the polygons. For any Cartesian data we use a non-zero
winding algorithm, which is quite fast. For geographic data we will also use this algorithm
as long as (1) the polygons do not include a geographic pole, and (2) the longitude extent
of the polygons is less than 360. If this is the situation we also carefully adjust
the test point longitude for any 360 degree offsets, if appropriate. Otherwise,
we employ a full spherical ray-shooting method to determine a points status.
Users may save their favorite operator combinations as macros via the
file grdmath.macros in their current or user (~/.gmt) directory. The file may contain
any number of macros (one per record); comment lines starting with # are
skipped. The format for the macros is name = arg1 arg2 … arg2
: comment where name is how the macro will be used. When this
operator appears on the command line we simply replace it with the
listed argument list. No macro may call another macro. As an example,
the following macro expects three arguments (radius x0 y0) and sets the
nodes that are inside the given Cartesian circle to 1 and those outside to 0:
INCIRCLE = CDIST EXCH DIV 1 LE : usage: r x y INCIRCLE to return 1
inside circle
Marine geophysicist often need to evaluate predicted seafloor depth from an age grid.
One such model is the classic Parsons and Sclater [1977] curve. It may be written
A good cross-over age for these curves is 26.2682 Myr. A macro for this system is a bit awkward due to the split but can be written
PS77 = STO@T POP 6400 RCL@T 62.8 DIV NEG EXP 3200 MUL SUB RCL@T 26.2682 GT MUL 2500 350 RCL@T SQRT MUL ADD RCL@T 26.2682 LE MUL ADD : usage: age PS77 returns depth
i.e., we evaluate both expressions and multiply them by 1 or 0 depending where they apply, and then add them. With this macro
installed you can compute predicted depths in the Pacific northwest via:
Note: Below are some examples of valid syntax for this module.
The examples that use remote files (file names starting with @)
can be cut and pasted into your terminal for testing.
Other commands requiring input files are just dummy examples of the types
of uses that are common but cannot be run verbatim as written.
To compute all distances to north pole, try:
gmtgrdmath-Rg-I1090SDIST=dist_to_NP.nc
To take \(\log_{10}\) of the average of 2 files, use:
gmtgrdmathfile1.ncfile2.ncADD0.5MULLOG10=file3.nc
Given the file ages.nc, which holds seafloor ages in m.y., use the
relation depth(in m) = 2500 + 350 * sqrt (age) to estimate normal seafloor depths, try:
gmtgrdmathages.ncSQRT350MUL2500ADD=depths.nc
To find the angle a (in degrees) of the largest principal stress from
the stress tensor given by the three files s_xx.nc s_yy.nc, and
s_xy.nc from the relation tan (2*a) = 2 * s_xy / (s_xx - s_yy), use:
To calculate the fully normalized spherical harmonic of degree 8 and
order 4 on a 1 by 1 degree world map, using the real amplitude 0.4 and
the imaginary amplitude 1.1, use:
Abramowitz, M., and I. A. Stegun, 1964, Handbook of Mathematical
Functions, Applied Mathematics Series, vol. 55, Dover, New York.
Holmes, S. A., and W. E. Featherstone, 2002, A unified approach to the
Clenshaw summation and the recursive computation of very high degree and
order normalized associated Legendre functions. J. of Geodesy,
76, 279-299.
B. Parsons and J. G. Sclater, 1977, An analysis of the variation of ocean
floor bathymetry and heat flow with age, J. Geophys. Res., 82, 803-827.
Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
1992, Numerical Recipes, 2nd edition, Cambridge Univ., New York.
Spanier, J., and K. B. Oldman, 1987, An Atlas of Functions, Hemisphere
Publishing Corp.
More on Reverse Polish Notation
Reverse Polish Notation (RPN) or postfix notation is widely used in computer science because of the simplicity of implementing a
stack-based computation system. The PostScipt language we used to make GMT graphics is postfix, for instance.
To learn more about RPN or postfix, watch this YouTube video explanation:
