2015年3月10日

Day 3 畫圖的基本元素:以全球氣溫分佈圖為例

接續 Day 2 的介紹,這裡使用NCEP Reanalysis 2014年的全球氣溫資料 (air.2014.nc) 示範如何用NCL畫圖並介紹其scrip的基本架構。方式如下:NCL101-3.ncl

;===== NCL101-3.ncl =====
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

; Work Station
        wks                     = gsn_open_wks( "png" , "AirTempMap" )

; Variable
        ncFile                  = addfile( "air.2014.nc" , "r" )
        Temp                    = ncFile->air

; Plot control
        res                     = True

; Plot command
        plot                    = gsn_csm_contour_map_ce( wks , Temp( 0 , 0 , : , : ) , res )

print("All Done.")
end
;=====

>>$ ncl NCL101-3.ncl

在終端機中輸入上面那段指令,執行 NCL101-3.ncl 後,圖就畫出來了!


這支超簡化的script清楚的展現畫圖所需的基本架構,包含:
1. 定義 WorkStation (wks),也就是畫布。
2. 定義變數,可以從檔案讀取,計算結果,或自己創造。
3. 作圖參數,設定各種畫圖的細節。
4. 把圖「畫」到WorkStation上。

在上面的範例中扣掉comments,只要10行的程式碼,就可以完成一張最基本的氣溫分佈圖。接下來為了更清楚的分項說明基本架構中的四個部分之用途與使用方法,請參考這個進階版的範例:NCL101-4.ncl

;===== NCL101-4.ncl =====
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

; Work Station
        wks_type                = "pdf"
        ; set the aspect ratio to 16:9, fitting your screen in full
        wks_type@wkPaperWidthF  = 8.0
        wks_type@wkPaperHeightF = 4.5
        wks                     = gsn_open_wks(wks_type,"AirTempMap2")

; Variable
        ncFile                  = addfile( "air.2014.nc" , "r" )
        ; load only part of variables, increase effeciency
        Temp                    = ncFile->air( 0 , : , : , : )

; Plot resources
        res                     = True
        res@gsnMaximize         = True
        res@gsnPaperOrientation = "portrait"

        ; zoom in to Asia
        res@mpCenterLonF        = 180.
        res@mpMinLatF           = 0.
        res@mpMaxLatF           = 80.
        res@mpMinLonF           = 60.
        res@mpMaxLonF           = 270.

        ; for shading plot
        res@cnFillOn            = True
        res@cnLinesOn           = False

; Plot command
        plot                    = gsn_csm_contour_map_ce(wks,Temp( 0 , : , : ),res)
        plot                    = gsn_csm_contour_map_ce(wks,Temp( 5 , : , : ),res)
        plot                    = gsn_csm_contour_map_ce(wks,Temp(10 , : , : ),res)

print("All Done.")
end
;=====

>>$ ncl NCL101-4.ncl

執行後會得到一份含有3頁圖片的PDF檔案 (AirTempMap2.pdf),每頁看起來像這樣:


Part 1 WorkStation
畫圖前必須先定義(創造)一個畫布,告訴NCL這個畫布的格式是PNG, JPG, PDF, EPS, ..., 還是直接透過X11視窗顯示。也可以在這裡指定頁面的大小、解析度、長寬比、橫/直幅、等等。個人覺得最方便於研究目的的格式是PDF檔案,因為 <1> 在做研究的過程中,常常需要一次畫大量的圖,不管是時間序列,或是多張垂直分層、一排 Cross Section 等等,為的是用平面(2D)圖了解多維結構。使用PDF的好處就是這一堆圖會被集合在單一的檔案中,在個人電腦上可以很方便的透過滑鼠滾輪或鍵盤方向鍵切換頁面,達成類似動畫的效果。<2> PDF和EPS一樣屬於向量圖檔,可以無限 zoom in 不失真 (無鋸齒)。這對於研究階段而言非常實用,可以只畫一張全區域圖,但需要 zoom in 時又不失細節。

在範例中可以看出,定義畫布的函數是 gsn_open_wks,第一個input是格式,第二個input是檔案名稱。對於PDF格式的workstation而言,總共有「這些」參數可以調整,設定的方法就是透過「@」指定第一個 input (e.g., wks_type) 的屬性。函數 gsn_open_wks 會依照第一個input的屬性去建立畫布規格。

Part 2 Variables
從NOAA下載下來的 air.2014.nc 裡頭總共有一整年的氣溫資料,假設現在只需要畫出第一天(1/1/2014)的氣溫資料,那麼若用指令 ncFile->air 一口氣讀取了一整年的資料將會拖慢工作效率。避免資源的浪費和增加執行速度的方法就是只讀取需要的部分,透過指定維度,可以只讀取第一天的三維(Z*Y*X)資料。因為 air 儲存在 .nc 檔案中有四個維度,分別為 ( time , level , lat , lon )。只讀取第一筆時間,並包含所有 X, Y, Z 的資料的方法是
Temp = ncFile->air( 0 , : , : , : )
「0」代表第一個index,也就是第一筆時間。「:」代表全部。如果只要最下面三層的資料:
Temp = ncFile->air( 0 , 0:3 , : , : )
讀完資料後,可以用 printVarSummary(Temp) 檢查所讀出來的變數資訊。

