Dr. Kuang's Utilities

本土化CALPUFF濃度預報系統之實現

Table of contents

背景

Leed大學CEMAC中心建置了Masaya 火山噴發SO2/SO4造成地面濃度的預報模式，範例如圖所示。由於氣象預報數據、大氣擴散模式等皆為官方作業系統與優選模式，因此具有高度的參考價值。

此處將其進行本土化，應用於空氣品質預報、緊急應變作業系統中。

該系統自動下載未來48小時MAM氣象數值預報的結果，轉成CALMET的輸入檔，進入CALPUFF模式進行地面空氣品質的模擬。

由於CALPUFF模式具有考慮3D風場、地形變化、化學反應等模擬能力，曾被美國環保署列為大氣擴散中的優選模式，用於大型固定污染源原生污染物以及二次氣膠長程傳輸的模式。

挑戰任務與成果

挑戰任務

系統工作流程由Run.sh批次檔所控制，重要節點及困難如下：

  1. 自動下載我國中央氣象局（CWB）數值預報（WRF）結果，並進行解讀。
    1. CWB會員登錄、檔案網址之定位與確認
    2. 氣象與空品模式的模擬範圍與解析度之決定
    3. 自動下載與儲存空間之預備
    4. 美國MAM預報內容與台灣WRF模式預報內容的差異與對照、
  2. 解讀程式（Create3DDAT.py）的安裝與修改
    1. Python環境之設定（內設為3.6版）
    2. grib-api模組之安裝（wasg不能import的問題）
    3. 座標系統的調整轉換（經緯度twd97 vs LCP）
  3. CALMET模式的偵錯
    1. Fortran程式編譯（內設使用pgf90，與ifort、gfortran版本的差異）
    2. 輸入資料間隔6小時程式的偵錯
  4. 排放量的設定
  5. CALPUFF模式輸出的修改
    1. Calpuff.inp中離散點的設定
    2. output.f 離散點之輸出
    3. PM2.5 的計算（結合銨鹽重量的計算）
  6. 等濃度圖與網站設定
    1. 模擬污染物項目、小時數、動畫、濃度等級的修改
    2. 指定外部IP（內設是本機）

執行方式

  1. 切換目錄到UNRESPForcastingSystem
  2. 編輯Run.sh打開要執行的程式
# Defaults that can be overwritten by editing HERE:
# Command line option m switches all to false
runTERREL=false (地形)
runCTGPROC=false（土地使用）
runMAKEGEO=false（綜合地理設定）
run3DDAT=false（下載WRF grib files）
runCALMET=false(執行calmet)
runCALPUFF=false（執行calpuff）
runmodel=true
  1. 將conda環境設成grippy: source ~/conda_ini grippy
  2. sh Run.sh -p -d yyyymmdd
    • -d 指定起始日期（之後48小時）
    • -p 繪製所有物質的圖形

排放位置等參數之改變

  1. 檔案模板：calpuff.inp
  2. 修改位置
    • 模式採用Lambert 投影座標系統，原點位在台灣中心點，北緯24.61，東經120.99.
    • (Taizhong PP)< 1 ! X = -49.873, 63.288, 150, 4.7, 11.00, 19.8, 363.0, .0, 5.440,
    • (Tongxiao PP)> 1 ! X = -32.064, 97.505, 150, 4.7, 11.00, 19.8, 363.0, .0, 5.440,

成果比較

unresp1.png
leed大學火山SO2預報結果畫面
unresp2.png
台中電廠PM2.5預報結果畫面

氣象數據下載與解讀

氣象數值預報的選擇

unresp3.png
CWB WRF 模式3KM預報結果風場

氣象與空品模式的模擬範圍與解析度之決定

unresp4.png unresp5.png
CWB_WRF 3維變數之高度分布 CWB_WRF數值預報地面項目 
conda activate ncl_stable
convert2nc M-A0064-000.grb2
  Python nc2db.py#將每一個格點的座標寫成csv檔案。

nc2db.py內容

import netCDF4
from pandas import *
import argparse
import numpy as np
parser = argparse.ArgumentParser(description="read NetCDF data and transform to DataFrame csvs. ",\
   epilog="Example of use: pythonb nc2df NC_FILENAME")
parser.add_argument("nc_filename", help="nc filebname",type=str)
args = parser.parse_args()
fname = args.nc_filename
nc = netCDF4.Dataset(fname)
cols=list(nc.variables)           #note 1
d={}
for c in ['gridlon_0','gridlat_0']:
    a=list(np.array(nc.variables[c]).flatten())           #note 3
    d.update({c: a})
