✍️ 朱邓达 • 📅 2025-05-12

grdgradient

官方文档:: grdgradient
简介:: 从网格文件中计算方向导数（光照强度）或梯度

语法

gmt grdgradient ingrid -Goutgrid [ -Aazim[/azim2] ] [ -D[a][c][n][o] ] [ -E[m][s][p]azim/elev[+aambient][+ddiffuse][+pspecular][+sshine] ] [ -N[e][t][amp][+aambient][+ssigma][+ooffset] ] [ -Qc|r|R[+ffile] ] [ -Rwest/east/south/north[/zmin/zmax][+r][+uunit] ] [ -Sslopfile ] [ -V[level] ] [ -fflags ] [ -nflags ] [ --PAR=value ]

描述

grdgradient 从网格文件中计算给定方向的方向导数(-A)，或梯度的方向(-D)和大小(-S)。

边界的梯度计算取决于网格的边界条件（详见中文文档 -n 选项或官方文档 -nflags ）。

必选选项

ingrid[=ID|?varname][+bband][+ddivisor][+ninvalid][+ooffset][+sscale]

输入网格名。通过追加 =ID 可指定网格格式 [默认为 =nf]；追加 ?varname 可指定 NetCDF 变量 [默认为 GMT 找到的第一个 2-D 网格]。对网格追加子选项可进行额外设置：

+b 选取一个波段 band（仅用于图片）[默认为 0]

+d 将网格除以一个数 divisor [默认为 1]

+n 将网格中的 invalid 替换为 NaN

+o 将网格中的值进行一定的偏移，即加一个数 offset [默认为 0]

+s 将网格中的值缩放，即乘以 scale [默认为 1]

注：offset 操作位于 scale 之前。

-Goutgrid[=ID][+ddivisor][+ninvalid][+ooffset|a][+sscale|a][:driver[dataType][+coptions]]

输出网格文件名。通过追加 =ID 可指定网格格式。对网格追加子选项可进行额外设置：

+d 将网格除以一个数 divisor [默认为 1]

+n 将网格中的 invalid 替换为 NaN

+o 将网格中的值进行一定的偏移，即加一个数 offset，或使用 a 自动对值进行调整以保证整数网格的精度 [默认为 0]

+s 将网格中的值缩放，即乘以 scale，或使用 a 自动对网格缩放以保证整数网格的精度 [默认为 1]

注：offset 操作位于 scale 之前；+sa 将会同时设置 +oa。如果需使用 GDAL 指定网格格式，则 ID 应设置为 gd，并指定 :driver 和可选的数据类型 dataType，以及 +coptions 选项传递给 GDAL 的 -co 选项。

可选选项

-Aazim[/azim2]: 给定光源方位角，计算方向导数。 azim 为方位角（度），在 x-y 平面从正北(+y)向着正东(+x)顺时针方向测量。实际计算的是方向导数的相反数， \(-\left( \frac{dz}{dx}\sin(a) + \frac{dz}{dy}\cos(a) \right)\)，其中 \(a\) 为方位角 azim 。这使得 \(z(x,y)\) 在沿着 \(a\) 方向下降时能输出正值，这符合在 x-y 平面上方从 \(a\) 方向对图像照明产生的阴影（见 grdimage 和 grdview ）。你也可以指定两个方位角， -Aazim[/azim2]，此时计算两个方向的方向导数，每个节点取方向导数幅值较大的方向导数。这对于照明具有两个方向的线性结构的数据是有用的，例如 -A0/270 表示从正北(top)和正西(left)两个方向照明。最后，如果 azim 是一个和 ingrid 相同区域、间隔、配准的网格文件，则计算方向导数时每个节点有各自的方位角。

-D[a][c][n][o]

不指定方向(-A)，而是计算每个节点的正梯度（上升）的方向。支持以下选项：

a - 计算负梯度（下降）的方向。
c - 使用传统的笛卡尔角度，从正东(+x)方向逆时针测量。默认从正北顺时针测量（ -A 中的 azim 一致）。
n - 对所有输出角度增加90度（例如获得表面的局部走向）。
o - 输出角度范围使用半圆定向（0~180度）而不是默认的全圆定向（0~360度）。例如全圆定向的30°和210°在半圆定向中均为30°，因为其位于同一延长线上。

