GrADS繪圖學(xué)習(xí)技巧與實(shí)例

1、實(shí)用標(biāo)準(zhǔn)文案以下技巧總結(jié)都是筆者從學(xué)習(xí)實(shí)踐過程中總結(jié)出來的，基本的問題。不求全面，希望對讀者學(xué)習(xí)有用，如果有問題，敬請留言指正，以促進(jìn)交流學(xué)習(xí)！1、軟件綜述：grads軟件是一款繪圖軟件除了繪制圖形，還可以提取數(shù)據(jù)，主 要應(yīng)用是在大氣科學(xué)中，當(dāng)然只要是數(shù)據(jù)處理成grads能夠讀取的數(shù)據(jù)文件就可以進(jìn)行相關(guān)繪圖。軟件版本問題，軟件本身不是很大，我接觸到18 1.9、 2.0版本的，1.8版本的安裝很多情況還要修改環(huán)境變量、1.9版本的不識別'sdfopen'命令，最穩(wěn)定的版本是2.0版本，所以筆者推薦學(xué)習(xí)者安裝2.0版本，選擇默認(rèn)安裝路徑就可以。2、文件類型簡述：grads處理的是

2、網(wǎng)格數(shù)據(jù)，可以處理的數(shù)據(jù)類型有：grd、grib、nc (海洋常用的數(shù)據(jù))，cdf (雷達(dá)衛(wèi)星數(shù)據(jù))，其中nc、cdf數(shù)據(jù)都是自 帶描述文件，不需要ctl, grib數(shù)據(jù)要通過命令生成ctk index數(shù)據(jù)才可以調(diào) 用，常用的是grd數(shù)據(jù)，需要ctl。3、數(shù)據(jù)文件轉(zhuǎn)換：grads軟件識別的數(shù)據(jù)是二進(jìn)制無格式數(shù)據(jù)，文件類型是'binary',寫入和生成時是不需要格式的如 read(20) sst(i,j,iz,it),20為文件號， 通常是十進(jìn)制數(shù)據(jù)與grd數(shù)據(jù)間轉(zhuǎn)換，這里給一個 grd轉(zhuǎn)換成txt數(shù)據(jù)的 fortran 程序：parameter(nx=56,ny=41,nz=1

3、,nt=360)dimension sst(nx,ny,nz,nt)real sstopen(15,file='sst.grd',form='binary') !固定的用 form= binary'就是二進(jìn)制數(shù)據(jù) open(16,file='sst.txt')!新建 txt 文件do it=1,ntdo iz=1,nzread(15) (sst(i,j,iz,it),i=1,nx),j=1,ny) ! read 后只有文件號，數(shù)據(jù)是無格式的 enddoenddodo it=1,ntdo iz=1,nzwrite(16,*) (sst(i,

