6.4 网格文件

GMT可以绘制2D网格数据。通常,2D网格文件的X方向代表经度、Y方向代表纬度,Z值可以表示高程、重力值、温度、速度等。也可以将 XYZ 格式的 1D 表数据网格化得到 2D 网格数据。

GMT可以使用的2D网格数据有三种格式:

  1. netCDF格式
  2. Sun光栅文件
  3. 自定义的二进制数据格式

最常见也最推荐的网格数据格式是 netCDF 格式。

6.4.1 网格文件格式

GMT支持多种格式的网格文件,包括netCDF格式、二进制网格文件以及8位Sun光栅文件。其中最常用的是netCDF格式。除此之外,GMT还可以读取用户自定义的网格文件,只要用户写好自定义网格文件的读写子程序,并将其与GMT函数库链接起来即可。详情参考源码中的 gmt_customio.c 文件。

NetCDF网格文件

GMT默认将2D网格保存为兼容COARDS的netCDF文件,一般以 .nc.grd 作为文件后缀。

COARDS是许多机构在分发网格文件时遵循的标准格式。GMT兼容该格式,因而GMT可以直接读取其他机构或程序提供的网格文件,GMT生成的网格文件也可以被其他程序读取。

netCDF格式目前有两个版本,netCDF3和netCDF4。GMT目前默认使用netCDF4版本的文件格式。

关于netCDF网格文件的详细信息,见本章的其他节。

Native二进制网格文件

旧版本的GMT不支持netCDF格式的文件,使用的是GMT自定义的二进制网格格式。新版本的 GMT依然支持这种二进制网格格式,但不建议使用。该文件格式由两部分组成:892个字节的头段区和长度不定的数据区。详细信息见 Native二进制网格格式

Sun光栅文件

Sun光栅文件格式包含了一个头段区以及一系列无符号一字节整型以表示bit-pattern。详细信息见 Sun光栅文件

6.4.2 写网格文件

上面已经说过,GMT可以读写三类格式的网格数据:netCDF格式、二进制格式和Sun光栅文件。进一步细分,GMT可以读取几十种不同格式的网格数据。

GMT所支持的网格文件格式在下表列出,每种网格文件格式都有一个对应的两字符ID。

表 6.1 网格文件格式ID
ID 说明
  GMT 4 netCDF standard formats
nb GMT netCDF format (8-bit integer, COARDS, CF-1.5)
ns GMT netCDF format (16-bit integer, COARDS, CF-1.5)
ni GMT netCDF format (32-bit integer, COARDS, CF-1.5)
nf GMT netCDF format (32-bit float, COARDS, CF-1.5)
nd GMT netCDF format (64-bit float, COARDS, CF-1.5)
  GMT 3 netCDF legacy formats
cb GMT netCDF format (8-bit integer, depreciated)
cs GMT netCDF format (16-bit integer, depreciated)
ci GMT netCDF format (32-bit integer, depreciated)
cf GMT netCDF format (32-bit float, depreciated)
cd GMT netCDF format (64-bit float, depreciated)
  GMT native binary formats
bm GMT native, C-binary format (bit-mask)
bb GMT native, C-binary format (8-bit integer)
bs GMT native, C-binary format (16-bit integer)
bi GMT native, C-binary format (32-bit integer)
bf GMT native, C-binary format (32-bit float)
bd GMT native, C-binary format (64-bit float)
  Miscellaneous grid formats
rb SUN raster file format (8-bit standard)
rf GEODAS grid format GRD98 (NGDC)
sf Golden Software Surfer format 6 (32-bit float)
sd Golden Software Surfer format 7 (64-bit float)
af Atlantic Geoscience Center AGC (32-bit float)
ei ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
ef ESRI Arc/Info ASCII Grid Interchange format (ASCII float)
gd Import/export via GDAL

