Dr. Kuang's Utilities

測站模擬比較之時序圖

 
  • Jul 21, 2022
  • views
>

背景

  • 測站模擬比較之時序圖是個簡單的XY-plot，用excel就可以輕鬆產出。何以會變得複雜，主要是因為污染物、空品站、月份等維度所造成。
  • 環保署公版模式提供有python版本的後製工具，此處介紹GrADS版本的作業方式。
  • 公版的後製工具可以產生個別測站項目之圖檔， 然而對撰寫模擬報告的作者而言，產生個別比較檔案還算容易，困難的是如將一個個的圖檔整併在頁面上，排列整齊、充分使用版面、還能保持各圖仍有一定的清晰度。這也是ovm2grads的工作重點。
  • ovm2grads的整體工作如下圖所示。
    • CAMx的後處理詳流程圖
    • GrADS安裝及資源詳參grads(also[[2022-07-21-grads]])
  • 其他CAMx後處理整體之作業方式可以參考早期版本的github wiki說明，此處著重時序圖的製作與示範。

ovm2grads整體作業流程

graph TD
    A(ovm.dat) --> G((ovm2gr.cs))
    G --> B(stn2grads.dat,ovm2grds.result, epa.ctl, taiwan.map)
    D(spec.gs.BLK) --> H((sss.cs))
    H --> C(spec.gs)
    B --> E((grads))
    C --> E
    B --> H
    H --> E
    E --> F(*.pdf)
  • ovm.dat為aok(for CAMx)或wsite(for CMAQ)的執行成果，為一個龐大的對照表格，格式為文字檔，必須先轉成GrADS可以讀取的二進位檔stn2grads.dat(little endian)。
  • ovm2grds系列程式有2支，先寫的stn2grads.f適用在氣象數據的比對，而ovm2grads.f則適用在空氣品質的比較。
  • ovm2gr.cs、sss.cs詳下述。

執行成果及檢討

a3
圖 GrADS執行spec.gs之結果範例

ovm2gr.cs作業流程

graph TD
    A(cwb_epa.csv) -->  D((ovm2grads.f)) 
    B(ovm.dat) --> D
    D --> C(stn2grads.dat)
    D --> E(ovm2grds.result)
    F(epa.ctl.BLK) -- sed --> G(epa.ctl)
    E -- cat, tail and awk --> G
    G --> H((stnmap))
    C --> H
    H --> I(taiwan.map)

ovm2grads.f程式說明

程式設計

  • GrADS平面資料檔與點狀資料檔的格式有很大的差異。主要在於存取的方式，
    • 平面資料檔為直接存取，一筆記錄的大小(recl長度)即為二維矩陣的大小：
      OPEN(12,file='...',action='write',form='unformatted',access='direct',    
     & recl=NOXG*NOYG,status='unknown')
  • 點狀數據則為循序存取：
      open(21,file='stn2grads.dat',form='unformatted',    
     &status='unknown' ,access='stream')
  • 此外在資料內容上也有很大的差異:
      nlev=1    
      nflag=1    
      Atime=0    
      itime=0    
      ifreq=3    
      DO it=1,ntime,ifreq    
      DO JS=1,NST    
      if(lstn(js).eq.0)cycle    
      write(stn,'(I3.3,A4)')isq(JS)!,A4(JS)    
      !if(it.eq.1)print*,stn    
      write(21)stn,y(JS),x(JS),Atime,nlev,nflag,((var(i,j,js,it),i=1,2),j=1,8)    
      ENDDO ! LOOP FOR STATION    
      write(21)stn,0.0,0.0,0.0,0,0    
      itime=itime+1    
      ENDDO !LOOP FOR DAY    
99    print*,itime, ifreq
  • 第三，端序的設定
    • 雖然GrADS可以在.ctl檔案內指定二進位檔是大端序檔案(byteswapped)，或是小端序檔案(內設)，
    • 平面資料檔avrg2grads即在程式中將byteswapped寫入.ctl檔案，不必特別轉換。
    • 然而在處理測站測值資料檔過程中，因必須使用stmap程式處理測站座標，因此資料檔必須以小端序方式儲存，ovm2grads程式不可以用ifort -convert big_endian編譯。

sss.cs流程

  • 腳本將會置換spec.gs模版中的污染物名稱、執行grads、一頁10個測站、共7頁檔案，並將分頁結果合併成1個完整圖檔。
graph TD
    J($PWD) -- date, cut -->D{loop for 8 spe's}
    D -- yes --> A(spec.gs.BLK)
    A -- cp & sed --> B(spec.gs)    
    B --> E((grads))
    C(stn2grads.dat) --> E
    I(taiwan.map) --> E
    H(epa.ctl) --> E
    E -- convert and combine 7 png pages --> F( 1 pdf file)
    F -- next spe --> D
    D -- ok --> G((done))

sss.cs內容說明

  • 年月資訊由所在目錄($PWD)中讀取
  • 使用date指令得到月份英文名稱(%b)，做為結果檔的主要檔名。
    • 早期使用name_of_month，有其便利性，但不如使用date指另單純。
    • date指令要注意locale的設定，%b結果不必然是英文的月份，要看LC_ALL的內容。
  • spe污染物項目，此處保留CO(cmo)，如程式沒有輸出CO，此部分結果只會有觀測值。
  • 批次方式執行grads，需使用 -b 指令，執行特定指令為 -c。-p 則為輸出結果圖紙面的方向(portrait)
yr=`echo $PWD|cut -d'/' -f4|cut -c-4`
mn=`echo $PWD|cut -d'/' -f6|cut -c4-5`
export MON=`LC_ALL=en_US.UTF-8 date -d "${yr}-${mn}-01" +%b`
for spe in so2 cmo ozn pmt nox p25 no2 voc;do
  cp ../spec.gs.BLK spec.gs
  for cmd in 3s/spe/$spe/ ;do
    sed -i $cmd spec.gs
  done
  rm ?.pdf
  grads -b -p -c 'run spec.gs'
  for j in {1..7};do
    convert $j.png $j.pdf
  done
  pdf_all "?" $MON$spe
done

gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=tmp.pdf result.pdf $i".pdf"
Search