img2grd
- 贡献者:
- 官方文档:
- 简介:
从墨卡托 img 格式文件中提取网格数据
img2grd 读取 img 格式的文件，提取数据并写入到网格文件中。-M 选项用来设置在提取的网格文件中保留球墨卡托投影，如果不使用该选项，则取消保留球墨卡托投影，结果为地理坐标网格。
img 格式是一种二进制网格文件，数据为球墨卡托投影后的结果。球墨卡托投影原理决定其纬度间隔在转换为地理坐标时不是等间距的，因此导致了下述两个坐标系统下的区域范围轻微的差异。球墨卡托投影原理还决定了 img 文件在纬度方向并不是全球覆盖的，在极区有空白，即 -D 选项设置的最大纬度 80.738 ；在经度方向则是全覆盖的。虽然该模块可以处理任何 img 格式文件，但目前较少的产品或模型以 img 格式发布。Sandwell 和 Smith 基于卫星测高数据创建的海洋重力场（包括重力异常、垂线偏差和垂直重力梯度）和海底地形模型采用了该格式发布。由于测高卫星轨道倾角的限制，在不同时期，上述模型的纬度范围不同，早期最大纬度为 72.006，当前最大纬度范围为 80.738 。
Sandwell 和 Smith 模型的大小和文件名有一定规律，用于在 GMT 中自动设置参数并处理，文件名在不同时期的命名方式不同，但规律类似，例如：world_grav.img.7.2 和 grav_28.1.img
grav 表示该文件为自由空气重力异常
img 表示文件格式为 img
7.2 表示该模型版本号为 7，分辨率为 2 分；28.1 表示版本号为 28，分辨率为 1 分
语法
gmt img2grd imgfile -Ggrdfile -Rregion [ -D[ minlat/maxlat ] ] [ -E ] [ -Iinc ] [ -F ] [ -M ] [ -Nnavg ] [ -S[scale] ] [ -Ttype ] [ -V[level] ] [ -Wmaxlon ] [ -nflags ] [ --PAR=value ]
必选选项
- imgfile
输入墨卡托 img 格式文件名，例如： Sandwell 和 Smith 创建的海洋重力场或海底地形文件。
- -Ggrdfile
输出网格文件名。
- -Rxmin/xmax/ymin/ymax[+r][+uunit] (more …)
指定数据范围
可选选项
- -D[minlat/maxlat]
设置纬度范围为 -80.738/+80.738。附加 minlat/maxlat 可以指定输入的 img 文件的纬度范围 [ 默认为 -72.006/72.006 ]。该选项通常不必使用，因为本模块可以自动确定输入文件的纬度范围。
- -E
当输出网格为地理坐标时，即未设置 -M 选项时，使用 -E 选项可以强制地理坐标格网符合纬度范围并强制经纬度间隔相等。默认情况下，本模块首先从输入文件中提取 -R 指定范围内的网格，网格的坐标形式和原始文件相同，即球墨卡托投影。因此，当输出网格为地理坐标时，需将球墨卡托投影转换为地理坐标。在转换过程中，通常会稍微超出指定的纬度范围，且在纬度方向的间隔和经度方向不匹配。但是，这种强制经纬度范围等间隔的重采样会导致内插误差。因此，只有强制输出网格为指定范围内，且经纬度间隔相同时才建议使用。在这种情况下，-R 设置的区域必须为经纬度间隔的整数倍。
- -F
平移 x 和 y 墨卡托坐标，使其位于左下角，即 (0,0)
- -I
指定输入的 img 文件纬度间隔 inc [ 默认为 2 ]。添加 m [ 默认 ] 或 s 来指定单位。该选项通常不必使用，因为本模块可以根据输入 img 文件大小自动判断。
- -M
输出球墨卡托网格 [ 默认为地理 lon/lat 网格 ]。输入文件的的球墨卡托投影被保留，因此 -R 指定的区域被轻微地调整，调整后的区域对应于像素的边界。设置输出网格文件的头部信息，x 和 y 轴方向的距离相对于使用 -Jm1 投影下的投影中心。使用 -F 选项可以将 x 和 y 轴的距离设置为相对于输出范围 -R 的西南角，距离单位为用户默认的单位。
- -Nnavg
对输入的墨卡托像素在 navg x navg*的范围内平均，并创建平均后的像素点。 如果和 *-T3** 同时使用，输出网格值为平均约束值，范围为 0 到 1。如果和 -T2 同时使用，输出网格值为像素平均值或 NaN，平均约束值 > 0.5 时，为像素平均值，平均约束值 < 0.5 时，为 NaN。 navg 的单位为像素 [ 默认值为 1，表示不进行平均]。
- -S[scale]
对输出的网格文件进行 scale 比例缩放，以保证结果单位正确[ 默认是 1.0 ]。对于目前的 img 文件：海底地形，单位为米，因此设置为 -S1 得到米 ；自由空气重力异常，单位为 mGal*10 ，因此设置为 -S0.1 得到 mGal ；垂线偏差，单位为 micro-radians*10 ， 因此设置为 -S0.1 得到 micro-radians ；垂直重力梯度，单位为 Eotvos*10 ， 因此设置为 -S0.1 得到 Eotvos 或者设置为 -S0.01 得到 mGal/km。如果不给定 scale，则将通过文件名来自动确定。
- -Ttype
type 处理约束信息的编码（译注：约束信息可理解为数据可靠性的一种衡量标准）。 type = 0 表明没有约束信息编码到 img 文件（1995 年前的重力场模型）中，并获取所有数据。 type > 0 表明有约束信息编码到 img 文件 (1995 以后以及当前版本的 img 文件)。 -T1 获取所有数据； -T2 获得约束点的数据值，内插点的数值设置为 NaN ； -T3 约束点的数据值设置为 1，内插点设置为 0 [ 默认值为 1]。
- -V[level] (more …)
设置 verbose 等级 [w]
- -Wmaxlon
指定输入 img 文件的最大经度 maxlon。从 1995 年开始，Sandwell 和 Smith 的模型的 maxlon = 360.0 ，但是某些早期的模型 maxlon = 390.0 [ 默认为 360.0 ]。
- -n[b|c|l|n][+a][+bBC][+c][+tthreshold] (more …)
设置网格文件的插值方式
- -^ 或 -
显示简短的帮助信息，包括模块简介和基本语法信息（Windows下只能使用 -）
- -+ 或 +
显示帮助信息，包括模块简介、基本语法以及模块特有选项的说明
- -? 或无参数
显示完整的帮助信息，包括模块简介、基本语法以及所有选项的说明
- --PAR=value
临时修改GMT参数的值，可重复多次使用。参数列表见 配置参数
地理坐标示例
如果需要输出地理坐标网格，则不能使用 -M 选项。从 grav_28.1.img 文件中提取数据并转换为地理坐标:
gmt img2grd grav_28.1.img -Gmerc_grav.nc -R-40/40/-70/-30 -V
由于在 img 文件的墨卡托单位中，纬度间隔是相等的，转换为地理坐标时导致纬度间隔不相等，因此最终提取的网格并不能严格匹配 -R 选项指定的纬度范围。使用 -E 选项可以精确的匹配 -R 指定的范围，并且强制经度和纬度的间隔相同:
gmt img2grd grav_28.1.img -Gmerc_grav.nc -R-40/40/-70/-30 -E -V
墨卡托示例
数据处理示例
由于 img 格式为墨卡托投影，因此如果最终数据处理结果为墨卡托地图，不建议使用该模块提取地理坐标，而是直接使用墨卡托投影坐标处理。如果强行从 img 提取到地理坐标再转换为墨卡托投影会导致丢失短波信息。因此，最好使用 -M 选项并与所需的墨卡托投影设置一定的线性比例缩放（见 GMT 例 29）。从 world_grav.img.7.2 中提取 -R-40/40/-70/-30 范围内的数据:
gmt img2grd -M world_grav.img.7.2 -Gmerc_grav.nc -R-40/40/-70/-30 -V
注意到 -V 选项告诉用户区域范围被调整为 -R-40/40/-70.0004681551/-29.9945810754 。如果在脚本编写中需要使用该范围，可以使用 grdinfo -Ii。此外，用户也可以从 grdinfo 的信息中看到，网格文件头部信息显示区域范围调整为 -R-40/40/-99.4333333333/-31.4666666667。该纬度范围为使用 -Jm1 和 -R-40/40/-70.0004681551/-29.9945810754 经过球墨卡托投影的结果。使用 ship.lonlatgrav 对 merc_grav.nc 进行采样的代码为:
gmt set PROJ_ELLIPSOID Sphere
gmt mapproject -R-40/40/-70.0004681551/-29.9945810754 -Jm1i -C ship.lonlatgrav | \
gmt grdtrack -Gmerc_grav.nc | \
gmt mapproject -R-40/40/-70.0004681551/-29.9945810754 -Jm1i -I -C > ship.lonlatgravsat
如上述代码，建议对数据投影和重投影，而不是对 img 文件进行类似的操作。因为对数据处理只进行一次插值（即 grdtrack）, 对 img 文件操作会进行两次插值（分别在转换和重采样中）。
如果想从上面的墨卡托网格中提取地理坐标网格，可以使用:
gmt grdproject merc_grav.nc -R-40/40/-70.0004681551/-29.9945810754 -Jm1i -I -D2m -Ggrav.nc
在某些情况下，上述操作中的两个坐标的范围不能很好的对齐。这种情况下，可以使用下面的代码:
gmt grd2xyz merc_grav.nc | \
gmt mapproject -R-40/40/-70.0004681551/-29.994581075 -Jm1i -I | \
gmt surface -R-40/40/-70/70 -I2m -Ggrav.nc
绘图示例
如果要获取上述区域的墨卡托地图，假定用户的 gmt.conf 中的 PROJ_LENGTH_UNIT 是 inch 。因为上面的 merc_grav.nc 是使用 -Jm1i 投影得到的，1度对应地图 1 英寸，因此为 80 英寸宽。可以通过使用 -Jx0.1i 投影来制作 8 英寸宽的地图以在其他程序中应用 （例如：grdcontour，grdimage，grdview），如果要叠加地理坐标图，可以使用上面的 -R 和 -Jm0.1 来使两个坐标系统匹配，避免导致移位。
除了上述的方法，还有更加方便的方法。注意到输入的 img 文件的像素宽为 2 分（可以通过 grdinfo mercgrav.nc 查看）,
merc_grav.nc 有 2400 x 2039 个像素。如果图的宽度为 8 英寸，每英寸中包含 300 个像素。但实际绘图中不需要这么高的分辨率，每英寸 100 个像素即可满足，因此可以把原始数据平均为 3 x 3 的像素。（如果想绘制等值线图的话，可能需要更多的平均保证等值线光滑，例如：6 x 6）由于 2039 不能被 3 整除，因此实际的 -R 会被调整:
gmt img2grd -M world_grav.img.7.2 -Gmerc_grav_2.nc -R-40/40/-70/-30 -N3 -V
调整以后的区域为 -R-40/40/-70.023256525/-29.9368261101 ，输出结果为 800 x 601 个像素。下面可以生成一个可以使用 grdgradient 生成一个光照文件:
gmt grdgradient merc_grav_2.nc -Gillum.nc -A0/270 -Ne0.6
如果已经有一个 “grav.cpt” 文件，可以使用下面的命令生成阴影图:
gmt begin
gmt grdimage merc_grav_2.nc -Iillum.nc -Cgrav.cpt -Jx0.1i
gmt basemap -R-40/40/-70.023256525/-29.9368261101 -Jm0.1i -Ba10
gmt end show
如果要从 img 文件中只获取地理坐标的约束数据值，则需要同时使用 -T2 选项，使用 grd2xyz 输出数据值，通过管道和``grep -v NaN`` 可以消除 NaN 值。