輸出 binary 檔案格式：跨平台應用

在 輸出 .nc 檔案格式 一文中介紹了如何輸出 netCDF 格式，筆者認為是NCL作業平台上檔案輸出格式的最佳選擇，但若輸出檔需要在其他平台上應用 (e.g., Matlab, Fortran, GrADS, GMT, Python, ...)，熟悉通用的 binary 格式是必要的。

本篇會示範 NCL - Fortran - Matlab 三者間如何互相讀寫檔案(傳送資料)。

窺探 binary 檔案格式*的奧秘應該要了解：
<1> **binary 格式和 ascii 格式有何不同？
<2> **計算機(電腦)如何儲存資料？
    -- precision: single v.s. double
    -- access: sequential v.s. direct
    -- endianness: big v.s. little (v.s. middle)
*註1：這邊強調的是浮點數(float number)的儲存。浮點數，又稱為科學記號，是一般研究資料最常見的形式(相較於整數(integer)和文字(string))。
**註2：之後再分享關於這兩個問題的心得 (另闢文章)。

如何讀寫 binary 檔案，筆者覺得可分為兩個層次來學習。

第一、追根究柢：理解了上述兩個問題的答案後，知道計算機如何儲存資料，然後學習每種型式的資料如何用不同的程式語言(e.g., NCL, Fortran, Matlab, ...)來表達。

第二、不求甚解：若一下子跳進「計算機概論」式的內容怕模糊了焦點：不是要示範讀寫檔案嗎？怎麼會在介紹 0110 1010 0101？因此秉持先求有再求精的精神，這邊將要以你應該要最熟悉的 Fortran 為主幹，示範NCL和Matlab要如何與之配合。

等習得「不求甚解式」後，再回頭看「追根究柢式」的祕笈，或許會突然打通任督二脈一切豁然開朗。

從 Fortran 開始

假設 var (nLon, nLat, nLevel, nTime) 是一個四維單精準浮點數的矩陣。
這樣宣告
REAL(KIND=4) var (nLon, nLat, nLevel, nTime)

這裡沒有要做 Fortran 的 step-by-step 教學，所以只簡單的複習。在 OPEN 檔案時最好將 ACCESS 和 ENDIANNES 都手動指定，確保知道自己在做什麼 (不手動指定就會變成預設值)。

將 var 存成一個名為 bFile_seq.dat 的檔案和一個 bFile_dir.dat 的檔案。
OPEN ( 51 , FILE = 'bFile_seq.dat' , FORM = 'UNFORMATTED' , CONVERT = 'LITTLE_ENDIAN' , &
        ACCESS = 'SEQUENTIAL' , STATUS = 'UNKNOWN' )
OPEN ( 61 , FILE = 'bFile_dir.dat' , FORM = 'UNFORMATTED' , CONVERT = 'LITTLE_ENDIAN' , &
        ACCESS = 'DIRECT' , RECL = 4*nTime*nLevel*nLat*nLon , STATUS = 'UNKNOWN' )

若要使用 ascii 檔案就改成 FORM = 'UNFORMATTED' 並刪除 ACCESS 和 CONVERT 。
若要使用 big endian 就改成 CONVERT = 'BIG_ENDIAN' 。

Sequential access 檔案寫入/讀出的方法如下
WRITE (51) var
READ (51) var

Direct access 檔案寫入/讀出的方法如下
WRITE ( 61 , rec = 1 ) var
READ ( 61 , rec = 1 ) var

回到 NCL

在 NCL 中，讀/寫檔案是兩個不同的函數，而兩種 access 也分別由兩種函數負責，共四種組合。指定 endianness 則是由 setfileoption 來處理。

Sequential access: fbinrecwrite / fbinrecread
Direct access: fbindirwrite / fbindirread
Endianness:
setfileoption ( "bin" , "WriteByteOrder" , "LittleEndian" )
setfileoption ( "bin" , "ReadByteOrder" , "LittleEndian" )

同樣拿 var (nTime, nLevel, nLat, nLon) 來舉例

寫入/讀出：
<<Sequential>>
fbinrecwrite ( "bFile_seq.dat" , -1 , var )
var = fbinrecread ( "bFile_seq.dat" , 0 , (/ nTime , nLevel , nLat , nLon /) , "float" )
<<Direct>>
fbindirwrite ( "bFile_dir.dat" , var )
var = fbindirread ( "bFile_dir.dat" , 0 , (/ nTime , nLevel , nLat , nLon /) , "float" )