df=DataFrame(d)
df.to_csv(fname+'.csv')

Create3DDAT.py的安裝與修改

Python環境之設定（內設為3.6版）

模擬範圍之選取

193,200c197,206
< latMinCP = 11.7  # Min lat of CALPUFF grid
< latMaxCP = 12.2  # Max lat of CALPUFF grid
< lonMinCP = 273.2  # Min lon of CALPUFF grid
< lonMaxCP = 274.1  # Max lon of CALPUFF grid
< inDir = '../NAM_data/raw/'+date  # Directory containing GRIB files
< nfiles = 17 #Number of GRIB files (files are 3 hourly, so 48 hours is 17 files including hours 0 and 48)
< outFile = '../NAM_data/processed/met_'+date+'.dat' #Output file path
< levsIncl = [1000,950,925,900,850,800,700,600,500,400,300,250,200,150,100,75,50,30,20,10,7,5,2] #pressure levels to include in output
---
> latMinCP = 21.4    # Min lat of CALPUFF grid
> latMaxCP = 25.7   # Max lat of CALPUFF grid
> lonMinCP = 119.4   # Min lon of CALPUFF grid
> lonMaxCP = 122.4   # Max lon of CALPUFF grid
> inDir =  '/home/kuang/MyPrograms/UNRESPForecastingSystem/CWB_data/raw/'+date  # Directory containing GRIB files
> nfiles = 15 #Number of GRIB files (files are 6 hourly, so 84 hours is 15 files including hours 0 and 48)
> outFile ='/home/kuang/MyPrograms/UNRESPForecastingSystem/CWB_data/processed/met_'+date+'.dat' #Output file path
> if os.path.isfile(outFile):
>     sys.exit('exist outFile: '+outFile)
> levsIncl = [1000,925,850,700,500,400,300,250,200,150,100] #pressure levels to include in output

座標系統的調整轉換（經緯度twd97 VS LCP）

< lats = gribapi.grib_get_array(gidPRMSL,'distinctLatitudes')
< lons = gribapi.grib_get_array(gidPRMSL,'distinctLongitudes')
< Ni = gribapi.grib_get(gidPRMSL,'Ni')
< Nj = gribapi.grib_get(gidPRMSL,'Nj')
---
> ishp=2
> CSV='/home/kuang/MyPrograms/UNRESPForecastingSystem/data/M-A0064nc2.csv'
> df=read_csv(CSV)
> #lats = gribapi.grib_get_array(gidPRMSL,'distinctygrid')#Latitudes')
> #lons = gribapi.grib_get_array(gidPRMSL,'distinctLongitudes')
> lats = np.array(list(df.gridlat_0))
> lons = np.array(list(df.gridlon_0))
> Ni = 1158
> Nj = 673

濃度場模擬範圍在氣象場中相對位置與切割

格點座標值由1維增加為(實質)2維

78,81c80,83
<     RXMIN = lons[iLonMinGRIB]  # W-most E longitude
<     RXMAX = lons[iLonMaxGRIB]  # E-most E longitude
<     RYMIN = lats[iLatMinGRIB]  # S-most N latitude
<     RYMAX = lats[iLatMaxGRIB]  # N-most N latitude
---
>     RXMIN = lons[ijLLMinGRIB]  # W-most E longitude
>     RXMAX = lons[ijLLMaxGRIB]  # E-most E longitude
>     RYMIN = lats[ijLLMinGRIB]  # S-most N latitude
>     RYMAX = lats[ijLLMaxGRIB]  # N-most N latitude
97c99
<         XLATDOT = lats[iLatMinGRIB+j]  # N latitude of grid point
---
>         XLATDOT = lats[ijLLMinGRIB+j*Ni]  # N latitude of grid point
100c102
<             XLONGDOT = lons[iLonMinGRIB+i]  # E longitude of grid point
---
>             XLONGDOT = lons[ijLLMinGRIB+i]  # E longitude of grid point

檔案名稱系統

累積降雨之轉接

...
import pygrib
...
filenames = [filePrefix+'18'+fileSuffix]
filePaths = [os.path.join(oldDir,filenames[0])]
...
# intialize rain by kuang
Ni = 1158
Nj = 673
fname=filePaths[0].replace('018','012')
rain_old=pygrib.open(filePaths[0])[62].values-pygrib.open(fname)[62].values

    RNgrd = np.zeros(shape=(Nj, Ni)) #by kuang
    for t in range(nfiles):
