6.2 读netCDF文件

netCDF格式的设计相当灵活，单个netCDF文件中可以包含多个多维变量。而GMT只能直接处理包含单个二维变量的netCDF文件。因而，对于单变量二维netCDF文件，直接给出文件名即可；而对于复杂的多变量多维netCDF文件，则需要用户在指定文件名时给出额外的信息。

6.2.1 读二维单变量netCDF文件

对于二维单变量netCDF文件，用户只需要直接给出网格文件的文件名，GMT可以自动检测 netCDF文件的格式并读入。具体来说，可以按照如下格式指定网格文件文件名:

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

其中

  • <name> 是网格文件名，必须指定，其它均是可选项
  • <ID> 显式告诉GMT当前文件的格式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

6.2.2 读取二维多变量netCDF文件

对于包含多个二维变量的netCDF网格文件，GMT默认会读取第一个二维变量作为Z值，并忽略其余的二维变量。用户可以通过在网格文件名后加上后缀 ?<varname> 的方式指定要读取某个特定的二维变量，即:

<filename>?<varname>

其中 <varname> 是netCDF文件中包含的变量名，其可以通过netCDF提供的命令 ncdump -c file.nc 得到。

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

gmt grdinfo "file.nc?slp"

注解

Linux下问号会被解析为通配符，因而在命令行或Bash中使用时需要将问号转义，或者将整个文件名放在单引号或双引号内。

6.2.3 读取三维单/多变量netCDF文件

最常见的三维单/多变量netCDF文件是地震成像得到的三维地球参考模型。三个维度分别是经度、纬度和深度，变量通常是P波速度、S波速度等。

在遇到三维单/多变量netCDF文件时，GMT默认只读取第一个变量的第一层数据（通常是深度值最小的那一层）。此时，可以将其当做一个变量数组，并通过如下两种方式指定读取特定层的数据。

  1. 变量名后加上 [<index>] 以指定某一层的索引值。第一层的索引值为0，第二层的索引值为1，依次类推。
  2. 变量名后加上 (<level>) 以指定获取深度为 <level> 处值。若网格文件中在 <level> 指定的深度处并不存在数据，则GMT会找到离 <level> 最近的有数据的那一深度的值，而不会去做插值。

假设有一个地球模型文件，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千米。则可以使用如下命令:

# 读取第二层（即深度100km）处的P波速度
gmt grdinfo "file.nc?vp[1]"

# 读取深度200千米处的P波速度
gmt grdinfo "file.nc?vp(200)"

注解

Linux下问号、中括号和小括号有特殊含义，因而在命令行或Bash中使用时需要进行转义，或者将整个文件名放在单引号或双引号内

6.2.4 读取四维单/多变量netCDF文件

对于四维单/多变量netCDF文件，处理方法类似。假设有一个四维单变量netCDF文件，四个维度分别是纬度、经度、深度、时间，变量为压强。利用 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个值

此时可以将变量 <pressure> 当成一个二维数组。

为了得到depth=10，time=24处的变量值，可以用:

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

或者:

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

在本例中，时间维度在前，深度维度在后。

6.2.5 读取一维单/多变量netCDF文件

一维单/多变量netCDF文件，即前面所说的以netCDF格式保存的表数据。即表数据中的每一列分别保存为netCDF文件中的一个变量。GMT自带的GSHHG数据和DCW数据就是一维多变量netCDF文件。

同样的，可以使用 ncdump -c file.nc 来查看netCDF文件所包含的变量名。然后即可通过在文件名后加上一系列用斜杠分隔的变量名来使用这些一维变量。例如:

# 将文件中的lon变量和lat变量作为输入数据的第1和2列
gmt psxy "file.nc?lon/lat" ...

# 将文件中的变量time、lat和lon分别作为输入数据的三列
gmt convert "file.nc?time/lat/lon" ...

如果要使用的变量是一个二维变量，并且其优先维度与其他被选变量相同，则该变量整体会被输出。例如，一个netCDF文件中包含6个时间步，其记录了4个点的温度。则变量 temp 是一个6x4的数组，因而使用如下命令会输出如下信息:

$ gmt convert "file.nc?time/temp
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
2012-06-27T12:00:00 27.2 27.2 27.5 27.5

如果只需要某个点的温度，例如第二列数据，则可以使用:

$ gmt convert "file.nc?time/temp[1]

6.2.6 修改坐标单位

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

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