在 Fortran 中只有讀寫 direct access 時需要指定位置編號 (i.e. rec)，sequential access 的時候都是用排隊的方式，第一個讀(寫)完讀(寫)第二個，一個接一個，因此不用指定位置號碼。但是在 NCL 中，範例中即顯示無論是 sequential 或 direct 的存取都需要指定位置 (範例中的那個 0 表示第一筆)。這其實是優點不是缺點...表示在NCL中即使讀取 sequential access 的檔案也不用每次都從頭讀到尾，可以只跳著讀取需要用到的部分(就像Fortran處理direct access一樣的概念)，當資料檔案很大，讀取很耗時間的時候，會非常實用。

Sequential / Direct Access 有何不同？

關於兩者的差別，應該連同著前面的兩個提問一起解釋會有更清楚的全盤了解。這裏先就概念式的描述來說明。

在binary格式中，單精準度的數字(無論是浮點數或整數)，每一個數字電腦使用了 4 byte 來儲存，2^10 (1024) 個 byte 就是「1 kb」，而 2^20 (1024^2) 個 byte 就是大家熟悉的「1 mb」。

假設現在有一筆2天的資料，每天有5個數字要記錄，所以總共是10個數字 (e.g., var(5,2))。

在 Fortran 中，無論是sequential或direct都可以一次寫入，亦即
WRITE(51) var
WRITE(61 , rec=1) var
也可以分批寫入 (例如每一筆時間寫入一次, 常用!)
WRITE(52) var(:,1)
WRITE(52) var(:,2)
WRITE(62 , rec=1) var(:,1)
WRITE(62 , rec=2) var(:,2)

對於direct access檔案來說，分次寫入是沒有差別的，也就是說檔案61號和檔案62號的大小都會是 10 * 4byte = 40 byte。如果寫入順序對的話(像這邊是正確的)，檔案61號和檔案62號兩個檔案會是identical的。

所謂的寫入順序正確，對於Fortran來說，維度的處理順序是由前而後，亦即檔案61號雖然只有一行指令寫入，但在寫入的時候這10個數字的順序是(1,1),(2,1),(3,1),...(5,1)，寫完第一個維度的5個數字後，再寫第二個維度 (NCL剛好相反!!!)。

而sequential access的方式，會在寫入的「資料塊」前後各加一個整數(佔用 4 byte 的硬碟空間)，記錄被夾起來的這一塊區塊總共佔用了幾個byte。舉例來說，檔案51號一次寫入了10個數字，所以檔案大小為 4 + (4*10) + 4 = 48 byte。前後兩個 4 byte 都存了「40」這個整數，一個單精準度的浮點數值需要 4 byte 來儲存，所以 40 代表這一筆資料總共有 10 個值。

這也是為什麼 Fortran 在讀取sequential資料的時候，不用事先指定變數的維度大小 (direct會在OPEN檔案的時候用 RECL=4*維度 來指定資料區塊大小，但sequential不用)，因為當讀取sequential檔案的時候，Fortran會先讀取第一個 4 byte，這樣就知道需要往下繼續讀多少個 byte，然後尾端的 4 byte 可以用來檢查資料長度是否正確 (因為頭的 4 byte 和尾的 4 byte 理當紀錄著同一個整數)。

舉一反三，檔案52號因為寫了兩次，每次5個值，所以檔案大小應為 [ 4 + (4*5) + 4 ]*2 = 56 byte。而前後的 4 byte 儲存的是「整數值 20」。

衍伸到 Matlab

Matlab 當然也可以讀寫 binary 檔案！方式跟 Fortran 中的 Direct Access 檔案非常類似。

寫入
fid1 = fopen ( 'for_write.bin' , 'w' )
fwrite ( fid1 , var , 'float' )
fclose ( fid1 )

讀出
fid2 = fopen ( 'for_read.bin' , 'r' )
var = fread ( fid2 , nX*nY*nZ*nT , 'float' )
fclose ( fid2 )

fid 就像 Fortran OPEN 檔案時指定的檔案編號。寫入資料的時候，不用給維度資訊，反正 var 有幾個值就寫入多少個值。讀出的時候，因為是 direct access 的格式(和概念)，頭尾並沒有 4 byte 的資訊說明後面有幾個值要讀，所以必須給定維度訊息。

那 Matlab 能夠存/取 sequential access 的檔案嗎？
答案是可以的，但是非常的自找麻煩！

由於 Matlab 的 binary 檔案讀寫屬於 direct access，不像 Fortran 或著 NCL 本來就有現成的「sequential 指令」可以用，所以如果要存取 sequential 檔案的話，使用者必須手動處理每個資料塊前後的那個 4 byte 整數值。單精準度 (用 4 byte 來儲存) 的整數在 Matlab 中以 'int32' 表示，其中的 32 代表這是 32 位元的系統所使用的整數儲存方式。