...
        #by kuang
        if t==0:
            RNgrd[:]=rain_old[:]
        else:
            RNgrd[:]+=pygrib.open(filePaths[t])[62].values[:]
...
        for j in range(NY):
            JX = j+1  # J-index of grid cell
            for i in range(NX):
                IX = i+1  # i-index of grid cell
...
                RAIN=int(RNgrd[iLatMinGRIB+j,iLonMinGRIB+i]) /10. #mm to cm
                fout.write(('{:4d}'+'{:02d}'*3+'{:3d}'*2+'{:7.1f}{:5.2f}{:2d}'+'{:8.1f}'*2+'\n').format(MYR,MMO,MDAY,MHR,IX,JX,PRES,RAIN,SC,RADSW,RADLW))

即時排放數據

數據來源及處理

項目 發電量即時資訊 CEMS數據 選擇考量
污染源完整性 所有發電機組都有 有的污染源沒有(可能因為規模太小未列管) 沒有數據就無從推估
污染項目 無(只有運轉率%) 煙氣流量、溫度、SOX、NOX及不透光率 後者太過複雜還有待QCQA
時間頻率 10min 15min 此處作業只需要逐時
資料穩定性 穩定 不甚穩定(CEMS允許一定時數之離線)  

前日運轉率之彙整與應用

$ head g20220412.csv
CP_NO,TEMP,VEL,CO_GPS,NMHC_GPS,NOX_GPS,PM25_GPS,PM_GPS,SOX_GPS,DateHr
F1700736P601,114,16.57785,0.0,0.0,21.48246822694089,1.1103541221509972,1.882682180377493,19.653922323336385,2022041200
F1700736P701,108,16.575588,0.0,0.0,16.499618542048232,1.1976660590277777,2.0307950954861114,12.677973511904762,2022041200

kuang@master /home/sespub/power
$ head p20220412.csv
CP_NO,HEI,DIA,TEMP,VEL,UTM_E,UTM_N
C1400170P301,202.0,6.4,127,20.0,321548.0,2777220.0
C1400170P401,205.0,6.4,118,20.0,321548.0,2777220.0
E5400878P001,60.0,6.2,150,9.1,180595.0,2503815.0

kuang@master /home/sespub/power
$ cat rd_today.py
#!/cluster/miniconda/envs/unresp/bin/python
from pandas import *
import subprocess, os
from pypinyin import pinyin, lazy_pinyin
cole=['CO_EMI', 'NMHC_EMI', 'NOX_EMI', 'PM25_EMI', 'PM_EMI', 'SOX_EMI']
colg=[i.replace('EMI','GPS') for i in cole]
hdtv=[ 'HEI', 'DIA', 'TEMP', 'VEL']

date=subprocess.check_output('date -d "-1 day" +%Y%m%d',shell=True).decode('utf8').strip('\n')
df=read_csv('names.csv')
unit2={}
col=['CP_NO']+hdtv+colg
for c in col:
  dd= {i:j for i,j in zip(list(df.UnitName),list(df[c]))}
  unit2.update({c:dd})

td=DataFrame({})
for t in range(24):
  tt='{:02d}'.format(t)
  fname='b'+date+tt+'.txt'
  if not os.path.exists(fname):continue
  with open(fname,'r') as f:
    lines=[i.strip('\n') for i in f]
  a=[i.split(',')[2] for i in lines[1:]]
  b=[]
  for i in a:
    ll=lazy_pinyin(i)
    s=''
    for ii in ll:
        s+=ii
    b.append(s)
  perc=[i.split(',')[5] for i in lines[1:] ]
  dft=DataFrame({'Name':b,'Perc':perc})
  dft=dft.loc[dft.Name.map(lambda x:x in set(df.UnitName))].reset_index(drop=True)
  dft.Perc=[float(i[:-1])/100. for i in dft.Perc]
  dft['DateHr']=[date+tt for i in range(len(dft))]
  dft['CP_NO']=[unit2['CP_NO'][i] for i in dft.Name]
  dft=dft.loc[dft.Perc>0]
  td=td.append(dft,ignore_index=True)
for c in colg:
  c0=np.array([unit2[c][i] for i in td.Name])
  td[c]=c0*td.Perc
for c in hdtv:
  c0=np.array([unit2[c][i] for i in td.Name])
  if c=='VEL':
    td[c]=c0*td.Perc
  else:
    td[c]=c0
td[col[:1]+col[3:]+['DateHr']].set_index('CP_NO').to_csv('g'+date+'.csv')

