讀netCDF文件

netCDF格式的設計相當靈活,單個netCDF文件中可以包含多個多維變量。而GMT只能直接處理包含單個二維變量的netCDF文件。因而,對於單變量二維netCDF文件,直接給出文件名即可;而對於複雜的多變量多維netCDF文件,則需要用戶在指定文件名時給出額外的信息。

讀二維單變量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

讀取二維多變量netCDF文件

對於包含多個二維變量的netCDF網格檔,GMT默認會讀取第一個二維變量作爲Z值,並忽略其餘的二維變量。用戶可以通過在網格檔名後加上後綴 ?<varname> 的方式指定要讀取某個特定的二維變量,即:

<filename>?<varname>

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

比如想要從文件中獲取名爲 slp 的二維變量的信息,可以用:

gmt grdinfo "file.nc?slp"

注解

Linux下問號會被解析爲通配符,因而在命令行或Bash中使用時需要將問號轉義,或者將整個文件名放在單引號或雙引號內。

讀取三維單/多變量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中使用時需要進行轉義,或者將整個文件名放在單引號或雙引號內

讀取四維單/多變量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)"

在本例中,時間維度在前,深度維度在後。

讀取一維單/多變量netCDF文件

一維單/多變量netCDF文件,即前面所說的以netCDF格式保存的表數據。即表數據中的每一列分別保存爲netCDF文件中的一個變量。GMT自帶的GSHHG數據和DCW數據就是一維多變量netCDF文件。

同樣的,可以使用 ncdump -c file.nc 來查看netCDF文件所包含的變量名。然後即可通過在文件名後加上一系列用斜槓分隔的變量名來使用這些一維變量。例如:

# 將文件中的lon變量和lat變量作爲輸入數據的第1和2列
gmt plot "file.nc?lon/lat" ...
gmt plot "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]

修改座標單位

某些GMT模塊要求網格中的兩個維度的單位必須是米,若輸入數據中的維度的單位不是米,則需要對網格座標做一些變換。

  1. 如果使用的是地理網格數據(即兩個維度是經度和緯度),可以加上 -fg 選項,則網格座標會根據Flat Earth近似,自動轉換成以米爲單位。
  2. 若使用的是笛卡爾座標下的網格,但維度的單位不是米(比如是千米),則可以在網格檔名後加上 +u<unit> 選項來指定當前網格的維度單位,程序會在內部自動轉換成以米爲單位。比如,要讀入一個維度單位爲千米的網格檔,可以通過 filename+uk 將其轉換成以米爲單位。在輸出網格時,會自動使用輸入數據的原始單位,除非輸出網格檔名中有額外的 +u 選項。也可以使用 +U<unit> 實現逆變換,將以米爲單位的網格座標變成以 <unit> 爲單位。