surface
- 官方文档:
- 简介:
使用可调节张力连续曲率样条插值法对数据进行网格化
surface 从标准输入或表数据文件 table 中读取离散的数据点 (x,y,z)， 通过解如下微分方程，得到网格化的数据 z(x,y) 并生成网格文件
其中 t 是张力因子，取值范围为 0 到 1，\(\nabla\) 是拉普拉斯算子。 若 t = 0, 则退化为最小曲率张力样条。最小曲率张力样条可能会导致数据中 出现虚假的震荡，产生错误的区域极值 [见 Smith and Wessel, 1990]。 大于 0 的张力可削弱这些效应。对于势能场数据，建议 t ~ 0.25 。对于陡峭 的地形数据，张力则建议 t ~ 0.35。当 t = 1 时，则上述微分方程退化为 Laplace 方程，其解即为调和面，除控制点外，不会出现最大值或最小值。 建议用户事先使用 blockmean， blockmedian，或 blockmode 预处理数据，以避免空间混叠（即数据点之间距离太小）与消除冗余数据。 用户也可设置输出解的范围以保证结果总在预定的范围内。解微分方程时，采用 自然边界条件。对于经度范围是 360 度的地理数据，将采用周期边界条件。
语法
gmt surface [ table ] -Goutgrid -Iincrement -Rregion [ -Aaspect_ratio|m ] [ -Cconvergence_limit[%] ] [ -Jparameters ] [ -Dbreakline_file[+z[level]] ] [ -Lllower ] [ -Luupper ] [ -Mmax_radius ] [ -Nmax_iterations ] [ -Q[r] ] [ -Ssearch_radius[m|s] ] [ -T[b|i]tension_factor ] [ -V[level] ] [ -W[logfile] ] [ -Zover-relaxation_factor ] [ -aflags ] [ -bibinary ] [ -dinodata[+ccol] ] [ -eregexp ] [ -fflags ] [ -hheaders ] [ -iflags ] [ -qiflags ] [ -rreg ] [ -wflags ] [ -:[i|o] ] [ --PAR=value ]
必选选项
- table
一个或多个ASCII或二进制表数据。若不提供表数据，则会从标准输入中读取。
- -Ixinc[+e|n][/yinc[+e|n]]
指定X和Y方向的网格间隔
xinc 和 yinc 为 X 和 Y 方向的网格间隔。对于地理坐标，可以指定网格间隔单位 [默认单位为度]
+e 微调X和Y方向范围的最大值，使得其是网格间隔的整数倍（默认会微调 网格间隔以适应给定的数据范围）
+n 表明 xinc 和 yinc 不是网格间隔，而是X和Y方向的节点数。 此时会根据节点数、网格区域范围以及网格配准方式重新计算网格间隔。
注意：
若 yinc 设置为0，则表示其与 xinc 相同
若使用 -Rgrdfile 选项，则网格间隔和配准方式已经根据网格文件自动初始化， 此时依然可以使用 -I 和 -r 覆盖相应的值
- -Rxmin/xmax/ymin/ymax[+r][+uunit] (more …)
指定数据范围
可选选项
- -Aaspect_ratio|m
设置网格纵横比 aspect_ratio ，其中 dy = dx / aspect_ratio。 对于地理数据，可以使用 -Am 设置网格纵横比为纬度平均值的余弦。 不设置本项时，默认网格纵横比为 1，即各向同性网格。
- -Cconvergence_limit[%]
设置收敛极限，单位与 z 的单位相同。当任意格点的绝对变化量小于该极限时， 即退出迭代。另外可给定 z 值的 RMS 的一定百分比作为收敛极限，后面必须追加 %。 收敛极限默认设置为 RMS 的 10e-4。该值为最终输出的网格间距下的值，但对于 GMT 内部，用于处理的中间临时网格的收敛值为最终网格的收敛值除以网格间距的倍数。
- -Dbreakline[+z[level]]
指定一个 breakline, 其中数据格式为 x, y, z。在该 breakline 上的网格点，将 不会做内插。海岸线或者湖岸线通常是理想的 breakline。breakline 文件可以为多段 文件。如果文件中不存在 z，或者用户想要覆盖文件中的 z 值，则可在文件后追加 +zlevel，如果不指定 level 值，则默认为 0。
- -Lllower and -Luupper
对微分方程的解加上限制。 -Lllower 设置下限，其中 lower 可以是一个 网格文件的文件名或一个固定值。设置以后，输出的网格的值将均大于等于设置的值。 -Luupper 设置上限，同样 upper 也可以是一个网格文件的文件名或一个 固定值，输出处网格的值将均小于等于该值。如果设置为网格, 其大小、间隔和配准 应与输出网格一致。如果网格其中包含 NaN，则 NaN 位置处的输出将不受约束。 另外， lower 或 upper 除可为网格或固定值外，也可以设置为如下字符：
u 不设置约束 [默认]
d 针对 l 和 u 分别表示设置为输入数据的最小值和最大值
- -Mmax_radius
设置一个半径大小，在求解结束后，位于输出离散点半径内的格点数据将被保留，半径 以外的数据将设置为 NaN。可追加距离单位。 或可使用 -Mn_cellsc 形式以不用距离，而是格点个数来保留数据，其中 n_cells 为离散点周围的格点个数。-M0c 表示只有离散点所处的格点被保留。 -M1c 将保留处其所处格点外的最近一圈格点，即 3x3 的范围，依此类推。
- -Nmax_iterations
迭代次数, 默认为 500。迭代将会在达到收敛极限或迭代次数达到 max_iterations 时停止。这个 最大迭代次数仅针对最终输出的网格，对于 GMT 内部的临时网格，其迭代次数放大为 网格间隔相应的倍数。
- -Q[r]
输出本模块建议运行时使用的网格维度。该选项主要目的为建议一个最优的网格维度以 允许本模块可使用多次中间步骤以加快计算并得到更好的结果。当用户得到此建议维度 时，可使用 -R 和 -I 设置新的结果，并在后续使用 grdsample 或 grdcut 以得到用户原本希望的维度和范围。 另外，可追加 r 使 surface 完全按照 -R 指定的范围计算结果 [默认会自动 寻找稍大的区域，然后在输出时将区域裁剪为 -R 设置的范围。]
- -Ssearch_radius[m|s]
设置搜索半径 search_radius，单位与 x 和 y 的单位相同，默认值为 0，即不做搜索； 追加 m 或者 s 分别表示角分和角秒。该选项用于在第一次迭代之前进行初始化。除非网格点个数为素数， 不存在中间阶段，否则无需设置该选项。
- -T[b|i]tension_factor
设置张力因子，必须在 0 到 1 之间。 张力可用于内部求解（见上面的方程）和边界条件（张力可使接近边缘的解区域平缓）。 其中 i 用于设置内部求解张力，b 用于设置边界条件张力。如果不给出子选项， 则二者同时设置为同一个值。若同时设置为 0，即退化为最小曲率样条解，这也是默认情况。
- -W[logfile]
将收敛信息输出到为 logfile 文件，默认为 surface_log.txt
- -Zover-relaxation_factor
Over-relaxation 因子，该参数用于加速收敛，取值应位于 1 到 2 之间。当其为 1 时，将精确地迭代，并始终保证稳定收敛。较大的值会增大收敛中的增量变化， 因此可以更快地收敛，但由于增量变大，最终的结果可能不稳定。如果使用较大的 值，最好使用 -Vi 选项观察迭代过程。默认值为 1.4，可以快速收敛且总能 得到稳定的解。
- -V[level] (more …)
设置 verbose 等级 [w]
- -acol=name[,…] (more …)
设置非空间数据项与数据列之间的对应关系
- -bi[ncols][type][w][+l|b] (more …)
设置二进制输入数据的格式
- -dinodata (more …)
将输入数据中值为 nodata 的列替换为 NaN
- -e[~]“pattern” | -e[~]/regexp/[i] (more …)
筛选或剔除匹配指定模式的数据记录
- -f[i|o]colinfo (more …)
指定输入或输出列的数据类型
- -h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle] (more …)
跳过或生成指定数目的头段记录
- -icols[+l][+sscale][+ooffset][,…][,t[word]] (more …)
设置输入数据列及简单变换（0表示第一列，t 表示文本列）
- -qi[~]rows[+ccol][+a|f|s] (more …)
筛选输入的行或数据范围
- -r[g|p] (more …)
设置网格配置方式 [默认为网格线配准]
- -wy|a|w|d|h|m|s|cperiod[/phase][+ccol] (more …)
将输入坐标转换为循环坐标
- -:[i|o] (more …)
交换输入或输出中的第一和第二列
- -^ 或 -
显示简短的帮助信息，包括模块简介和基本语法信息（Windows下只能使用 -）
- -+ 或 +
显示帮助信息，包括模块简介、基本语法以及模块特有选项的说明
- -? 或无参数
显示完整的帮助信息，包括模块简介、基本语法以及所有选项的说明
- --PAR=value
临时修改GMT参数的值，可重复多次使用。参数列表见 配置参数
距离单位
示例
ASCII 数据 hawaii_5x5.xyg 为经 blockmean 处理后的数据，使用 surface 网格化为 5 分，设置 tension_factor = 0.25，convergence_limit = 0.1 mgal， 将结果输出到 havaii_grd.nc 并使用 -Vi 选项查看每个迭代过程
gmt surface hawaii_5x5.xyg -R198/208/18/25 -I5m -Ghawaii_grd.nc -T0.25 -C0.1 -Vi
注意事项
尽管 -R 选项控制了最终输出网格的范围，但在 GMT 内部会使用一个稍大一些的区域， 这种做法有利于最终结果的收敛。
格网化地理坐标数据：边界条件
surface 内部执行的有限差分算法基于笛卡尔坐标，因此设置了 -A 选项以适应 不同纬度下不同方向（即经度和纬度方向）的纵横比。当输出地理坐标网格的经度范围 为 0～360 范围时，将在经度方向实施周期边界条件。但是，在极点处，却不能使用该 边界条件，因为在极点处的有限差分解时无效的（在极点处，网格间隔为 0）。如果 用户在此情况下运行该模块，GMT 将给出警告，但不会阻止程序运行。 对于地理坐标，位于极点的值可能不一致，因此将被替换为极点附近的平均值，这可能 会导致在极点附近区域的结果存在问题。 鉴于上述问题，用户在处理全球数据时，应考虑使用球面网格化，例如：greenspline 或 sphinterpolate；或者使用 stereographic 投影数据并对投影后的数据 做 surface。
格网化地理坐标数据：设置间隔
对地理网格设置带有距离单位的网格间隔时，将会激活距离转换，即将间隔转换为等效 的度，在经纬度方向是不同的并且与椭球设置有关，但不论使用什么椭球，最终仍 使用球近似，即大圆距离（译注：椭球的大小是影响该距离的主要因素，球近似不是 主要因素）。对于纬度，GMT 将 y 方向的距离除以每维度所对应的距离，对于经度，则 除以平均纬度对应的距离。因此，转换后的度与用户输出的增量可能不完全一致。这 其中所涉及到的舍入是用户所不希望的，并会导致网格的维度为素数。-Q 选项可用来 处理这种情况，但间隔转换导致的无穷小数仍然难以有效解决。因此，不要使用长度 单位设置间隔是一种比较好的方法，设置为角度、角分、角秒可以得到更加干净的结果， 因为这样网格的维度是可控的并且可避免舍入误差。 另外，如果计算区域远离赤道，用户若担心经纬度方向不同的增量，则可选择对数据作 投影，这样就可以对投影坐标作网格化，间隔也可设置为米/千米等。 任意一种方法均可以避免类似 0.161697s 这种丑陋的间隔。如果用户所需结果却实为 以米/千米等距离为单位的网格，则投影方法是更适合的。
已知 Bug
在输出网格的任意节点存在多个输入数据，surface 会提示用户使用 blockmean， blockmedian，或 blockmode。 如果已经使用了上述模块，但仍存在警告，则意味着设置的网格间隔太小， 输出数据需要更多的有效数字。用户可在运行该模块前编辑 gmt.conf 中的 FORMAT_FLOAT_OUT 选项以增加有效数字个数，或使用单精度或者 双精度的二进制文件以保证其有效数字个数。
参考文献
Smith, W. H. F, and P. Wessel, 1990, Gridding with continuous curvature splines in tension, Geophysics, 55, 293-305, https://doi.org/10.1190/1.1442837.
相关模块
blockmean, blockmedian, blockmode, grdcut, grdsample, greenspline, nearneighbor, triangiulate, sphinterpolate