df=read_csv('names.csv')
df=df.loc[df.CP_NO.map(lambda x:x in set(td.CP_NO))].reset_index(drop=True)
df[col[:5]+['UTM_E','UTM_N']].set_index('CP_NO').to_csv('p'+date+'.csv')

CALPUFF逐時排放量檔案之準備

...
echo ${rundate} ${prevdate} ${enddate}> ptem_PWR.inp
/home/cpuff/2018/ptem/ptem_PWR<ptem_PWR.inp
if ! [ -e ptemarb_pwr.dat ]; then echo -n "fail ptemarb_PWR !";exit 0;fi
NPT2=$(head -n11 ptemarb_pwr.dat|tail -n1| awk '{print $1}')
...
  sed   -e ...
        -e "s/?NPT2?/$NPT2/g"  \
    ./CALPUFF_INP/calpuff_template.inp > ./CALPUFF_INP/calpuff.inp

CALPUFF系統的更新

CALMET模式的偵錯

  1. Fortran程式編譯（內設使用ifort、gfortran版本的差異）
    • dec舊式fortran可以接受未定寬度的輸出格式，如(i)、(f)等
    • 如read(cfver,'(f)') cfdataset
    • ifort/pgi皆可接受，但gfortran無法接受
      • 過去可以使用-fdec-format-defaults選項控制，但新版gfortran已經取消這個選項，不能接受未定寬度的輸出格式，詳gcc官網
      • 解決方法：因錯誤並不多，直接將程式碼中未定長度格式，改成自由格式(*)
  2. 輸入資料間隔(nlevag11)6小時程式的偵錯
    • 基本版本：CALMET_v6.5.0_L150223
    • ifort 新舊版本的差異
(unresp)
kuang@master /cluster/CALPUFF6/CALMET
$ diff CALMET_v6.5.0_L150223/CALMET.FOR new/calmet.for
21043c21045
<       save it1,it2,isec1,isec2
---
>       save it1,it2,isec1,isec2,nlevag11

kuang@master ~/MyPrograms/UNRESPForecastingSystem/CALPUFF_MODS/CALMET
$ diff calmet.for /cluster/CALPUFF6/CALMET/new/calmet.for
21045c21045
<       save it1,it2,isec1,isec2,nlevag1,msec
---
>       save it1,it2,isec1,isec2,nlevag11

CALPUFF模式輸出的修改

  1. Fortran程式編譯（內設使用ifort、gfortran版本的差異）
    • 與前述CALMET相同，也具有未定寬度的輸出格式的問題
      • dec舊式fortran可以接受未定寬度的輸出格式，如(i)、(f)等
      • in readcf.f and runsize.f  read(cfver,'(f)') cfdataset
    • 解決方法：直接將未定長度格式改成自由格式(*)
  2. calpuff.inp中離散點的設定
    • 研究團隊並沒有使用CALPOST做為後處理程式，而是利用CALPUFF模式逐時輸出離散點濃度的做法，將各污染物的逐時濃度逐一開啟檔案輸出。
    • CALPUFF.INP須輸入所有接受點的總數，以及各點的座標(LCP)，
    • 先由前述csv檔讀入lat/lon值，經由twd97模組將經緯度轉成twd97座標，扣除Xcent、Ycent之offset即成為LCP座標
    • CALPUFF為KM單位
    • 檔案另輸出成data/xy_taiwan3.dat
Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
Xcent, Ycent = twd97.fromwgs84(Latitude_Pole, Longitude_Pole)
  1. output.f 離散點之輸出
$ diff CALPUFF_v7.2.1_L150618/output.f new/output.f
125a126
>       character*150 infile
373a375,392
> c--------------------------------------------------------------
> ckuang
> c --- Sara 08/10/2004: Just to have an output of the concentration
> c --- directly readable.
> c --- Concentration in g/m^3 is written for each receptor
>             do ispec=1,nspec
>               write(infile,10000) ispec, istep
> c     &                     I2.2,I4.4,'.vtk')
> 10000         format('concrec',I2.2,I4.4,'.dat')
> c             print *, nspec,',',ispec,',',nn,',',infile,',',
> c     +       nrec,',',chirec(1,ispec)
>               open(50,file=infile)
>               do i=1,nrec
>                 write(50,*) chirec(i,ispec)
>               end do
>               close(50)
>             end do
  1. PM2.5 的計算（結合銨鹽重量的計算）
    • CALPUFF之CSPEC共有8種，利用其中的成份計算總PM2.5，包括結合銨鹽
    • 以DataFrame架構整理檔案名稱(spec, hour)，再讀入個別檔案的內容成為3維矩陣C。
    • 輸出檔案為00-hhhh.dat。
    • 程式下載點conc2pm25.py