除了上面列出的网格文件格式之外,有C编程经验的高级用户还可以自己自定义网格文件格式,并将读写该格式的子程序链接到GMT函数库中,使得GMT可以支持自定义网格文件格式的读取。详情见源码中的 gmt_customio.c

GMT在读网格数据时,可以自动检测网格文件的格式,所以不用担心。而在写网格数据时,GMT会默认使用 nf 格式(可以修改 IO_GRIDFILE_FORMAT 以设置其他网格文件格式为默认值)。如果想要以其他格式写网格数据,则可以在网格文件名后加上 =<ID> 参数以指定要使用的网格文件格式。

无论是读还是写网格文件,都可以按照如下格式指定网格文件的文件名:

<name>[=<ID>][+s<scale>][+o<offset>][+n<nan>]

其中

  • <name> 是网格文件名
  • <ID> 是写网格文件时要使用的网格文件格式
  • <scale> 将所有数据乘以比例因子 <scale> ,默认值为1
  • <offset> 将所有数据加上一个常数 <offset> ,默认值为0
  • <nan> 表明将文件中值为 <nan> 认为是NaN

<scale><offset> 都可以取为 a ,表明由程序自动决定比例因子和偏移量的值。在读网格文件时,会先乘以比例因子再加上偏移量;在写网格文件时,会先加上偏移量,再乘以比例因子。

举几个例子:

  1. 读入Golden软件公司的surfer软件生成的网格文件,GMT可以自动识别,故而直接用 file.grd
  2. 读一个二进制短整型网格文件,先将所有值为32767的值设置为NaN,再将数据乘以10 并加上32000,可以用 myfile.i2=bs+s10+o32000+n32767
  3. 将一个二进制短整数网格文件减去32000再除以10,然后写到标准输出,可以用 =bs+s0.1+o-3200
  4. 读一个8字节标准Sun光栅文件(其原始范围为0到255),并将其归一化到正负1范围内,可以用 rasterfile+s7.84313725e-3+o-1 ,即先乘以因子使得数据范围从0到255 变成0到2,再减去1,则数据范围变成-1到1
  5. 写一个8字节整型netCDF网格文件,偏移距由GMT自动决定,可以用 =nb+oa

GMT还支持通过网格文件后缀自动识别网格文件格式,详情见 网格文件后缀 一节。

6.4.3 读netCDF文件

netCDF格式的设计相当灵活,可以包含多个多维变量。而GMT中与网格相关的模块,只能直接处理包含一个二维变量的netCDF文件。因而,GMT在读取包含了多个多维变量的 netCDF文件时,可以做一些特殊的处理。

多个二维变量的处理

当netCDF网格文件中包含多个二维变量时,GMT默认会读取第一个二维变量作为Z值,并忽略其余的二维变量。如果用户想要自己指定读取某个特定的二维变量,可以在网格文件名后加上后缀 ?<varname> 来实现,其中 <varname> 是netCDF文件中包含的变量名。

比如想要从文件中获取名为 slp 的二维变量的信息,可以用:

gmt grdinfo "file.nc?slp"

两点说明:

  1. netCDF中包含的变量名 <varname> 可以用 ncdump -c file.nc 得到
  2. Linux下问号会被解析为通配符,因而在命令行或Bash中使用时需要将问号转义,或者将整个文件名放在单引号或双引号内

三维变量的处理

偶尔会遇到三维网格文件,比如地球参考模型,三个维度分别是经度、纬度和深度,模型中的速度和密度等则是一个三维变量。

在遇到多维变量时,GMT默认会读取第一层(即深度值最小的那一层)数据。可以通过如下两种方法来读取特定层的数据。

  1. 文件名后加上 [<index>]

    <index> 是第三维度变量(比如深度)的索引值,第一层的索引值为0

  2. 文件名后加上 (<level>)

    <level> 是要获取数据的那一层的深度值。若 <level> 指定的深度与网格不重合,则GMT会找到离其最近的深度,而不会去做插值

假设有一个地球模型文件, ncdump -c file.nc 的结果为(只列出与深度有关的部分):

