Day 2 方便的netCDF格式和NCL的變數結構

何謂netCDF？
相較於ascii格式，一般科學資料大多以binary格式儲存，因為可以有效的縮小檔案以利於節省硬碟空間及資料傳輸時間。如果你還不曉得ascii/binary兩種格式的區別，請容我偷吃步的解釋：Windows底下的記事本、Mac底下的TextEdit、或Linux底下的VI能直接打開看得到內容的檔案都是ascii，反之亦然。舉例來說，強制用vi開啟一個binary格式的檔案，會像下面這樣出現亂碼：

由於ascii資料佔空間，科學上處理龐大的資料量時一般都採用binary檔案。存取和寫入 (I/O) binary檔案變成科學研究的一項基本門檻。有鑒於此重要性，將會有另外一篇專門的binary格式介紹，包含如何跨平台、跨軟體的I/O等等。

作為NCL101入門篇的第二章，這裏先介紹一種非常方便的binary格式：netCDF。之所以會從這開始介紹，是因為使用NCL讀取netCDF格式的步驟非常簡單，而netCDF也是一般大氣模式常見的輸出檔之一；再者，很多衛星資料採用的是一種叫做HDF (Hierarchical Data Format) 的格式儲存，而對於NCL來說讀取netCDF的方法和HDF大同小異。因此第二章最主要的目的，就是示範如何在NCL中讀取netCDF檔案中的一筆資料。

Step 1. 取得一份netCDF檔案
一般來說，netCDF檔案的副檔名為.nc。首先，到美國NOAA網站下載NCEP Reanalysis資料，下面的說明以連結中的2004年氣溫資料 (air.2014.nc) 為例。

Step 2. 查看資料內容
在安裝NCL時，有一個好用的小工具 ncl_filedump 也會一併被安裝。這個工具可以用來快速檢視一個.nc檔案的內容。

>>$ ncl_filedump air.2014.nc

現在就可以知道，剛剛下載的air.2004.nc檔案裡面共有五個變數，分別是日均溫 (air) 和紀錄座標資訊的 level, lat, lon, time。再仔細看可以發現這是一筆經度從0到357.7，緯度從-90到90，高度從1000到10的全球資料，網格共有144*73*17個點，時間上有365筆資訊。知道所需變數在.nc檔案裡面的名稱後 (e.g., air)，就可以繼續進行下一步，編寫script了。

Step 3. 在script中讀取.nc
方式如下：NCL101-2.ncl

;===== NCL101-2.ncl =====
begin
        ncFile          = addfile("air.2014.nc","r")
        Temp            = ncFile->air
        printVarSummary(Temp)
        print(Temp@units)
        print(Temp!0)
        print(Temp!1)
        print(Temp!2)
        print(Temp!3)
;       print(Temp&lon)
;       print(Temp)
print("All Done.")
end
;=====

把下載下來的 air.2014.nc 和 NCL101-2.ncl 放在同一個資料夾內，並執行

>>$ ncl NCL101-2.ncl

ncFile = addfile("air.2014.nc","r")
告訴電腦 air.2014.nc 是一個已存在的netCDF檔案，並準備讀取 ("r", read) 它。

Temp = ncFile->air
要從檔案中讀取變數，必須先知道變數的名稱 (e.g., air)，可以用Step 2的方法先找出所需的變數名稱。而在NCL中讀取變數，就是這一行，就是這麼簡單！

printVarSummary(Temp)
這是一個非常實用的函數，可以知道 Temp 這個變數的結構和屬性。Temp 總共有4個維度 (dimensions) 和15個屬性 (attributes)。

從這裡就可以知道netCDF方便的地方。一行指令讀取進來的變數已經包含了它的座標軸和無效值 (_FillValue) 了。之後在畫圖的時候，不用另外指定，NCL就知道每一筆溫度應該點在地圖上的哪個點；而且缺資料的時候，也就是當某個網個的溫度值等於無效值的時候，NCL會自動地跳過不畫，計算時 (e.g., 平均)，NCL也會自動略過無效點。

變數的結構 ～ !, &, @
<1> 座標
NCL中座標的習慣用法，跟"一般認知"的順序相反。由例子中可看出(NCEP Reanalysis是個常見的資料庫)，變數的維度是 time * level * lat * lon，而一般對於四維坐標常見的順序應該是 X*Y*Z*T。一個變數的座標資訊，分成兩個層次：

    一、維度(座標軸)的名稱
    座標軸的名稱用「!」來存取。舉例來說，時間(time)是Temp的第一個座標，因此可由 
    print(Temp!0)
    來取得。可以用同樣的方法來改變或指定一個變數的座標軸名稱。例如：
    Temp!0 = "day"
    Temp!1 = "pressure"

    二、座標軸的值
    座標軸的值用「&」來存取。舉例來說，Temp在經向上面有144個點，
    但是資料是0至360度還是-180至180度呢？每格有等間距嗎？此時可用
    print(Temp&lon)
    來把經度座標全部顯示出來。相同的，可以用此方法來給定座標的數值。例如：
    Temp&lon = fspan(0 , 360-2.5 , 144)

在讀取.nc檔案的時候，多半不用理會座標軸名稱和值，因為都已經事先設定好了。但是當在做計算有新變數產生時，或著有時候讀取一些"純"binary檔案，必須手動設定其變數的座標資訊時，「!」和「&」就是非用不可的功能了。

<2> 陣列 (array) 的index
從範例 NCL101-2.ncl 中使用的是 Temp!0 到 Temp!3 可以看出在NCL中，陣列第一項的index不是1而是0。一個含有17個 elements 的陣列 (e.g., level)，index 是0-16，而不是1-17。level是Temp的第二個維度，但是Temp!2是lat，Temp!1才是level。

這很重要！真的很重要！超級！所以講三次，要特別記住。因為當你跟Matlab等軟體交互使用的時候，有時候會錯亂 (Matlab的index從1開始，而且Matlab的座標順序是"一般"的那種：X*Y*Z*T)，等之後介紹binary格式與平台的轉換時，會有更詳細的資訊。

<3> 屬性 (attributes)
大部份的屬性都是給使用者看的，提醒使用者變數的來源、處理過程、精準度...等等。NCL只會拿少部分的特殊屬性用，像是無效值(_FillValue)不會參與計算；或是在畫圖的時候，long_name 和 units 被預設用來放在圖的左上和右上角。

就像範例中示範的，屬性可以用「@」來存取。舉例來說
print(Temp@units)
可以看到Temp這個變數的單位。也可以透過這個方式新增/修改某項屬性
Temp@long_name = "Air Temperature from NCEP"
Temp@today = "2015/03/09"
也可以將屬性刪除
delete(Temp@var_desc)





待續... Day 3