背景
- 此一腳本將在一頁A4版面內繪製8個測站的濃度時序圖，並集合所有頁面在同一pdf檔案內。
- 各污染項目的製圖腳本是以模版形式儲存。
- 模版內可更動的變數為第三行的spe：污染物名稱
- 可能的污染物名稱：包括在epa.ctl檔案中（
so2 cmo ozn pmt nox p25 no2 voc）
- 控制修改模版與執行grads的(上層)程式：sss.cs
- NMHC只有部分測站有觀測值。作圖時仍然保留沒有測值之測站模擬值。
- CMAQ類似的功能詳見[[2022-07-21-wsite]]1
- GrADS安裝及資源詳參[[2022-07-21-grads]]2
分段說明
讀取epa.ctl檔案
- 由所在目錄下的epa.ct檔案中讀取小時頻率(
HHH)，其值可能是1或3小時。視資料量而定。 - 依序讀取epa.ct中的變數定義段落(第7行以後)：
- 每污染物有2行，分別為觀測(o)、或模擬(m)。
- 如果變數名稱與要求繪製的內容一致，則記下該變數名稱(
spnm)及單位(spnun)，以便寫在畫面上。
1 'reinit'
2 'clear'
3 SPE=oz8
4 "open epa.ctl"
5 "q file 1"
6 rec = sublin(result,5)
7 HHH = subwrd(rec,3)
8 i=1
9 while(i<=8)
10 ii=(i*2-1)+6
11 rec = sublin(result,ii)
12 sso = subwrd(rec,1)
13 sss = substr(sso,1,3)
14 if(sss=SPE)
15 spnm=subwrd(rec,4)
16 spun=subwrd(rec,5)
17 endif
18 i=i+1
19 endwhile
迴圈之控制
- 控制字的大小：20~22、顏色(24)
- x/y/t的初始值與間距(前述
HHH) - 頁面序號：
ipage、共7頁。測站序號：is，共每頁10站，最多65站(第7頁只有5站)。
20 'set string 1 bc'
21 'set strsiz 0.15'
22 'set csmooth on'
23 'set grads off'
24 'set ccolor 2'
25 "set x 1"
26 "set y 1"
27 "set t 1 "HHH
28 ipage=1
29 while (ipage <=7)
30 is=(ipage-1)*10+1
31 isb=is
32 isend=isb+9
33 if (isend > 65) ;isend=65;endif
34 while (is <= isend)
35 "set digsiz 0.05"
- 由副程式找到對應的測站編號及站名(羅馬拼音)
- GrADS的程式語言並沒有序列或矩陣、DICT等功能，因此如果要類似功能，只能使用副程式
- 測站圖面控制
- 因為是垂直頁面，尺寸為8.5吋×11吋。
vv為控制每站繪圖垂直位置的原點。每站圖的高度為2吋。
36 ist=num_of_st(is)
37 jst=ist-1+1
38 ast=nam_of_st(jst)
39 iss=is-isb+1
40 vv=10-iss*1
41 tt=vv+2.0
42 *"set parea 0 8.5 0 11"
43 "set vpage 0. 8 "vv" "tt
44 "set grads off"
- 由數據檔統計得到時間範圍的極值，做為濃度刻度的計算
- 觀測值的極值(
ymx,ymn) - 模擬值的極值(
zmx,zmn) - 取這4值的範圍(
ymx,ymn)
- 觀測值的極值(
45 'set gxout stat'
46 "d "SPE"o(stid="ist")"
47 data = sublin(result,8)
48 ymx = subwrd(data,5)
49 ymn = subwrd(data,4)
50 "d "SPE"m(stid="ist")"
51 data = sublin(result,8)
52 zmx = subwrd(data,5)
53 zmn = subwrd(data,4)
54 if (zmx > ymx) ; ymx = zmx ; endif
55 if (zmn < ymn) ; ymn = zmn ; endif
- 最大最小間至少要有4個刻度
- 歷線型態
- 觀測值(黃點：
line 1 1 6) - 模擬值(紅點線：
line 1 2 3)
- 觀測值(黃點：
- 在左測寫上測站編號及名稱
56 "set gxout line"
57 di1=(ymx-ymn)/4/10
58 di2=int(di1)
59 if (di2 < 1) ; di2 = 1 ; endif
60 div=di2*10
61 "set ylint "div
62 "set xlint 24"
63 'set vrange 'ymn' 'ymx
64 'set line 1 1 6'
65 'set strsiz 0.03 0.03'
66 "set cthick 15"
67 "d "SPE"o(stid="ist")"
68 'set line 1 2 3'
69 "d "SPE"m(stid="ist")"
70 'set strsiz 0.1 0.08'
71 "draw string 1.2 1.2 " ist
72 "draw string 1.2 1 " ast
73 is=is+1
74 endwhile
分頁圖名之填寫
- draw 污染名稱、單位、起迄站名
- 輸出單頁結果到png檔案
75 ast0=num_of_st(isb)
76 ast1=num_of_st(isend)
77 'set strsiz 0.15 0.1'
78 "draw string 4.5 0.2 "spnm" "spun" of Station " ast0 " to " ast1 "(o for obs)"
79 *pull dummy
80 "printim "ipage".png x1700 y2200 white"
81 clear
82 ipage=ipage+1
83 endwhile
84 quit
85 return
86 end
副程式
測站代碼與序號之對照
88 function num_of_st(is)
89 nst.1=001
90 nst.2=002
...
153 nst.65=084
154 ist=nst.is
155 return ist
測站名稱與序號之對照
156 function nam_of_st(is)
157 nst.1=KeeLong
158 nst.2=Xizhi
159 nst.3=WanLi
...
219 nst.83=MaiLiao
220 nst.84=FuGuiJiao
221 ast=nst.is
222 return ast
實數整數之交換
234 * Take a floating-point (decimal) number as a parameter and returns a truncated integer (with no decimal point)
235 function int(num)
236 outnum = ''
237 i = 1
238 while(i <= strlen(num))
239 char = substr(num,i,1)
240 if(char = ".")
241 break
242 else
243 outnum = outnum%char
244 i = i+1
245 endif
246 endwhile
247
248 return outnum
結果
 |
|---|
| 圖 GrADS執行spec.gs之結果範例 |