img2grd

贡献者:

周茂


官方文档:

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

相关模块

gmt, img2google