Grads 使用手冊

全新改板中,預計 8 /12 前完成

歡迎拷貝,提供建議,心得或

在本手冊加入您的使用方法及心得

任何建議請e-maildpengcm@rainbow.atm.ncu.edu.tw 彭啟明

 

Index:

  1. Grads 的介紹
  2. Grads 的硬體需求及基本設定
  3. 基本指令
  4. 繪製範例1
  5. 繪製範例2
  6. 繪製範例3
  7. 繪製範例4
  8. 斜溫圖
  9. Grads 的網路資源
  10. 台灣地區使用 Grads 的人士

 

 

 

 

  1. Grads 的介紹

    當您在分析資料時,清晰的圖表,是表現成果最有效的方式,雖然市面上已有許多相當方便的科學繪圖軟體,但對大氣科學上而言,其仍有所不足之處,例如等值線很難看、風標及流線場無法畫、通用性及價錢等考量。目前氣象上最通用的就屬NCAR GRAPHIC,其雖可滿足我們的需求,但每次都要寫程式,又有一堆堆的設定,對一般的"懶人"及要求時效的人來說,是相當麻煩的。而Grads 是一套專門為氣象方面所用的繪圖軟體,優點為:

    1. 簡單,只要給一定的BINARY DATA,不須寫程式,只需線上下指令即可。
    2. 方便,目前在各種平台工作站如 DEC, IBM, HP, Linux 等都可支援(X-WINDOWS)PCDOS亦可,目前亦準備有Windows的版本。
    3. 易學,您也許不相信,看這demo,把您的常用資料放進來試試看,一小時以內即可畫出來。
    4. 經濟,目前是freeware

    雖有以上的便利,可能是過於便利,因此缺少了一些圖形的細節控制,例如許多都已是default,要改很難,筆者建議您若是做研究分析,時間就是一種競爭,用Grads是很好的工具來協助您,但是若成果要發表成論文時時,NCAR GRAPHIC還是較好的選擇,除非你不大計較這些東西。

     

  2. Grads 的硬體需求

    目前各工作站及PC幾乎皆可支援,詳情可見

    http://grads.iges.org/grads

    或中央大學大氣系軟體資料庫

    ftp://www.atm.ncu.edu.tw/grads/

    內有最新訊息及hypertext,ps,text,doc檔的說明,約150頁。(為節省資源,用WWW看就好)

    其實根本不需要印,只要看看 grads reference card version就可了,

    ftp://www.atm.ncu.edu.tw/grads/manual/

     

    ftp: grads.iges.org

    您要最新的版本可在此發現。

     

     

  3. 基本指令及基本概念

    開始學 GRADS 可分為六個東西

    1. 資料處理程式: 將資料格式轉成 BINARY檔 。
    2. OPEN(8,FILE='samp.dat',FORM='UNFORMATTED',ACCESS='DIRECT', RECL=100)

      WRITE (8,REC=1) X 

    3. Binary : 可為 *.bin
    4. 描述 data 格式的檔: *.ctl
    5. 繪圖控制檔: *.gs
    6. 轉圖成為 postscript。(gxps, gxpsc, gxpscw)
    7. 印 postscript圖。(用postscript印表機或用ghostview)

     

     

    若您在PC下請在autoexec.bat內加入兩個路徑

    path .........;c:\grads

    append ......;c:\grads

     

    若您使用 PC 版本,妳可利用網路芳鄰連至 ATMOS 網域內,選擇 SOFTWARE-FREE當成妳的 H 碟,此時在 H:\下鍵入 RUN 此時包括 FORTRAN 設定及 GRADS 設定皆已完成。

     

    執行grads就可以進入 Grads繪圖模式內,出現Grads的提示符號鍵入

    ga>run maps.gs

    你就可以看見一些DEMO的圖,您要看更多可打test即可(PC版才有)

     

    您在IBM上以X-windows內打grads即可。 

     

     

     

    GRADS 繪圖樣式有兩種選擇,landscape(-A4) portait(A4)在打入GRADS 後即要妳選擇繪圖樣式,妳也可加入一些指令簡化流程,

    -b run batch mode 此時不會秀出任何圖形

    -l run in landscape mode

    -p run in portrait mode

    -c 執行後者所提供的 command 指令,例如 *.gs

    ex:

    grads -c "run batch.gs"

    grads -blc "run batch.gs"

     

    原始範例

     

      我們看看 maps.gs 內容,'reinit'reinit是對系統的設定值重新定義,而 ’ ’內則表示這是在Grads繪圖狀態下的指令,相對於若你不以run maps.gs的方式進行,直接在 ga>reinit 一行一行鍵入也是可以。接著在maps.gs內開了三個檔案

      'open sa.ctl'

      'open ob.ctl'

      'open wx.ctl'

      在*.ctl 檔內包含了二個很重要的部分,一是定義圖形的格式,另一個是定義所要繪出的變數,我們來看 sa.clt 內容:

    DSET sa.dat

    TITLE 0.5 degree analyzed grids

    UNDEF -9.9E33

    XDEF 81 LINEAR -110 0.5

    YDEF 51 LINEAR 25 0.5

    ZDEF 1 linear 1 1

    TDEF 1 LINEAR 16z16feb1993 1hr

    VARS 5

    slp 0 99 sea level pressure

    ts 0 99 surface temp

    ds 0 99 surface dwpts

    us 0 99 surface u winds

    vs 0 99 surface v winds

    ENDVARS

     

    一開始的 DSET sa.dat 是定義你所要用的資料檔檔名,

     

    VARS 5

    slp 0 99 sea level pressure

    ts 0 99 surface temp

    ds 0 99 surface dwpts

    us 0 99 surface u winds

    vs 0 99 surface v winds

    ENDVARS

     

    表示有五個變數,內含的就是定義所有變數的名稱和格式,例slp 0 99 sea level pressure 其中099表示時間和空間都只有一層,而 sea level pressure是變數的說明(若時間和空間不止一層則可用相關loop 指令做成動畫,在此不做介紹)。實際上應寫成

    Time 1,Level ? ,Variable slp

    即第一個是代表時間,第二個是代表層數(編者不知99的意義,Grads.doc內沒有說明詳細,大概是預留的未知個數)

      TITLE這僅是對此檔案的說明,幫助了解這一個檔案的作用。

      UNDEF無效值或可忽略的資料值,在此為-9.9E33

      XDEF 表示經度上X軸的點數,在此有81個點,起始值由-110開始,間值是0.5,增加的方式是線性(Linear)。若不是線性,則用 LEVELS -110,-100,80,......,逐一列出。

      YDEF 緯度上......和上面相同。

      ZDEF 高度上......和上面相同。

      TDEF 時間上......和上面相同。但時間表示的格式是 hh:mmZddmmmyyyy

          hh:hour

          mm:minutes

          dd:day

          mmm:month

          yyyy:year

    例如:12Z1JAN1990,14:20Z22JAN1987等。

     

    做完*.ctl後,如何把資料變數放進去,把變數場繪出來。底下先說明一維陣列的方式

    REAL X(100)

    DO 10 I=1,100

    X(I)=I

    10 CONTINUE

    OPEN(8,FILE='samp.dat',FORM='UNFORMATTED',ACCESS='DIRECT', RECL=100)

    WRITE (8,REC=1) X

    STOP

    END

     

    必須注意的是RECL,在PC中實數是4byes因此必須乘上4,才是所須要的大小。底下則是一個三維矩陣

    三維矩陣的排列方式:

    REAL Z(72,46,16)

    OPEN (8,FILE='grads.dat',FORM='UNFORMATTED', ACCESS='DIRECT',RECL=72*46)

    IREC=1

    DO 10 I=1,16

    WRITE(8,REC=IREC) ((Z(J,K,I),J=1,72),K=1,46)

    IREC=IREC+1

    10 CONTINUE

     

     

    ob.ctl,在Grads內對於測站圖形的處理相當的美觀,因此如何將測站的資料融入圖形內是很重要,對此Grads自有一套格式來解決。

      dtype station

      stnmap ob.map

     

     

    Year Month Stid Lat Lon Rainfall

     

    1980 1 QQQ 34.3 -85.5 123.3

    1980 1 RRR 44.2 -84.5 87.1

    1980 1 SSS 22.4 -83.5 412.8

    1980 1 TTT 33.4 -82.5 23.3

    1980 2 QQQ 34.3 -85.5 145.1

    1980 2 RRR 44.2 -84.5 871.4

    1980 2 SSS 22.4 -83.5 223.1

    1980 2 TTT 33.4 -82.5 45.5

    .

    .

    .

    A FORTRAN program in DEC FORTRAN to write this data set

    in GrADS format might be:

     

     

    CHARACTER*8 STID

    C

    OPEN (8,NAME='rain.ch')

    OPEN (10,NAME='rain.dat',FORM='UNFORMATTED',

    & RECORDTYPE='STREAM')

    C

    IFLAG = 0

    C

    C Read and Write

    C

    10 READ (8,9000,END=90) IYEAR,IMONTH,STID,RLAT,RLON,RVAL

    9000 FORMAT (I4,3X,I2,2X,A8,3F8.1)

    IF (IFLAG.EQ.0) THEN

    IFLAG = 1

    IYROLD = IYEAR

    IMNOLD = IMONTH

    ENDIF

    C

    C If new time group, write time group terminator.

    C Assuming no empty time groups.

    C

    IF (IYROLD.NE.IYEAR.OR.IMNOLD.NE.IMONTH) THEN

    NLEV = 0

    WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG

    ENDIF

    IYROLD = IYEAR

    IMNOLD = IMONTH

    C

    C Write this report

    C

    TIM = 0.0

    NLEV = 1

    NFLAG = 1

    WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG

    WRITE (10) RVAL

    GO TO 10

    C

    C On end of file write last time group terminator.

    C

    90 CONTINUE

    NLEV = 0

    WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG

    STOP

    END

     

     

     

      另外一個檔案wx.ctl 則依此類推。

     

     

    常用功能摘要(請見附表)

    Graphic Type

    set gxout bar BAR

    set gxout winbarb winbarb

    set gxout contour 畫等值線

    set gxout shaded 陰影

    set gxout vector 風向量圖

    set gxout sream 氣流線

    set gxout line 直條圖

    set gxout fgrid Shaded grid boxes

     

    Graphic options

     

    set csmooth on|off 是否用cubic interpolation內差

    set cmin 0 最小值

    set cmax 20 最大值

    set cint 5 距隔值

    set clevs lev1 lev2 .... 等值線值

    Page control 圖形大小控制

    Map projections/drawing 地圖投影

    'set mproj latlon' 定義值,經緯度直角座標。

    scaled 緯度投影。

    nps north ploar stereographic

    sps south polar sterographic

    off 不做地圖投影。

    'set mpvals -105 -75 25 50' 設定繪出座標的範圍,在此表示105W~75W,25N~50N

    'set mpdset lowers' 內定值,

    mres 表示有州界

    hires 表示含國界

    nmaps

    'set poli on|off 內定值是on ,表示上一項有否要作用。

     

     

     

  4. 繪製範例 1

    筆者範例 ftp://www.atm.ncu.edu.tw/grads/example/1

    以上是照原文翻過來的方法,在此筆者用一個自己嘗試的範例來說明,筆者的資料是環保署桃園站一個月的資料,目標是看其一個月的每日變化,筆者用上述 z(i,j,k)的寫法將原始環保署data成了一個 hs9417.bin檔因此描述此檔的file可寫為 epa.ctl

    DSET hs9417.bin

    TITLE model data

    UNDEF -999.0

    XDEF 24 LINEAR 0 24 -----> x軸為24小時,24

    YDEF 31 LINEAR 0 31 -----> y軸為31日,31

    ZDEF 1 linear 1 1 -----> 沒有z故以1 表示

    TDEF 1 LINEAR 12z4feb1994 12hr

    VARS 16 -----> 16個變數

    so2 0 99 precipitable water ----->依順序為16個變數,內容懶得寫故都一樣

    co 0 99 lifted index

    o3 0 99 sea level pressure

    pm 0 99 precip

    nox 0 99 z1000

    no 0 99 z850

    no2 0 99 t850

    thc 0 99 u850

    nmhc 0 99 v850

    wspeed 0 99 rh850

    wdir 0 99 z500

    U 0 99 z500 ----> 畫風標其只認得分量

    V 0 99 z500

    temp 0 99 u500

    rainfall 0 99 rain

    ch4 0 99 ch4

    ENDVARS

     

    接下來可畫圖了,進入grads

    ga>open plot.ctl

    您可看到 lon .....,lat.....等,因為其會將其當成經緯度來看,我們先不管,座標軸示等一下說明。可以先用query file來顯示plot.ctl內宣告的內容,若此時要畫 O3的月內日變化圖可直接 ga>d o3,此即可繪出,但如前所述,其會將其當成經緯度來看,故座標會怪怪的,因此可ga> set xaxis 1 24 1,表x軸圖示為124間隔1,同時還有世界地圖在內,故須 set mpdraw off,再將網格去掉 set grid offclear一下,再d 03就可畫出一張不錯的圖,若要看每日中午12時的 O3日變化,只需set x 12,再 d O3 即可看到一個剖面圖。當要畫好幾個圖時,這樣太麻煩了,可以寫成批次檔如下所示

     

    epa.gs

     

    'reinit'

    'clear'

    'open plot.ctl'

    'enable print test' ---->將圖形初出至test

    'set xaxis 0 24 1'

    'set yaxis 0 31 1'

    'set mpdraw off'

    'set csmooth on' -----> smooth

    'set grid off'

    'draw title EPA 17 O3'

    'draw xlab Hour'

    'draw ylab Day(May)'

    'set t 1'

    'set z 1'

    'd o3'

    'print' ----> 印在test

    pull dummy ----> 有點像call wait

    'clear'

    -----> 下一張

    'set xaxis 0 24 1'

    'set yaxis 0 31 1'

    'set mpdraw off'

    'set csmooth on'

    'set grid off'

    'draw title EPA 17 O3'

    'draw xlab Hour'

    'draw ylab Day(May)'

    'set t 1'

    'set z 1'

    'd u;v' -----> 風的 vector

    'print'

    pull dummy

    'clear'

    'disable print' -----> 結束輸出圖形

     

    故只要 ga> run plot.gs即可做完一堆工作,此時會產生test,將其用gxps轉成postscript檔,在pc上用ghostview,工作站用xv ghostview即可看出。

     

  5. 繪製範例2

    ftp://www.atm.ncu.edu.tw/grads/example/2

    這是幫我老闆投到 JGR 的PAPER所弄的,妳可去最近的1997 PAPER上找到,特殊的是為非等間距,無效值為陰影區及三張圖疊在一起,不過要加一些特殊符號還是很麻煩的,所以GRADS還是有其侷限性!

     

     

  6. 繪製範例3

    ftp://www.atm.ncu.edu.tw/software/grads/example/3

    由於我不是做大尺度的,這是謝榮傑所提供的,主要是在 ibm與 dec機器在存取binary 時之差異,妳可傳至 IBM試試

    在 IBM, HP, SUN 的位元數值,是採取 4BYTES=32BIT的存取位元資料,是IEEE的標準,其是屬於 big endian,若是在 ms-dos下開資料時用 *4則可適用於 big endian.

    在 DEC, PC linux則是 little endian,與big endian 不同,但在dec 上,也可讀 big endian code,只要在 open 這個敘述中加入 convert='big_endian' 即可

    open(10,file='10.bin',status='old',form='unformatted',convert='big_endian')

     

     

  7. 繪製範例4

    ftp://www.atm.ncu.edu.tw/grads/example/4

    這是東亞地區硫化物排放之模擬。

     

  8. 斜溫圖

    ftp://www.atm.ncu.edu.tw/grads/example/skew

    斜溫圖,是氣象很獨特的一種圖形,一般而言,目前除了 NCAR GraphicPC上的繪圖軟體,一般是相當難繪製的。但皇天不負苦心人,賓州大學氣象系研究生Bob Hart (Graduate Student, Penn State University Meteorology, http://www.ems.psu.edu/~hart/skew.html )終於將 Grads 版斜溫圖寫出來,並公開予大眾使用。

    目前經測試,其函數及各項數值均為正確,但目前只能在工作站上 runpc部份仍有些問題待解決。

    關於斜溫圖的程式範例詳見 ftp://www.atm.ncu.edu.tw/grads/example/skew

    須注意的是風速可以設一個很小值,絕不能為0,這是許多錯誤嘗試的經驗!

     

     

  9. Grads 的網路資源

   GrADS Users Group List Server

   email to:

   listserv@icineca.cineca.it

   To subscribe, send mail to that address containing the following line:

   subscribe gradsusr #####