這本章的兩個範例中,變數 Temp 都是直接從 .nc 檔案裡面讀取出來的,變數本身就已經包含了座標的資訊在裡頭,因此畫圖的時候無需額外指定X軸與Y軸,NCL便會自動的判斷並做出最佳化的設定。然而,如果要畫圖的變數是經過一連串計算的結果 (e.g., 前後兩個時間點的溫度變化),或著變數是從其他格式的binary檔案中讀取出來的,在此狀況下,變數本身不一定含有座標軸的資訊,必須透過手動的方式指定,方法在 Day 2 中提過,需利用「!」和「&」。判斷一個變數是否有合適的座標軸資訊,一樣可以用 printVarSummary 檢查。

Part 3 Plot Resources
在NCL中,所有畫圖的設定 (resources) 都是透過一個邏輯變數 (logical variable, true or false) 來控制。
res = True
res = False
可以控制這張圖要不要使用使用者客製化的設定,如果是 False 就全部都用預設 (default) 設定。而一個個的控制參數就用屬性的方式透過「@」加上去。

舉例來說,res@mpMinLatF, res@mpMaxLatF, res@mpMinLonF, res@mpMaxLonF 可以控制地圖的經緯度範圍。這裏示範了一個常見的問題,就是NCL在預設上,地圖是以0度經度當做中心,亦即經度的範圍是-180至180度。如果畫一張大西洋的圖沒有問題,可以設定mpMinLonF為-100、mpMaxLonF為0。但是若要畫一張太平洋的圖,由於畫面會橫跨換日線 (180度),因此必須告訴NCL要把經度-180度和經度180度的資料拼起來,經度的部分要看成0至360度,否則範例中的設定就會變得不合邏輯:從60度畫到-90度?mpMinLonF 比 mpMaxLonF 還要大?到底是要畫 60~270 度還是 -90~60 度?因此得將地圖中心移至180度。
res@mpCenterLonF = 180.

[屁話區] 既然要畫 60~270 度,怎麼不是把 mpCenterLonF 設定在畫面正中央的 165 度,而設 180 度?當然要設成 165 也是可以... 是可以拉... 就是懶而已。設中心為 180 度的意思是告訴NCL現在經度從0算到360,左 (右) 邊界是 0 (360) 度。而設中心為 165 度則是將左 (右) 邊界設為 -15 (345) 度。只要畫圖範圍不要「跨越」邊界就可以了。

另一個設定是將 NCL101-3.ncl 所畫的 contour 塗滿變成 NCL101-4.ncl 所畫的 shading,並將輪廓線刪除,只留下色塊的部分。
res@cnFillOn = True
res@cnLinesOn = False
關於 contour 圖和 shading 圖的作圖細節會在之後進一步討論,這章要強調的是整體大架構的概念。還記得在NCL中,大寫和小寫是被當作不同的東西在看待嗎?因此這邊的各項設定例如 mpMinLonF 不能隨意改成 MPminLONf。而其實在這些大小寫交錯出現的字串中,隱含著很多有用的訊息。

從前面介紹的部分不難發現,mp 開頭的設定都跟地圖 (map) 有關,要設定等值線 (contour) 的相關參數就會是 cn 開頭。還有兩個沒介紹到的由 gsn 開頭的設定,是跟整張圖有關的一些進階設定。所以第一個子字串 (sub-string) 是大分類,接下來的子字串,長度不依,但都是由大寫字母作為區隔,都有其代表的意思。例如
mpCenterLonF:地圖,中心,經度,數值
cnFillOn:等值線,填滿,開關
一般來說,以 F 結尾的參數,其值多為浮點數;以 On 結尾的參數,是開關,所以通常是邏輯變數 True / False。NCL以這樣的方式設計並分類其數以千計的參數,不但利於讀者看懂程式碼,也便使用者搜尋所需的參數。

Part 4 Plot Command
當設定好了 WorkStation 和 Resources 後,最後真正要「畫」圖的動作其實非常簡單,通常只需要一行指令就可以解決了。在NCL中把各種種類的圖,都各自獨立出一個專屬的函數,好處就是執行「畫」這個動作的時候,作圖函數的 input 通常只需要三個:Where (WorkStation),What (Variables),How (Resources)。像此範例中的全球氣溫分佈:
gsn_csm_contour_map_ce( wks , Temp( 0,:,:) , res )

範例中的 gsn_csm_contour_map_ce 用來處理平面地圖,ce 代表 "Cylindrical Equidistant",所謂的等圓柱投影。

如果要畫箭頭向量 (e.g., 風速風向),可以用
gsn_csm_vector_map_ce

如果要換一種投影方式,例如極投影,可以用
gsn_csm_contour_map_polar

如果要畫垂直頗面,譬如想看 Hadley Cell
gsn_csm_pres_hgt

如果資料不是經緯度座標,或不想要背景地圖,只想畫某種等值線
gsn_csm_contour

想畫折線圖?
gsn_csm_xy

不要傻傻的作筆記把上述的每種類型的函數抄下來,這邊只有隨機列出幾種,目的在於展現NCL對於這類型函數的命名邏輯,就像 Resources 一樣是有跡可循的。有發現吧?除非能過目不忘,否則與其把這些函數名稱抄下來背下來,不如學學如何查詢,要用的時候再請谷歌大神幫忙就好拉。





待續... Day 4

沒有留言:

張貼留言