Question

贡献者

Caesar（自然资源部第一海洋研究所），周茂

---

官方文档

grdmath

简介

网格数据的逆波兰表示法（RPN）计算

grdmath 使用逆波兰表示法对网格文件或常量进行逐元素运算，最终结果写入到新的网格文件。若不给定网格文件，则必须使用

gmt grdmath [

operand

如果 operand 是一个文件，GMT 将试着以网格文件的形式读取该文件。如果非文件，将被解析为

Operator 为运算符名称； args 分别表示输入参数和输出参数个数； Returns 中的 A B C … 等字符表示输入参数。

Operator	Args	Returns
ABS	1 1	abs (A)
ACOS	1 1	acos (A)
ACOSH	1 1	acosh (A)
ACOT	1 1	acot (A)
ACSC	1 1	acsc (A)
ADD	2 1	A + B
AND	2 1	B if A == NaN, else A
ARC	2 1	Return arc(A,B) on [0 pi]
AREA	0 1	Area of each gridnode cell (in km^2 if geographic)
ASEC	1 1	asec (A)
ASIN	1 1	asin (A)
ASINH	1 1	asinh (A)
ATAN	1 1	atan (A)
ATAN2	2 1	atan2 (A, B)
ATANH	1 1	atanh (A)
BCDF	3 1	Binomial cumulative distribution function for p = A, n = B, and x = C
BPDF	3 1	Binomial probability density function for p = A, n = B, and x = C
BEI	1 1	Kelvin function bei (A)
BER	1 1	Kelvin function ber (A)
BITAND	2 1	A & B (bitwise AND operator)
BITLEFT	2 1	A << B (bitwise left-shift operator)
BITNOT	1 1	~A (bitwise NOT operator, i.e., return two’s complement)
BITOR	2 1	A | B (bitwise OR operator)
BITRIGHT	2 1	A >> B (bitwise right-shift operator)
BITTEST	2 1	1 if bit B of A is set, else 0 (bitwise TEST operator)
BITXOR	2 1	A ^ B (bitwise XOR operator)
BLEND	3 1	Blend A and B using weights in C (0-1 range) as A*C + B*(1-C)
CAZ	2 1	Cartesian azimuth from grid nodes to stack x,y (i.e., A, B)
CBAZ	2 1	Cartesian back-azimuth from grid nodes to stack x,y (i.e., A, B)
CDIST	2 1	Cartesian distance between grid nodes and stack x,y (i.e., A, B)
CDIST2	2 1	As CDIST but only to nodes that are != 0
CEIL	1 1	ceil (A) (smallest integer >= A)
CHICRIT	2 1	Chi-squared critical value for alpha = A and nu = B
CHICDF	2 1	Chi-squared cumulative distribution function for chi2 = A and nu = B
CHIPDF	2 1	Chi-squared probability density function for chi2 = A and nu = B
COMB	2 1	Combinations n_C_r, with n = A and r = B
CORRCOEFF	2 1	Correlation coefficient r(A, B)
COS	1 1	cos (A) (A in radians)
COSD	1 1	cos (A) (A in degrees)
COSH	1 1	cosh (A)
COT	1 1	cot (A) (A in radians)
COTD	1 1	cot (A) (A in degrees)
CSC	1 1	csc (A) (A in radians)
CSCD	1 1	csc (A) (A in degrees)
CUMSUM	2 1	Cumulative sum of each row (B=±1|3) or column (B=±2|4) in A. Sign of B gives direction of summation
CURV	1 1	Curvature of A (Laplacian)
D2DX2	1 1	d^2(A)/dx^2 2nd derivative
D2DY2	1 1	d^2(A)/dy^2 2nd derivative
D2DXY	1 1	d^2(A)/dxdy 2nd derivative
D2R	1 1	Converts Degrees to Radians
DDX	1 1	d(A)/dx Central 1st derivative
DAYNIGHT	3 1	1 where sun at (A, B) shines and 0 elsewhere, with C transition width
DDY	1 1	d(A)/dy Central 1st derivative
DEG2KM	1 1	Converts spherical degrees to kilometers
DENAN	2 1	Replace NaNs in A with values from B
DILOG	1 1	dilog (A)
DIV	2 1	A / B
DOT	2 1	2-D (Cartesian) or 3-D (geographic) dot products between nodes and stack (A, B) unit vector(s)
DUP	1 2	Places duplicate of A on the stack
ECDF	2 1	Exponential cumulative distribution function for x = A and lambda = B
ECRIT	2 1	Exponential distribution critical value for alpha = A and lambda = B
EPDF	2 1	Exponential probability density function for x = A and lambda = B
ERF	1 1	Error function erf (A)
ERFC	1 1	Complementary Error function erfc (A)
EQ	2 1	1 if A == B, else 0
ERFINV	1 1	Inverse error function of A
EXCH	2 2	Exchanges A and B on the stack
EXP	1 1	exp (A)
FACT	1 1	A! (A factorial)
EXTREMA	1 1	Local Extrema: +2/-2 is max/min, +1/-1 is saddle with max/min in x, 0 elsewhere
FCDF	3 1	F cumulative distribution function for F = A, nu1 = B, and nu2 = C
FCRIT	3 1	F distribution critical value for alpha = A, nu1 = B, and nu2 = C
FISHER	3 1	Fisher probability density function at nodes for center lon = A, lat = B, with kappa = C
FLIPLR	1 1	Reverse order of values in each row
FLIPUD	1 1	Reverse order of values in each column
FLOOR	1 1	floor (A) (greatest integer <= A)
FMOD	2 1	A % B (remainder after truncated division)
FPDF	3 1	F probability density function for F = A, nu1 = B, and nu2 = C
GE	2 1	1 if A >= B, else 0
GT	2 1	1 if A > B, else 0
HSV2RGB	3 3	Convert h,s,v triplets to r,g,b triplets, with h = A (0-360), s = B and v = C (both in 0-1)
HSV2RGB	3 3	Convert h,s,v triplets to r,g,b triplets, with h = A (0-360), s = B and v = C (0-1)
HSV2XYZ	3 3	Convert h,s,v triplets to x,t,z triplets, with h = A (0-360), s = B and v = C (0-1)
HYPOT	2 1	hypot (A, B) = sqrt (A*A + B*B)
I0	1 1	Modified Bessel function of A (1st kind, order 0)
I1	1 1	Modified Bessel function of A (1st kind, order 1)
IFELSE	3 1	B if A != 0, else C
IN	2 1	Modified Bessel function of A (1st kind, order B)
INRANGE	3 1	1 if B <= A <= C, else 0
INSIDE	1 1	1 when inside or on polygon(s) in A, else 0
INV	1 1	1 / A
ISFINITE	1 1	1 if A is finite, else 0
ISNAN	1 1	1 if A == NaN, else 0
J0	1 1	Bessel function of A (1st kind, order 0)
J1	1 1	Bessel function of A (1st kind, order 1)
JN	2 1	Bessel function of A (1st kind, order B)
K0	1 1	Modified Kelvin function of A (2nd kind, order 0)
K1	1 1	Modified Bessel function of A (2nd kind, order 1)
KEI	1 1	Kelvin function kei (A)
KER	1 1	Kelvin function ker (A)
KM2DEG	1 1	Converts kilometers to spherical degrees
KN	2 1	Modified Bessel function of A (2nd kind, order B)
KURT	1 1	Kurtosis of A
LAB2HSV	3 3	Convert l,a,b triplets to h,s,v triplets
LAB2RGB	3 3	Convert l,a,b triplets to r,g,b triplets
LAB2XYZ	3 3	Convert l,a,b triplets to x,y,z triplets
LCDF	1 1	Laplace cumulative distribution function for z = A
LCRIT	1 1	Laplace distribution critical value for alpha = A
LDIST	1 1	Compute minimum distance (in km if -fg) from lines in multi-segment ASCII file A
LDIST2	2 1	As LDIST, from lines in ASCII file B but only to nodes where A != 0
LDISTG	0 1	As LDIST, but operates on the GSHHG dataset (see -A, -D for options).
LE	2 1	1 if A <= B, else 0
LOG	1 1	log (A) (natural log)
LOG10	1 1	log10 (A) (base 10)
LOG1P	1 1	log (1+A) (accurate for small A)
LOG2	1 1	log2 (A) (base 2)
LMSSCL	1 1	LMS (Least Median of Squares) scale estimate (LMS STD) of A
LMSSCLW	2 1	Weighted LMS (Least Median of Squares) scale estimate (LMS STD) of A for weights in B
LOWER	1 1	The lowest (minimum) value of A
LPDF	1 1	Laplace probability density function for z = A
LRAND	2 1	Laplace random noise with mean A and std. deviation B
LT	2 1	1 if A < B, else 0
MAD	1 1	Median Absolute Deviation (L1 STD) of A
MAX	2 1	Maximum of A and B
MEAN	1 1	Mean value of A
MEANW	2 1	Weighted mean value of A for weights in B
MEDIAN	1 1	Median value of A
MEDIANW	2 1	Weighted median value of A for weights in B
MIN	2 1	Minimum of A and B
MOD	2 1	A mod B (remainder after floored division)
MODE	1 1	Mode value (Least Median of Squares) of A
MODEW	2 1	Weighted mode value (Least Median of Squares) of A for weights in B
MUL	2 1	A * B
NAN	2 1	NaN if A == B, else A
NEG	1 1	-A
NEQ	2 1	1 if A != B, else 0
NORM	1 1	Normalize (A) so max(A)-min(A) = 1
NOT	1 1	NaN if A == NaN, 1 if A == 0, else 0
NRAND	2 1	Normal, random values with mean A and std. deviation B
OR	2 1	NaN if B == NaN, else A
PCDF	2 1	Poisson cumulative distribution function for x = A and lambda = B
PDIST	1 1	Compute minimum distance (in km if -fg) from points in ASCII file A
PDIST2	2 1	As PDIST, from points in ASCII file B but only to nodes where A != 0
PERM	2 1	Permutations n_P_r, with n = A and r = B
PLM	3 1	Associated Legendre polynomial P(A) degree B order C
PLMg	3 1	Normalized associated Legendre polynomial P(A) degree B order C (geophysical convention)
POINT	1 2	Compute mean x and y from ASCII file A and place them on the stack
POP	1 0	Delete top element from the stack
POW	2 1	A ^ B
PPDF	2 1	Poisson distribution P(x,lambda), with x = A and lambda = B
PQUANT	2 1	The B’th Quantile (0-100%) of A
PQUANTW	3 1	The C’th weighted quantile (0-100%) of A for weights in B
PSI	1 1	Psi (or Digamma) of A
PV	3 1	Legendre function Pv(A) of degree v = real(B) + imag(C)
QV	3 1	Legendre function Qv(A) of degree v = real(B) + imag(C)
R2	2 1	R2 = A^2 + B^2
R2D	1 1	Convert Radians to Degrees
RAND	2 1	Uniform random values between A and B
RCDF	1 1	Rayleigh cumulative distribution function for z = A
RCRIT	1 1	Rayleigh distribution critical value for alpha = A
RGB2HSV	3 3	Convert r,g,b triplets to h,s,v triplets, with r = A, g = B, and b = C (all in 0-255 range)
RGB2LAB	3 3	Convert r,g,b triplets to l,a,b triplets, with r = A, g = B, and b = C (in 0-255 range)
RGB2XYZ	3 3	Convert r,g,b triplets to x,y,x triplets, with r = A, g = B, and b = C (in 0-255 range)
RINT	1 1	rint (A) (round to integral value nearest to A)
RMS	1 1	Root-mean-square of A
RMSW	1 1	Root-mean-square of A for weights in B
RPDF	1 1	Rayleigh probability density function for z = A
ROLL	2 0	Cyclicly shifts the top A stack items by an amount B
ROTX	2 1	Rotate A by the (constant) shift B in x-direction
ROTY	2 1	Rotate A by the (constant) shift B in y-direction
SDIST	2 1	Spherical (Great circle|geodesic) distance (in km) between nodes and stack (A, B)
SDIST2	2 1	As SDIST but only to nodes that are != 0
SAZ	2 1	Spherical azimuth from grid nodes to stack lon, lat (i.e., A, B)
SBAZ	2 1	Spherical back-azimuth from grid nodes to stack lon, lat (i.e., A, B)
SEC	1 1	sec (A) (A in radians)
SECD	1 1	sec (A) (A in degrees)
SIGN	1 1	sign (+1 or -1) of A
SIN	1 1	sin (A) (A in radians)
SINC	1 1	sinc (A) (sin (pi*A)/(pi*A))
SIND	1 1	sin (A) (A in degrees)
SINH	1 1	sinh (A)
SKEW	1 1	Skewness of A
SQR	1 1	A^2
SQRT	1 1	sqrt (A)
STD	1 1	Standard deviation of A
STDW	2 1	Weighted standard deviation of A for weights in B
STEP	1 1	Heaviside step function: H(A)
STEPX	1 1	Heaviside step function in x: H(x-A)
STEPY	1 1	Heaviside step function in y: H(y-A)
SUB	2 1	A - B
SUM	1 1	Sum of all values in A
TAN	1 1	tan (A) (A in radians)
TAND	1 1	tan (A) (A in degrees)
TANH	1 1	tanh (A)
TAPER	2 1	Unit weights cosine-tapered to zero within A and B of x and y grid margins
TCDF	2 1	Student’s t cumulative distribution function for t = A, and nu = B
TCRIT	2 1	Student’s t distribution critical value for alpha = A and nu = B
TN	2 1	Chebyshev polynomial Tn(-1<t<+1,n), with t = A, and n = B
TPDF	2 1	Student’s t probability density function for t = A, and nu = B
TRIM	3 1	Alpha-trim C to NaN if values fall in tails A and B (in percentage)
UPPER	1 1	The highest (maximum) value of A
VAR	1 1	Variance of A
VARW	2 1	Weighted variance of A for weights in B
VPDF	3 1	Von Mises density distribution P(x,mu,kappa), with x = A, mu = B, and kappa = C
WCDF	3 1	Weibull cumulative distribution function for x = A, scale = B, and shape = C
WCRIT	3 1	Weibull distribution critical value for alpha = A, scale = B, and shape = C
WPDF	3 1	Weibull density distribution P(x,scale,shape), with x = A, scale = B, and shape = C
WRAP	1 1	wrap A in radians onto [-pi,pi]
XOR	2 1	0 if A == NaN and B == NaN, NaN if B == NaN, else A
XYZ2HSV	3 3	Convert x,y,x triplets to h,s,v triplets
XYZ2LAB	3 3	Convert x,y,x triplets to l,a,b triplets
XYZ2RGB	3 3	Convert x,y,x triplets to r,g,b triplets
Y0	1 1	Bessel function of A (2nd kind, order 0)
Y1	1 1	Bessel function of A (2nd kind, order 1)
YLM	2 2	Re and Im orthonormalized spherical harmonics degree A order B
YLMg	2 2	Cos and Sin normalized spherical harmonics degree A order B (geophysical convention)
YN	2 1	Bessel function of A (2nd kind, order B)
ZCDF	1 1	Normal cumulative distribution function for z = A
ZPDF	1 1	Normal probability density function for z = A
ZCRIT	1 1	Normal distribution critical value for alpha = A

其他符号

以下符号都具有特殊意义：

PI	3.1415926…
E	2.7182818…
EULER	0.5772156…
PHI	1.6180339… (golden ratio)
EPS_F	1.192092896e-07 (single precision epsilon
XMIN	Minimum x value
XMAX	Maximum x value
XRANGE	Range of x values
XINC	x increment
NX	The number of x nodes
YMIN	Minimum y value
YMAX	Maximum y value
YRANGE	Range of y values
YINC	y increment
NY	The number of y nodes
X	Grid with x-coordinates
Y	Grid with y-coordinates
XNORM	Grid with normalized [-1 to +1] x-coordinates
YNORM	Grid with normalized [-1 to +1] y-coordinates
XCOL	Grid with column numbers 0, 1, …, NX-1
YROW	Grid with row numbers 0, 1, …, NY-1
NODE	Grid with node numbers 0, 1, …, (NX*NY)-1
NODEP	Grid with node numbers in presence of pad

上述符号均可以作为变量使用，当其为多个数时，逐元素操作。

运算符注意事项

1. 对于笛卡尔坐标，运算符 MEAN，MEDIAN，MODE，LMSSCL，MAD， PQUANT，RMS，STD 和 VAR 将返回基于给定网格的值。对于地理坐标，GMT 将施行球面加权运算，其中每个节点的值都由该点代表的地理区域加权得到。

2. SDIST 计算堆栈上的点 (lon,lat) 与所有网格点的球面距离，单位为 km。网格范围和点的坐标单位均为度。类似地，SAZ 和 SBAZ 分别用来球面方位角以及球面反方位角，单位也为度。如果设置了 -fg 或隐含该选项，LDIST 和 PDIST 以 km 为单位计算球面距离，否则返回笛卡尔距离。注 : 如果 PROJ_ELLIPSOID 为椭球，则计算结果为大地线（测地线）的长度。这可能导致计算比较慢，但可以通过 PROJ_GEODESIC 中设置使用其他算法以平衡精度和效率。

   LDISTG 和 LDIST 作用类似，但是 LDISTG 使用 GSHHG 数据运算。其中 GSHHG 数据来自

   用户可以将中间计算结果储存到一个变量中，并在后续计算中调用该变量。这在需要对某部分进行多次重复计算时可以提高效率和可读性。保存结果需要使用特殊的运算符 STO@label ， 其中 label 是变量的名称。调用该变量时，使用 [RCL]@label ， RCL 是可选的。使用后要清除该变量，可以使用 CLR@label ， STO 和 CLR 均不影响计算中的堆栈。

   GSHHG 信息

   GSHHG (The Global Self-consistent, Hierarchical, High-resolution Geography Database) 最初为 GSHHS，是一个海岸线数据库，主要有三个来源：World Vector Shorelines (WVS，不包括南极洲), CIA World Data Bank II (WDBII), 和 Atlas of the Cryosphere (AC，只包含南极洲)。除了南极洲，所有的 1 级多边形（海陆边界）都来自更准确的 WVS，所有的更高级别的多边形（2-4 级，湖泊边界、湖中岛边界和湖中岛中湖边界）来自 WDBII。南极洲的海岸线有两种：冰盖的边界和陆地的边界，可以使用 -A 选项选择。为了将 WVS 、WDBII 和 AC 数据转换为 GMT 可用的格式，GMT 已经进行了多种处理，包括：由线段组建多边形，去重，以及校正多边形之间的交叉等。每个多边形的面积已经被计算出来，因此，用户可以自己选择不绘制小于最小面积的多边形，还可以限制绘制的多边形的级别，见 -A 选项。绘制海岸线时，还可以使用不同的精细程度，其中低分辨率的海岸线是由 Douglas-Peucker 算法简化得到的。河流和边界的分类遵循 WDBII。详细细节见 GSHHG 。

   点位于多边形内/外

   为了确定点在多边形内，外或在边界上，GMT 会平衡数据类型以及多边形形状等因素来确定算法。对于笛卡尔坐标，GMT 使用 non-zero winding 算法，该算法非常快。对于地理坐标，如果多边形不包括两极点且多边形的经度范围不超过 360 度，同样使用该算法。否则，GMT 会采用 full spherical ray-shooting 方法。

   宏

   用户可以将特定的运算符组合保存为宏文件 gmtmath.macros 。文件中可以包含任意数量的宏， # 开头的行为注释。宏的格式为

   name = arg1 arg2 ... arg2[ : comment]
   

   其中，name 是宏名，当此运算符出现在命令中时，则将其简单替换为参数列表。宏不可以互相调用。下面给出一个宏例子：INCIRCLE 宏需要三个参数：半径，x0 和 y0，并将在圆内的点设置为 1，外部设置为 0

   INCIRCLE = CDIST EXCH DIV 1 LE : usage: r x y INCIRCLE to return 1 inside circl
   

   由于在宏中可能使用地理或时间常数，因此可使用 : 后加一个空格的形式作为注释的开端。

   示例

   计算网格点到北极的距离

   gmt grdmath -Rg -I1 0 90 SDIST = dist_to_NP.nc
   

   求两个文件的均值以及结果的 log10 函数值

   gmt grdmath file1.nc file2.nc ADD 0.5 MUL LOG10 = file3.nc
   

   给定含有海底地形年代（单位 m.y.）和深度（单位 m）网格文件 agrs.nc ，使用公式 depth(单位为 m) = 2500 + 350 * sqrt (age) 估计正常深度

   gmt grdmath ages.nc SQRT 350 MUL 2500 ADD = depths.nc
   

   从 s_xx.nc，s_yy.nc 和 s_xy.nc 三个文件包含的应力张量中，根据 tan (2*a) = 2 * s_xy / (s_xx - s_yy) 求出最大主应力的角度，单位为度

   gmt grdmath 2 s_xy.nc MUL s_xx.nc s_yy.nc SUB DIV ATAN 2 DIV = direction.nc
   

   计算 1 度分辨率的网格上完全正则化的 8 阶 4 次球谐函数，实数和虚数的振幅分别为 0.4 和 1.1

   gmt grdmath -R0/360/-90/90 -I1 8 4 YLM 1.1 MUL EXCH 0.4 MUL ADD = harm.nc
   

   提取文件:file:faa.nc 中超过 100 mgal 的局部最大值的位置

   gmt grdmath faa.nc DUP EXTREMA 2 EQ MUL DUP 100 GT MUL 0 NAN = z.nc
   gmt grd2xyz z.nc -s > max.xyz
   

   变量的使用: consider this radial wave where we store and recall the normalized radial arguments in radians

   gmt grdmath -R0/10/0/10 -I0.25 5 5 CDIST 2 MUL PI MUL 5 DIV STO@r COS @r SIN MUL = wave.nc
   

   创建一个保存 32 位浮点型 GeoTiff 文件

   gmt grdmath -Rd -I10 X Y MUL = lixo.tiff=gd:GTiff
   

   计算地理网格 data.grd 中网格点到 trace.txt 的距离，单位为 km

   gmt grdmath -Rdata.grd trace.txt LDIST = dist_from_line.grd
   

   -S 选项的使用：计算所有以 model_*.grd 为名的网格对应节点的标准差

   gmt grdmath model_*.grd -S STD = std_of_models.grd
   

   创建 0.5 度分辨率的 geotiff 网格文件，网格值为离海岸线的距离，单位为 km

   gmt grdmath -RNO,IS -Dc -I.5 LDISTG = distance.tif=gd:GTIFF
   

   参考文献

   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. Journal of Geodesy, 76, 279-299.

   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.

   相关模块

   gmtmath, grd2xyz, grdedit, grdinfo, xyz2grd