dimensions:
    depth = 32 ;
variables:
    float depth(depth) ;
    depth:long_name = "depth below earth surface" ;
    depth:units = "km" ;
    depth:positive = "down" ;
data:
    depth = 50, 100, 200, 300, 400, 400, 500, 600, 600, 700, 800, 900, 1000,
        1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200,
        2300, 2400, 2500, 2600, 2700, 2800, 2850 ;

从中可以看到,该模型在深度方向上有32层,分别对应50千米、100千米,一直到2850千米。 file.nc?vp[1] 会读取第二层(即深度100 km处)的P波速度;而 file.nc?vp(200) 会读取深度200千米处的P波速度。

说明:

  1. ncdump -c file.nc 命令可以查看netCDF网格文件中的信息
  2. Linux下中括号和小括号有特殊含义,因而在命令行或Bash中使用时需要进行转义,或者将整个文件名放在单引号或双引号内

四维变量的处理

对于四维变量,方法类似。假设有一个四维网格文件,四个维度分别是纬度、经度、深度、时间,变量为压强。利用 ncdump 可以查看四个纬度的取值范围:

lat(lat): 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
lon(lon): 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
depth(depth): 0, 10, 20, 30, 40, 50, 60, 70, 80, 90
time(time): 0, 12, 24, 36, 48
pressure(time,depth,lat,lon): 共10x10x10x5=5000个值

为了得到depth=10,time=24处的变量信息,可以用:

gmt grdinfo "file.nc?pressure[2,1]"

或者:

gmt grdinfo "file.nc?pressure(24,10)"

需要注意,时间在前,深度在后。

一维变量的处理

包含一维变量的netCDF文件,也就是前面所说的netCDF表。可以通过在文件名后加上变量名来使用一个一维变量,比如:

gmt psxy "file.nc?lon/lat" ...
gmt convert "file.nc?time/lat/lon"

If one or more of the selected variables are two-dimensional, and have the same leading dimension as the other selected variables they will be plotted in their entirety. For example, if a netCDF files contains 6 time steps recording temperature at 4 points, and the variable temp is a 6 by 4 array, then the command gmt convert "file.nc?time/temp" can result in:

2012-06-25T00:00:00 20.1 20.2 20.1 20.3 2012-06-25T12:00:00 24.2 23.2 24.5 23.5 2012-06-26T00:00:00 16.1 16.2 16.1 16.3 2012-06-26T12:00:00 22.1 23.0 23.9 23.5 2012-06-27T00:00:00 17.5 16.9 17.2 16.8

If, for example, only the second temperature column is needed, use gmt convert "file.nc?time/temp[1]" (indices start counting at 0).

修改坐标单位

某些GMT工具要求网格中的两个维度的单位必须是米,若输入数据中的维度的单位不是米,则需要对网格坐标做一些变换。

  1. 如果使用的是地理网格数据(即两个维度是经度和纬度),可以加上 -fg 选项,则网格坐标会根据Flat Earth近似,自动转换成以米为单位。
  2. 若使用的是笛卡尔坐标下的网格,但维度的单位不是米(比如是千米),则可以在网格文件名后加上 +u<unit> 选项来指定当前网格的维度单位,程序会在内部自动转换成以米为单位。比如,要读入一个维度单位为千米的网格文件,可以通过 filename+uk 将其转换成以米为单位。在输出网格时,会自动使用输入数据的原始单位,除非输出网格文件名中有额外的 +u 选项。也可以使用 +U<unit> 实现逆变换,将以米为单位的网格坐标变成以 <unit> 为单位。

6.4.4 网格配准

GMT中的2D网格文件,在确定了网格范围和网格间隔后,网格线会出现在 \(x = x_{min}, x_{min} + x_{inc}, x_{min} + 2 \cdot x_{inc}, \ldots, x_{max}\)\(y = y_{min}, y_{min} + y_{inc}, y_{min} + 2 \cdot y_{inc}, \ldots, y_{max}\) 处。而节点的位置有两种选择,即网格线配准(gridline registration)和像素配准(pixel registration)。GMT默认使用的是网格线配准方式。