-E[m][s][p]azim/elev[+aambient][+ddiffuse][+pspecular][+sshine]

计算朗伯辐射(Lambertian radiance)以用于 grdimage 和 grdview 。朗伯反射(Lambertian Reflection)又称理想散射，其假设一个理想的表面，它能反射所有照射到它的光，并且从各个方向看，表面亮度不变。azim 和 elev 分别为光照的方位角和仰角。以下选项可控制表面反射性质，中括号内为默认值，

+aambient - 环境光系数 [0.55]。
+ddiffuse - 漫反射系数 [0.6]。
+pspecular - 镜面反射系数 [0.4]。
+sshine - 镜面高光度 [10]。

s|p|m 指令可使用简化算法：

s - 简易朗伯算法（只要求 azim 和 elev ）。
p - 普克(Peucker)分段线性近似（简易但高效），此时 azim 和 elev 分别固定为315和45，忽略用户设置。
m - 与ESRI的hillshade算法类似，但更快。

-N[e][t][amp][+aambient][+ssigma][+ooffset]

归一化（默认不做归一化）。实际梯度 \(g\) 被偏移和缩放成 \(g_n\) ，幅值为 amp 。如果 amp 未给定，默认 amp = 1。如果偏移 offset 未给定，则设置为 \(g\) 的平均值。支持以下选项，其中 \(o\) 为 offset ， \(a\) 为 amp ，

-N - 根据 \(g_n = a\cdot \frac{g-o}{\max{|g-o|}}\) 进行归一化。
-Ne - 使用类似于拉普拉斯分布的累积分布函数 \(g_n = a\cdot\left( 1 - \exp{\left( \sqrt{2} \frac{g-o}{\sigma}\right)} \right)\) 进行归一化，其中 \(\sigma\) 为 sigma ，如果未给定则使用 \((g-o)\) 的L1范数。
-Nt - 使用类似于柯西分布的累计分布函数 \(g_n = \frac{2a}{\pi} \tan^{-1} \left(\frac{g-o}{\sigma}\right)\) 进行归一化，其中 \(\sigma\) 为 sigma ，如果未给定则使用 \((g-o)\) 的L2范数。

为了继续使用之前运行该命令计算的 \(o\) 和/或 \(\sigma\)，可不指定对应的具体数值，详见 -Q 。最后，可使用 +aambient 给归一化之后的梯度整体加上 ambient 。

-Qc|r|R[+ffile]: 控制 -N 的归一化。当需要对多个网格进行相同的归一化处理（即相同的偏移 offset 和标准差 sigma），用户必须使用 -N 传入这些值。然而，用户从网格中计算出这些值并不方便。可使用 -Qc 将计算的 offset 和 sigma 保存到一个统计文件中，如果这时不需要输出网格文件，可不指定 -G 。后续运行该命令时，可只使用 -Qr 来读取这些值。使用 -QR 则读取后删除该统计文件。也可以加上 +ffile 来指定读/写特定文件（默认为当前系统TMP目录下的grdgradient.stat文件）。

-Rxmin/xmax/ymin/ymax[+r][+uunit] (more …): 指定数据范围

-Sslopfile: 将梯度矢量的大小输出到网格文件 slopfile 。要求和 -D 使用，此时 -G 变为可选项。

-V[level] (more …): 设置 verbose 等级 [w]

-f[i|o]colinfo (more …): 指定输入或输出列的数据类型

-n[b|c|l|n][+a][+bBC][+c][+tthreshold] (more …): 设置网格文件的插值方式

-^ 或 -: 显示简短的帮助信息，包括模块简介和基本语法信息（Windows下只能使用 -）
-+ 或 +: 显示帮助信息，包括模块简介、基本语法以及模块特有选项的说明
-? 或无参数: 显示完整的帮助信息，包括模块简介、基本语法以及所有选项的说明
--PAR=value: 临时修改GMT参数的值，可重复多次使用。参数列表见配置参数

网格距离单位

详见官方文档关于 +u 的说明。