4、j,iz,it),i=1,nx),j=1,ny) !輸出日中是 txt 文件可直接看的數(shù)據(jù)， 有格式輸出，有 *enddoenddoclose(15)close(16)end寫程序時：注意格點(diǎn)數(shù)要與數(shù)據(jù)對應(yīng)，如：上程序?qū)?yīng)的數(shù)據(jù)是 經(jīng)度90200, 緯度-2060,時間：1971.012000.12共360個月的海面溫度數(shù)據(jù)，數(shù)據(jù)格點(diǎn)精 度 2*2 , nx=(200-90)/2+1,ny=(60-(-20)/2+1,nt=360, nz=1,大氣的數(shù)據(jù)要根據(jù) 數(shù)據(jù)的層次確定幾層。4、 grd、ctl、gs nc 詳述grd文件：grd數(shù)據(jù)不可直接看，為二進(jìn)制無格式數(shù)據(jù)，簡單的說只有1和0,而

5、 且數(shù)據(jù)間沒有間隔，grads識別grd是根據(jù)ctl進(jìn)行劃分白1根據(jù)ctl中的經(jīng)度、 緯度、層次、時間，精度進(jìn)行數(shù)據(jù)分塊。ctl實(shí)例：dset C:datasst.grdundef -9.99E+33title sea surface tempturexdef 56 linear 90 2ydef 41 linear -20 2tdef 360 linear jan1971 1mozdef 1 levels 0vars 1sst 0 99 surface sea tempture endvars上例數(shù)據(jù)的數(shù)據(jù)順序是（以下是數(shù)據(jù)對應(yīng)的經(jīng)緯度）緯度經(jīng)度（t=1）- 2090200的 緯度20S的從

6、90E200的56個數(shù)據(jù)- 1890200的 緯度18S的從90E200的56個數(shù)據(jù)6090200的 緯度60N的從90E200的56個數(shù)據(jù)以上為一層的數(shù)據(jù)，接下來是t=2, t=3t=360的數(shù)據(jù)，每個時間點(diǎn)的每一層 是如上格式，編程時讀取和寫如的數(shù)據(jù)循環(huán)順序依次是：時間、層次、緯度、經(jīng)度，讀者參照3中的fortran程序加以理解。ctl文件：具體其他指導(dǎo)書上都有，我這里強(qiáng)調(diào)的是 sst后面的0表示一層，如 果是兩層以上則是2,3，1層是固定用法，sst后面的99是默認(rèn)設(shè)置；undef- 9.99E+33此處的數(shù)值決定了軟件將文件中的那些值 認(rèn)定為不繪制的數(shù)據(jù)，所以 這個值一定要與數(shù)據(jù)對應(yīng)。

7、gs文件，批命令文件，與ctl 一樣是用記事本編寫，另存為.gs文件，文件都是 命令，方便大段的命令編寫、修改，很常用，尤其時繪制的圖要求比較多時必 須用，免得在命令窗口重復(fù)輸入命令浪費(fèi)時間，如下例子：'reinit''open c:datauv.ctl''set t 7'set lev 850''set grid off''set vpage 0 8.5 0 3.6''set parea 0.3 8.5 0.3 3.5''set gxout vector''set

8、grads off''d u;v''set vpage 0 8.5 3.6 7.2''set grads off''set gxout barb''d u;v''set vpage 0 8.5 7.2 11''set grads off''set gxout grid''d u;v''printim c:imagesgxout3uv850.png white'Gs文件以；結(jié)尾，分號后面不能有空格，除了循環(huán)命令外都需要單引號將命令

9、引起來。nc文件：此類數(shù)據(jù)不需要ctl可直接用sdfopen命令直接打開，往往需要知道數(shù) 據(jù)文件中的各個緯度特征，可以用：q ctlinfo命令查詢該文件的內(nèi)置 ctl,這樣 一切都會很清楚，尤其是數(shù)據(jù)精度：截圖如下（下例即2*2的網(wǎng)格經(jīng)度）：sdfopen c：datasst.nnmean _ncSeann insr self-describing file： c ：st. Rnme an«n cSDF file c ： Xdatasst. mnmean -nc is open as file 1LOH set to 0 36日 LAT set to -88 88 LEU set

10、to O 0Tine ua.lues set： 1854:1 ：1 -9 1854:1 ：1 ：0 E set toll q ctlinfodset c ； datasst. mnmean ,nctitle NOAA Extended Reconstructed SSI U3 Lin def -9 .99e+33 dtype netedf xdeF 180 lineup 8 2 ydef 8? linear -88 2 zdef 1 linear 0 1 tdef 1869 linear 001JAH1854 imo uai*s 2tine_bnds =>t iu»e_bncl

11、s 0 t I ine Boundar-iessst,=>sst 0 t, y,x Montliil Means of Sea Surf ace Tenperature endvarsSR-以上為基礎(chǔ)知識介紹，以下為技巧命令1、繪圖時，圖的時間下標(biāo)可用set grads off命令關(guān)掉，網(wǎng)格用set grid off命令 關(guān)掉，需注意的是，網(wǎng)格關(guān)一次就一直有效，而下標(biāo)則是每次繪完圖就自動開啟，所以建議讀者在每次的 繪圖命令之前加set grads off。2、3、reinit命令是讓窗口恢復(fù)到剛打開時的界面，會恢復(fù)所有的sei:,同時關(guān)掉所有的打開文件，為為防止前面打開文件的干擾，建議讀

12、者在每個gs文件的第 一條就加上這條命令(參考gs描述的例子)。4、5、數(shù)據(jù)維度設(shè)定：set lon 90 set x 1set lat 80 set y 20set lev 1000 set z 1set t 1 12 set time jan1970 dec1970上面兩種設(shè)定等價，左邊是實(shí)際維度設(shè)定法，右邊是給點(diǎn)設(shè)定法，讀者需 要確定具體的格點(diǎn)數(shù)。需要注意的是，在設(shè)定全球尺度時，經(jīng)度0和360是同一個格點(diǎn)，所以set 0 360會出錯，這時設(shè)定格點(diǎn)的方法比較好：set x 1 180。6、 vpage和 parea 的區(qū)另U:vapge是對整個繪圖區(qū)分塊，需要幾張圖就劃分成幾塊，給英寸時

13、不需要給 標(biāo)注、標(biāo)題預(yù)留空間，如：要橫著繪制兩張圖 那么就是set vpage 0 5.5 0 8.5(11*8.5)繪制第一張圖set vpage 5.5 11 0 8.5繪制第二張圖parea是描述的繪圖時圖形的四根邊線的大小，而且是 虛頁的尺寸如：畫一張圖 set parea 1.0 10.2 0.8 7.8四個值的范圍 取決于vapge的長 度0v1.0v10.2v11,如果是上例中的兩張圖，那么： set parea 0.5 4.8 0.8 7.90<0.5<4.8<5.5 ,數(shù)值范圍取決于水平、垂直的長度，與起點(diǎn)無關(guān)如 set vpage 0 11 3.5 7se

14、t parea 0.8 10 0.3 3.2y上滿足0<0.3<3.2<(7-3.5)即可，實(shí)際y長度決定范圍，與起點(diǎn)無關(guān)。5、'set annot 5 8''set xlopts 3 5 0.18''set ylopts 3 5 0.18'上面這兩條命令是設(shè)定的x、y軸下標(biāo)數(shù)字的顏色、粗細(xì)、字號'set annot 5 8'是設(shè)定坐標(biāo)軸線、標(biāo)題的顏色粗細(xì)，會重置xlopts的部分設(shè) 定，希望讀者注意，命令間的互相干擾6、cbar cbarn 命令 cabr03&912151321242730這里只講解c

15、barn命令，cbarn sf vert xmid ymid其中sf為標(biāo)尺,1為全尺寸0.5為半尺寸；vert為放置位置，0為水平，1為 垂直；xmind、ymid為色標(biāo)的中心位置如：cbarn 0.5 0 2.5 1.5色標(biāo) 半尺度長，水平放置，色標(biāo)的中心英寸坐標(biāo)(2.5,1.5)。如果只是cabrn或cbar命令，會按照默認(rèn)的全長，水平或數(shù)值取決于圖的哪 邊空位大繪制。強(qiáng)調(diào)一點(diǎn)，加色標(biāo)和加標(biāo)題都是在繪制出圖形后才能加 ，因?yàn)樯珮?biāo)是根據(jù)陰影 圖確定對應(yīng)色值，標(biāo)題根據(jù)圖形大小確定標(biāo)題位置。7、cmin、cmax是命令是用于繪制大于或小于某數(shù)值的線或區(qū)域圖形，在每次 繪圖之后會重設(shè)，所以如果多次