../../_images/GMT_grid_registration.png

图 6.1 GMT网格配准方式

(左)网格线配准;(右)像素配准。

注解

大多数原始观测数据都采样网格线配准方式,而有时经过处理的数据会以像素配准方式发布。尽管两种配准方式可以互相转换,但转换过程中会降低Nyquist采样率,阻尼一些高频信息。因而如果你可以控制,应尽量避免配准转换。

网格线配准

在网格线配准方式下,节点(图中黑色圆圈)中心位于网格线的交叉点处,节点的值代表了长宽为 \(x_{inc} \cdot y_{inc}\) 的单元(图中红色区域)内的平均值。这种情况下,节点数目与网格范围和间隔的关系为:

\[\begin{split}\begin{array}{ccl} nx & = & (x_{max} - x_{min}) / x_{inc} + 1 \\ ny & = & (y_{max} - y_{min}) / y_{inc} + 1 \end{array}\end{split}\]

左图中nx=ny=4。

像素配准

在像素配准方式下,节点(图中黑色圆圈)位于网格单元的中心,即网格点之间的区域,节点的值代表了每个单元(图中红色区域)内的平均值。在这种情况下,节点数目与网格范围和间隔的关系为:

\[\begin{split}\begin{array}{ccl} nx & = & (x_{max} - x_{min}) / x_{inc} \\ ny & = & (y_{max} - y_{min}) / y_{inc} \end{array}\end{split}\]

因而,对于相同的网格区域和网格间隔而言,像素配准比网格线配准要少一列和一行数据。右图中nx=ny=3。

6.4.5 边界条件

GMT中的某些模块在对网格文件做某些操作(比如插值或计算偏导)时,在网格边界处需要指定网格的边界条件。边界条件的选取会影响到区域边界处的计算结果。GMT中可以通过 -n 选项指定网格的边界条件。

GMT中网格文件的边界条件有三类:

默认边界条件

默认的边界条件是:

\[\nabla^2 f = \frac{\partial}{\partial n} \nabla^2 f = 0\]

其中 \(f(x, y)\) 是网格文件内的值, \(\partial/\partial n\) 是垂直于这个方向的偏导。

\[\nabla^2 = \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)\]

是二维Laplace操作符。

周期边界条件

X方向的周期边界条件表明数据是以周期 \(x_{max} - x_{min}\) 重复的,数据每 \(N = (x_{max} - x_{min})/x_{inc}\) 个点重复一次。Y方向同理。

  • 对于网格线配准的网格文件,共N+1列数据。第一列数据位于 \(x = x_{min}\) 处,最后一列(N+1列)数据位于 \(x = x_{max}\) 处,周期边界条件意味着数据的第一列和最后一列是完全相同的
  • 对于像素配准的网格文件,有N列数据,第一列位于 \(x_{min} + x_{inc}/2\) ,最后一列(第N列)位于 \(x_{max} - x_{inc}/2\) ,第一列和最后一列的数据是不同的。

地理边界条件

地理边界条件表明:

  1. \((x_{max} - x_{min}) \geq 360\) 且180是 \(x_{inc}\) 的整数倍,则在X方向使用周期为360的周期边界条件,否则使用默认边界条件
  2. 若条件1为真,且 \(y_{max} = 90\) 则Y方向上使用“北极边界条件”,否则使用默认边界条件
  3. 若条件1为真,且 \(y_{min} = -90\) 则Y方向上使用“南极边界条件”,否则使用默认边界条件

6.4.6 查看netCDF文件

某些软件可以直接用于查看netCDF文件的内容:

更多相关工具,见 netCDF网站上的列表

注解

尽管大多数程序都可以读取 netCDF 文件,但某些不支持 netCDF4 格式。