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 值。