16、使用一定要每次 d之前加上，這點(diǎn)與set grads off相似。8、求12個月每個月的海溫距平值 'set t 1 12''asst=ave(sst,t+1404,t=1764,12)''modify asst seasonal' 'set time jan1971' 'd asst''set t 1405 1764''nasst=sst-asst''set time JAN1998' 'd nasst'set t 1 12是設(shè)定asst變量有12個

17、時間序列，每個時間格點(diǎn)放一個平面的平 均值，asst=ave(sst,t+1404,t=1764,12淀義變量asst放每個月的平均值，隨著 t從1 變化到12, 一次求的每個月的平均值，起始時間是t=1405,終止時間是1764,t=1時，相當(dāng)于 1+13+25+37，即每年的 1月份的值求平均(這里t=1+1404=1405是1971年1月，因此是求19712000年每年1月份的平均值), t=2,3.4,12與1同理。如此將12個月的平均值都放在了 asst里。Modify asst seasonal本來12個月的平均值知識放在了 112的時間序列里， (這里調(diào)用的文件數(shù)據(jù)開始時間是 1

18、854年1月)，所以asst的值只是在設(shè)定112 以內(nèi)才能畫出，假了這條命令，可以將 asst的時間序列擴(kuò)展到所有時間里，使 得每一年的每個月對應(yīng)都是該月30年的平均值，方便后面求距平。需要注意的，如果set t 1 4 ,那么這條命令的作用是每年的 asst的14月值是一樣的，是該 月的平均值。由于grads本身軟件有一定的問題，當(dāng)你不是設(shè)定一段時間而是一 個時間點(diǎn)如：sett 1,步、不用modify命令，所有的時間序列都會有該值'sett 1405 1764''nasst=sst-asst'這兩句是求30年（1971 012000 12每個月的海溫距平值，

19、共360各月的距 平值。這個不難理解，不做解釋。以下為實(shí)例：1、利用所提供的數(shù)據(jù)文件，繪制出 2003年7月60-150E、0-40N區(qū)域內(nèi)700hPa 流線圖，且地形高度場超過2000米以上用黑色陰影顯示，并且給出相應(yīng)標(biāo)題（請 包含姓名拼音與學(xué)號），最終將圖形保存。所有命令編寫于.gs文件中。'reinit''open c:datadxgd.ctl''open c:datauv.ctl''set grid off'set grads off'set map 3 1 5''set xlopts 3 5 0.

20、18''set ylopts 3 5 0.18''set parea 0.7 10.4 0.6 7.8''set t 1''set lon 60 150''set lat 0 40''set gxout shaded''set rbcols 1 1 1 1 1 1 1''set cmin 2000''d h''set gxout stream''set time jul2003''set grads of

21、f'set lev 700''d u.2;v.2''set annot 5 8''draw title DingXiaoli 20081331001''printim c:imagesuv200307850.png white'2、利用所提供的數(shù)據(jù)文件，畫出 2002年1-12月120 E、0-40N 200hPa緯向風(fēng) 的緯度-時間剖面圖。圖型要求：(1)緯向風(fēng)為西風(fēng)時填色，東風(fēng)繪制等值線，給出色標(biāo)， 0值線加粗；(2) X軸標(biāo)注為“time ”，Y軸標(biāo)注為“ lat ”，標(biāo)題標(biāo)注為“ u 100-120E”。(

22、3)最終將圖形保存為gmf格式。所有命令編寫于.gs文件中'reinit''open c:datauv.ctl''set grid off'set grads off'set xlopts 3 5 0.16''set ylopts 3 5 0.16''set parea 1 10.4 1.2 7.9''set lon 120''set lat 0 40''set lev 200''set t 1 12''set annot 5

23、8''enable print c:imagesuv200.gmf'set xyrev on''set gxout shaded''set cmin 0''d u''cbarn 1 0 5.5 0.3''set grads off'set gxout contour''set cmax 0''set grads off'd u''set clevs 0''set ccolor 2''set cth

24、ick 10''set grads off'd u''draw xlab time''draw ylab lat''draw title u 100-120E''print''disable print'3、十二個月的nc數(shù)據(jù)整合到一個文件中(fortran)這個程序可用于整合數(shù)據(jù) parameter(nx=360,ny=181,nz=26,nt=12)dimension temp(nx,ny,nz,nt)open(1,file='argo_200901.grd',f

25、orm='binary')open(2,file='argo_200902.grd',form='binary')open(3,file='argo_200903.grd',form='binary')open(4,file='argo_200904.grd',form='binary')open(5,file='argo_200905.grd',form='binary')open(6,file='argo_200906.grd',f

26、orm='binary')open(7,file='argo_200907.grd',form='binary')open(8,file='argo_200908.grd',form='binary')open(9,file='argo_200909.grd',form='binary') open(10,file='argo_200910.grd',form='binary') open(11,file='argo_200911.grd

27、9;,form='binary') open(12,file='argo_200912.grd',form='binary') open(37,file='argo_2009.grd',form='binary') do iz=1,nzread(1) (temp(i,j,iz,1),i=1,nx),j=1,ny)read(2) (temp(i,j,iz,2),i=1,nx),j=1,ny)read(3) (temp(i,j,iz,3),i=1,nx),j=1,ny)read(4) (temp(i,j,iz,4),i

28、=1,nx),j=1,ny)read(5) (temp(i,j,iz,5),i=1,nx),j=1,ny)read(6) (temp(i,j,iz,6),i=1,nx),j=1,ny)read(7) (temp(i,j,iz,7),i=1,nx),j=1,ny)read(8) (temp(i,j,iz,8),i=1,nx),j=1,ny)read(9) (temp(i,j,iz,9),i=1,nx),j=1,ny)read(10) (temp(i,j,iz,10),i=1,nx),j=1,ny)read(11) (temp(i,j,iz,11),i=1,nx),j=1,ny)read(12) (temp(i,j,iz,12),i=1,nx),j=1,ny)enddodo it=1,ntdo iz=1,nzwrite(37) (temp(i,j,iz,it),i=1,nx),j=1,ny) enddoenddoclose(1)close(2)close(3)close(4)close(5)close(6)close(7)close(37)end4、豎著繪制三張不同輸出形式的uv風(fēng)場圖(parea可保證圖形達(dá)到盡可能的大)'reinit''open c:da