对于笛卡尔坐标系的网格，如果网格在水平方向的单位不为米，则使用 +uunit 指定输入单位，并在计算中自动转换为米。
对于地理坐标系的网格，使用 -fg 根据当前的椭球参数的“展平地球”将网格单位转换为米。

提示

如果你不知道怎么设置 -N 来为 grdimage 和 grdview 创建光照强度文件， -Ne0.6 是一个不错的尝试。
通常，255个色调对于可视化的目的来说已经足够了。可在输出文件名 outgrid 后面加上 =nb/a ，这可以节省75%的磁盘空间。 =nb/a 效果是输出网格文件格式为8-bit的整数，包含缩放因子。
如果你想为一个大数据集的子区域制作多个光照图，且要求在所有地图上保持一致的照明效果，使用 -N 选项并为每个地图提供相同的 sigma 和 offset 。一个不错的尝试是 offset = 0， sigma 通过 grdinfo -L2 或 -L1 应用于未归一化的梯度网格得到。
如果您只需要网格的x或y导数，请使用 grdmath 。

分块处理

针对超大规模数据集（或超大尺寸绘图），您可能需要将任务分解为多个区块处理。此时需特别注意：每个区块的光照强度归一化处理应保持一致。默认情况下，程序会为每个区块单独重新计算偏移量( offset )和标准差( sigma )。这会导致同一大型网格中的不同区块计算出不同的偏移量和标准差，最终使得相同方向的光照强度值在整张地图上呈现不一致。此类问题可能导致区块拼接处出现肉眼可见的色差/明暗变化。有以下两种解决方法：

手动指定参数。使用 -N 强制指定相同的 offset 和 sigma ，这适用于已知标准化参数的情况。
自动传递参数（推荐工作流）。首先对完整网格（或计算机能处理的最大网格分区）运行 grdgradient -Qc ，该操作会生成包含全局 offset 和 sigma 的隐藏统计文件。然后在后续各区块计算时，在 -N 后添加无参数的 +o 和/或 -s ，同时指定 -Qr 以从统计文件中读取参数。在处理最后一个区块时可改用 -QR 使程序读取统计文件后自动删除。

环境光

环境光 ambient 偏移用于使光照强度整体变暗或变亮。这个设置通常用于通过从所有光照强度中减去一个常数来使整个图像变暗。例如，如果你使用 +a-0.5，那么你从所有光照强度中减去0.5，使它们更负，从而使图像变暗。

特效

如果用户希望只突出显示网格的一个区域，我们鼓励创造性地使用 grdmath 中的大量运算符和 grdmask 和 grdclip 中的网格裁剪方案。通过 grdmath 中的各种网格代数和布尔运算符，可以将多个独立的光照强度网格组合成一个网格。天空才是极限。

示例

注: 下面是该模块的一些有效语法示例。使用远程文件（文件名以 @ 开头）的示例可以剪切并粘贴到您的终端中进行测试。其他需要输入文件的命令只是常见的使用类型的虚拟示例，但不能按编写的原样运行。

在远程文件 @earth_relief_05m 中使用指数归一化梯度在[-0.6,0.6]范围内模拟北、西方向的光源，制作一个文件来照亮部分数据:

gmt grdgradient @earth_relief_05m -R0/20/0/20 -A0/270 -Ggradients.nc -Ne0.6 -V

基于地形文件 topo.nc 确定海底构造(seafloor fabric)的走向:

gmt grdgradient topo.nc -Dno -Gazimuths.nc -V

保存对 topo.nc 的光照强度值进行归一化的 offset 和 sigma

gmt grdgradient topo.nc -A30 -Nt0.6 -Qc -V

其中没有指定 -G ，不会有网格文件输出，仅创建隐藏的统计文件。

使用先前确定的 offset 和 sigma 对 tile_3.nc 的光照强度进行归一化:

gmt grdgradient tile_3.nc -A30 -Nt0.6+o+s -Qr -V -Gtile_3_int.nc

参考文献

Horn, B.K.P., Hill-Shading and the Reflectance Map, Proceedings of the IEEE, Vol. 69, No. 1, January 1981, pp. 14-47. (https://www.researchgate.net/publication/2996084_Hill_shading_and_the_reflectance_map)