EARR 數據批次轉檔
背景
- 由於EARR1數據檔案是按小時切割，因此檔案格式繁多，如果一一轉檔將會造成檔案管理上的困難。
- 策略上將全年檔案讀過一遍，寫出平均結果，這樣應該較為單純。
- 檔案期間跨越2010～2019，因此如果可以平行運作，將可大幅提高處理效率。
- 目標雖然是鄉鎮區平均，但因為EARR解析度12公里，如果沒有內插，很多鄉鎮區將會有一樣的值、或找不到值。新網格系統參照tempTW.nc,公版模式範圍、解析度為1Km。
- 內插採用
scipy.interpolate.griddata模組（參考Near Real Time的檔案處理），使用cubic spline，以使空間的變化趨勢較為平緩。 - GFS再分析數據的年期範圍較小(2015~)，解析度也較差(0.25度約20Km)。然而因具有行星邊界層數據，是EARR沒有的，因此，也須加以處理。因程式邏輯完全一樣，也就在此處一併介紹。
EARR2csv程式說明
執行方式
for yr in 20{10..19};do sub python EARR2csv.py $yr;done
IO’s
- 引數：
yr4碼年代（2010～2019） - EARR檔案：按照內容分別建次目錄(
./slev,./precip) - 鄉鎮區檔案目錄：
root=/nas2/cmaqruns/2019TZPP/output/Annual/aTZPP/LGHAP.PM25.D001 - 1公里解析度模版：
$root/tempTW.nc - 網格-鄉鎮區對照表：
$root/gridLL.csv，此檔案由mk_gridLL所產生。 - 結果檔案：
'EARR'+yr+'.csv'
網格內插之準備
#1X1 TWN domain
root='/nas2/cmaqruns/2019TZPP/output/Annual/aTZPP/LGHAP.PM25.D001'
fname=root+'/tempTW.nc'
nc = netCDF4.Dataset(fname,'r')
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[V[3][0]].shape[i] for i in range(4))
pnyc = Proj(proj='lcc', datum='NAD83', lat_1=nc.P_ALP, lat_2=nc.P_BET,lat_0=nc.YCENT, lon_0=nc.XCENT, x_0=0, y_0=0.0)
#new coordinates
x1d=[nc.XORIG+nc.XCELL*i for i in range(ncol)]
y1d=[nc.YORIG+nc.YCELL*i for i in range(nrow)]
x1,y1=np.meshgrid(x1d,y1d)
ncXCELL=nc.XCELL
ncYCELL=nc.YCELL
#old coordinates
fname='slev/EARR.slev.20100606.t12z.nc'
nc = netCDF4.Dataset(fname,'r')
V0=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
for v in V0[1]:
exec(v+'=nc["'+v+'"][:]')
V0=[i for i in V0[1] if i[0] != 'l']+['precip_acc6h']
nc.close()
x,y = pnyc(lon, lat, inverse=False)
#neglect the old points outside new boundaries
maxx,maxy=x1[-1,-1],y1[-1,-1]
minx,miny=x1[0,0],y1[0,0]
boo=(x<=maxx+ncXCELL*10) & (x>=minx-ncXCELL*10) & (y<=maxy+ncYCELL*10) & (y>=miny-ncYCELL*10)
idx = np.where(boo)
mp=len(idx[0])
xyc= [(x[idx[0][i],idx[1][i]],y[idx[0][i],idx[1][i]]) for i in range(mp)]
鄉鎮區及網格之對照
#township mapping
df=read_csv(root+'/gridLL.csv')
df.TOWNCODE=['{:08d}'.format(i) for i in df.TOWNCODE]
df.COUNTYCODE=['{:05d}'.format(i) for i in df.COUNTYCODE]
tn={i:j for i,j in zip(df.TOWNCODE, df.TOWNNAME)}
cn={i:j for i,j in zip(df.COUNTYCODE, df.COUNTYNAME)}
時間及檔案管理
#temporal range
yr=sys.argv[1]
bdate=datetime(int(yr),1,1)
edate=datetime(int(yr),12,31)
ndate=bdate
df0=DataFrame({})
# file directories
dirs=['slev','precip']
tail={'slev':'.nc','precip':'.f06h.nc'}
ndiv={'slev':4,'precip':1}
逐日解讀與內插
4個迴圈的順序
- 使用while loop進行日的迴圈
- 依序讀取地面氣象要素、以及降雨量檔案
- 小時的迴圈
- 變數類別的迴圈
因相同邏輯套用在不同變數，除數(分母ndiv[d])需有所不同：
- 降雨量：全日之加總，除數為1。
- 其他項目：4次之平均值，除數為4。
#daily loop
while ndate <= edate:
nowd=ndate.strftime("%Y%m%d")
zz={v:np.zeros(shape=(nrow,ncol)) for v in V0}
for d in dirs:
for h in range(0,24,6):
fname=d+'/EARR.'+d+'.'+nowd+'.t{:02d}z'.format(h)+tail[d]
nc = netCDF4.Dataset(fname,'r+')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
for v in V0:
if v not in V[1]:continue
c=np.array([nc[v][idx[0][i], idx[1][i]] for i in range(mp)])
zz[v][:,:]+=griddata(xyc, c[:], (x1, y1), method='cubic')/ndiv[d]
每日結果寫成DataFrame
- 使用
pandas.pivot_table來進行鄉鎮區內網格之平均。 - 因大多數網格的累積雨量是0，cubic spline結果會出現負值，因此在平均之前需將其歸零。
- 山地鄉的範圍大、高程差異也大，氣壓、水汽濃度、氣溫等等可能會差異很大，其平均值的解釋需小心。
- 使用
pandas.DataFrame.append來附加每日轉換結果
for v in V0:
df[v]=zz[v][:,:].flatten()
precip=df['precip_acc6h']
precip=np.where(precip>=0,precip,0)
df['precip_acc6h']=precip
df_tm=pivot_table(df,index='TOWNCODE',values=V0,aggfunc=np.mean).reset_index()
df_tm['COUNTYCODE']=[i[:5] for i in df_tm.TOWNCODE]
df_tm['COUNTYNAME']=[cn[i] for i in df_tm.COUNTYCODE]
df_tm['TOWNNAME']=[tn[i] for i in df_tm.TOWNCODE]
df_tm['YMD']=nowd
df0=df0.append(df_tm,ignore_index=True,sort=False)
print(nowd)
ndate=ndate+timedelta(days=1)
結果輸出
df0.set_index('YMD').to_csv('EARR'+yr+'.csv')
程式下載
Download: EARR2csv.py