kuang@master ~/MyPrograms/UNRESPForecastingSystem/Python
$ cat conc2pm25.py
import numpy as np
import os,sys,subprocess
from pandas import *
CSPEC='SO2 SO4 NOX HNO3 NO3 PMS1 PMS2 PMS3'.split()
fnames=list(subprocess.check_output('ls concrec*dat',shell=True).split(b'\n'))
fnames=[i.decode('utf8') for i in fnames if len(i)>0 ]
if len(fnames)==0:sys.exit('concrec not found')
wc=int(subprocess.check_output('cat '+fnames[0]+'|wc -l',shell=True).split(b'\n')[0])
jt=[int(fname.split('/')[-1].replace('.dat','')[-4:]) for fname in fnames if len(fname)>0]
js=[int(fname.split('/')[-1].replace('.dat','')[-6:-4]) for fname in fnames if len(fname)>0]
df=DataFrame({'hr':jt,'spec':js,'fname':fnames})
df=df.loc[df.spec>0].reset_index(drop=True)
C=np.zeros(shape=(max(js)+1,max(jt)+1,wc))
for i in range(len(df)):
  with open(df.loc[i,'fname'],'r') as f:
    tmp=[float(l.strip('\n')) for l in f]
  C[df.loc[i,'spec'],df.loc[i,'hr'],:]=tmp[:]
so4 =C[1,:,:]
hno3=C[3,:,:]
no3 =C[4,:,:]
p25 =C[5,:,:]
nh4=so4*(36./96.)+no3*(18./62.)+hno3*(18./63.)
total=so4+no3+hno3+nh4+p25
fnRoot=fnames[0].replace('.dat','')[:-6]+'00'
for it in range(1,max(jt)+1):
  fname=fnRoot+'{:04d}'.format(it)+'.dat'
  with open(fname,'w') as f:
    for ic in range(wc):
      f.write(str(total[it,ic])+'\n')

等濃度圖與網站設定(格柵底圖方案)

32d31
< import utm
35a35,36
> import twd97
123,128c128,131
>     Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
>     Xcent, Ycent = twd97.fromwgs84(Latitude_Pole, Longitude_Pole)
>     ll = np.array([twd97.towgs84(i*1000+Xcent,j*1000+Ycent) for i, j in zip(x, y)])
>     lat, lon = (ll[:, i] for i in [0, 1])
166,169c169,172
>     Latitude_Pole, Longitude_Pole = 23.61000, 120.9900
>     Xcent, Ycent = twd97.fromwgs84(Latitude_Pole, Longitude_Pole)
>     ll = np.array([twd97.towgs84(i*1000+Xcent,j*1000+Ycent) for i, j in zip(x2, y2)])
>     lat, lon = (ll[:, i] for i in [0, 1])

273,274c276,277
<         self.binLims = [10, 350, 600, 2600, 9000, 14000]  # SO2 bin limits
<         self.binLimsSO4 = [1E-8, 12, 35, 55, 150, 250]  # SO4 bin limits from:
---
>         self.binLims = [1E-3,1E-2,1,5, 10, 350]  # SO2 bin limits
>         self.binLimsSO4 = [1E-3, 1E-2,1E-1, 1, 5, 15]  # SO4 bin limits from:

399,400c403,404
>         latTicks = np.arange(round(latMin, 1), round(latMax, 1) + 0.1, 0.5)
>         lonTicks = np.arange(round(lonMin, 1), round(lonMax, 1) + 0.1, 0.5)

cd /home/cpuff/UNRESPForecastingSystem/VIZ_SITE_CODE/public_html
/cluster/anaconda3/bin/python3.7/python -m http.server --bind 200.200.12.191 8030 >&/dev/null &

VERDI向量底圖方案

向量底圖

...
    for i in $(ls ${s}${static}??.PNG);do
      /usr/bin/convert ${i} -bordercolor white -trim tmp.png
      /usr/bin/convert tmp.png -bordercolor white -border 5%x5% ${i}
    done
    /usr/bin/convert -dispose 2 -coalesce +repage -background none ${s}${static}??.PNG ${s}.gif
    cp ${s}.gif /var/www/html/LC-GIF-Player/example_gifs
...

網站與播放器

Download Run.sh

Reference

This site is open source. Improve this page.