同樣的以 var(5,2) 為例，如果要用 Matlab 寫成一個 sequential 的 binary 檔案，看起來會像這樣
fid3 = fopen ( 'for_write_sequential.bin' , 'w' )
fwrite ( fid3 , 40 , 'int32' )
fwrite ( fid3 , var , 'float' )
fwrite ( fid3 , 40 , 'int32' )
fclose ( fid3 )

必須手動在資料前後加上 40 這個數字，以整數 (int32) 方式儲存。40 代表 40 byte，代表總共有 10 個單精準度 (4 byte) 的值。

照過來～範例在這裡！

如果你一行一行的仔細瀏覽到這裡還沒被搞糊塗，那就算輸給你了...因為連筆者都快被自己弄糊塗了。接下來為了更清楚的表達 binary 的讀寫方法，下面提供 Fortran, NCL, Matlab 三種程式各兩支範例，一支示範 write binary data，另一支示範 read binary data，總共有 6 支範例。

每支 write binary 的範例都會建立 4 個 binary 檔案，分別是單次和多次寫入的 sequential 和 direct 檔案。
Fortran_write_example.f90 執行完會產生
    fortran_seq_single.dat
    fortran_dir_single.dat
    fortran_seq_multiple.dat
    fortran_dir_multiple.dat
NCL_write_example.ncl 執行完會產生
    ncl_seq_single.dat
    ncl_dir_single.dat
    ncl_seq_multiple.dat
    ncl_dir_multiple.dat
Matlab_write_example.m 執行完會產生
    matlab_seq_single.dat
    matlab_dir_single.dat
    matlab_seq_multiple.dat
    matlab_dir_multiple.dat

這三組資料是 identical 的 (e.g., fortran_seq_single.dat，ncl_seq_single.dat，matlab_seq_single.dat 三個檔案一模一樣)。在 Linux 底下可以用 diff 指令確認兩個檔案是否完全相等：

>>$ diff fortran_seq_single.dat ncl_seq_single.dat

如果兩個檔案不是 identical 的內容，電腦會告訴你這兩個檔案不相同。

另外三支範例檔案，都是用來讀取這一組 (四個) binary 檔案用的。
Fortran_read_example.f90
NCL_read_example.ncl
Matlab_read_example.m

因為這三組資料一模一樣，所以每一支讀檔範例都可以自己改檔名套用到任何一組資料上面。

這六支範例檔案一次打包下載。

在範例中，使用的變數 (var) 有四個維度，假設為常用的 X (lon) , Y (lat) , Z (level) , T (time)，每個維度的大小分別為 5, 4, 3, 2，所以 var 總共有 5*4*3*2 = 120 個值。數值設計了一些玄機，X 的大小表示個位數字，Y 的大小表示十位數字，Z 的大小表示百位數字，T 的大小表示千位數字。在這樣的設計下 var 的 120 個數值都不會重複，舉例來說 var(2,1,3,1) 的值是 1312。這樣設計的好處是方便除錯，舉例來說當讀到 2233 這個值的時候，位置應該在 var(X=3,Y=3,Z=2,T=2)，如果不是，那就是某個地方錯了！

這裡強調如何給定 var 的值，目的在於凸顯 Fortran, NCL, Matlab 三者對於矩陣的排列方式的差異。簡單的說，Fortran 和 Matlab 類似，座標排列方式為 (X, Y, Z, T)；NCL 和前面兩者相反，座標排列方式為 (T, Z, Y, X)。這點，在跨平台轉換資料的時候尤其重要，不可忘記。至於處理上的細節，直接參考範例中的操作方式應該更為容易理解。

最後提醒一點，在 Matlab 中讀取 binary 檔案中的變數後，都是一維的向量 (長度就會是 nX*nY*nZ*nT)，必須再透過 reshape 去把這個一維向量轉換成四維矩陣。

var_after = reshape ( var_before , [nX, nY, nZ, nT] )

在看完六支範例之後，如果能再重新從頭到尾瀏覽一次這篇文章，應該會更理解筆者想要表達什麼... (Sorry, 看起來很亂 但筆者已經盡力了...)

Take Home Message

若有跨平台需求，不要鐵齒，請選用 Direct Access 就好。

Fortran
    OPEN
    WRITE
    READ
NCL
    setfileoption
    fbindirwrite
    fbindirread
Matlab
    fopen
    fwrite
    fread
    fclose