Dr. Kuang's Utilities

跨批次ncf檔之串連與均勻化

 
  • Oct 05, 2022
  • views
>

背景

  • 將CCTM_ACONC或者是combine結果檔，按照時間軸以ncrcat串連後，因在換日時更換了起始濃度與執行批次，致使等濃度結果(gif檔)產生跳動、停頓等等在時間軸上不連續、不合理情形。
  • 這一題類似intp_json，差異在於程式作用在nc檔案還是json檔案、以及漸變的時間周期。nc檔是為做m3nc2gif.py（see [[4.m3nc2gif]] ) 此處漸變週期24小時），json檔則是為earth顯示系統（週期6小時）。
  • 跨批次nc檔案的辨識，可以利用nc檔全域屬性中的CDATE及CTIME(current date/time)標籤：
    • 同一執行批次檔案CDATE/CTIME時間會彼此相近，後產生的檔案略大於前者，因此可以直接串連並無問題。
    • 不同批次結果之間的CDATE/CTIME時間會有較大的差距、甚至負值(重新執行舊的案例)。因此以datetime進行時間差異的計算，即可予以辨識、並進行必要之修正。
  • 在整體預報系統中([[2022-08-20-CMAQ_fcst]])，本程式屬於後處理部分

join_NCs.py程式設計

IO

  • 引數：具時間（日期）順序之nc檔、以及結果檔名（如已存在，將會被覆蓋）
    • python join_NCs.py pmt_1.nc pmt_2.nc pmt_3.nc pmt_4.nc pmt_5.nc pmt_6.nc pmt_7.nc pmt.nc
  • 各濃度檔

是否同一批次的辨識

  • 時間差在1小時之內，為同一執行批次的結果。0<=delt<3600.
  • datetime的時間差可以輸出年、日、時，但是都是整數，如果要精確計算，只能輸出秒數(total_seconds())再予以轉換。

時間軸之漸變

  • 周期：nt=24
  • 加權方式：
    • 線性加權平均
    • var0[:,t,:,:,:]=var[:,t0+t,:,:,:]*(nt-t)/nt + var[:,t1,:,:,:]*t/nt

檢討

  • 由於每天更新初始濃度，對CWBWRF_45k等大範圍網格而言算是很頻繁，因為GFS與ECMWF氣象場的差異還不小，致使煙團的移動也顯著的停頓、甚至倒退的情形。
    • GFS經過WRF的四階同化之後，在海面上的移動速度較ECMWF快一些。
    • REAS雖是較舊的月均值排放量，但比起CAMS排放模式略低一些，初始化後濃度消散很快。
  • 解決方案
    1. 開啟風蝕揚塵機制（CTM_WB_DUST)，雖然公版模式考量到CMAQ揚塵機制可能高估而並未建議開啟，在REAS整體低估的情況下，開啟風蝕揚塵機制有顯著提升濃度的效果。
    2. 以CAMS初始濃度修正REAS排放量。
    3. 重新以事件之角度重新進行批次之切割與組合、並重新執行模式模擬。

程式碼

#kuang@DEVP /nas2/cmaqruns/2022fcst/grid45/cctm.fcst/daily
#$ cat join_NCs.py
#!/cluster/miniconda/envs/unresp/bin/python
import os, sys
import netCDF4
import numpy as np
from dtconvertor import dt2jul, jul2dt

nf=len(sys.argv)-1
fnames=[sys.argv[i+1] for i in range(nf)]
cdates, tflags=[],[]
for n in range(nf-1):
  nc = netCDF4.Dataset(fnames[n],'r')
  t0,t1=n*24,(n+1)*24
  if n==0: #first time declarations
    V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
    nt,nlay,nrow,ncol=nc.variables[V[3][0]].shape
    nv=len(V[3])
    var=np.zeros(shape=(nv,(nf-1)*24,nlay,nrow,ncol))
  for iv in range(nv): #store the variables and current time/flags of files
    var[iv,t0:t1,:,:,:]=nc.variables[V[3][iv]][:,:,:,:]
  cdates.append(jul2dt([nc.CDATE,nc.CTIME]))
  tflags.append([jul2dt(nc['TFLAG'][t,0,:]) for t in range(24)])
nt=24 #time interplation period, in hours
for n in range(nf-2):
  delt=(cdates[n+1]-cdates[n]).total_seconds()
  if 0<=delt<3600.:continue # if conc are from same batch, skip interpolation
  t0,t1=(n+1)*24-nt,(n+1)*24
  var0=np.zeros(shape=(nv,nt,nlay,nrow,ncol))
  for t in range(nt): # transition in nt hours
    var0[:,t,:,:,:]=var[:,t0+t,:,:,:]*(nt-t)/nt +  var[:,t1,:,:,:]*t/nt
  var[:,t0:t1,:,:,:]=var0[:,:,:,:,:]

os.system('cp '+fnames[-2]+' '+fnames[-1])
nc = netCDF4.Dataset(fnames[-1],'r+')
#expand the nc file
for n in range(nf-1):
  t0,t1=n*24,(n+1)*24
  for t in range(t0,t1):
    nc.variables['TFLAG'][t,0,:]=dt2jul(tflags[n][t-t0])
var0=np.zeros(shape=nc.variables['TFLAG'].shape)
var0[:,:,:]=np.array(nc.variables['TFLAG'][:,0,:])[:,None,:]
nc.variables['TFLAG'][:,:,:]=var0[:,:,:]

for iv in range(nv):
  nc.variables[V[3][iv]][:,:,:,:]=var[iv,:,:,:,:]
nc.close()
Search