ensemble output 4 grads

pythonで作成したデータを，データ可視化アプリケーションであるgradsでも読めると作業を効率的に行うことができる．ここで問題になるのが，ensembleである．ensembleの無いデータは多くの場合，netcdfで出力すれば，sdfopenで読むことができる．しかしensembleはgradsでその開発終盤に取り込まれた機能であるため，その対応が弱い．

gradsでensembleを含むnetcdf dataを読み込むのは，全アンサンブルを一つのnetcdf file にまとめて，xdfopenでopenすることで可能である．

コントロールファイルの例は以下の通りで，これをgradsセッション内で，xdfopen xxx.ctl でオープンする．

DSET ^./nc/m001-100.nc

DTYPE netcdf

TITLE ensemble test data

undef -1.e+10

XDEF lon

YDEF lat

TDEF time

EDEF ensemble

VARS 1

data=>data 0 data

ENDVARS

これでオープンするためのnetcdf file (また参考のためにbinary fileも)の作成用python scriptは以下の通りである．

import xarray as xr

import numpy as np

import pandas as pd

lons=np.arange(2.5,360,5)

lats=np.arange(-87.5,90,5)

ysz=len(lats)

xsz=len(lons)

esz=100

tsz=20

xydata=np.cos(lats/90*np.pi*2).reshape(ysz,1) @ np.cos(lons/180*np.pi*2).reshape(1,xsz)

xytedata=np.zeros([esz,tsz,ysz,xsz])

#xytedata=np.zeros([esz,tsz,ysz,xsz],dtype='float32')

for e in range(esz):

for t in range(tsz):

xytedata[e,t,:,:]=np.cos(t/10*np.pi*2) * xydata.reshape(1,ysz,xsz) + e

time_list=[]

for t in range(tsz):

strtime='%4.4i-01-15'%(t+2000)

time_list.extend([strtime])

datetime_array=pd.to_datetime(time_list)

es=np.arange(esz)+1

xr_data=xr.DataArray(xytedata, name='data',

coords={'ensemble':('ensemble',es,{'units':'number'}),

'time':datetime_array,

'lat':('lat',lats,{'units':'degrees_north'}),

'lon':('lon',lons,{'units':'degrees_east'})},

dims=['ensemble','time','lat','lon'])

xr_data.to_netcdf('./nc/m001-100.nc')

# ---------- binary outputs for grads -------------

for e in range(esz):

f=open('./bin/m%3.3i.bin'%(e+1),'wb')

f.write(xytedata[e,:,:,:])

f.close()

f=open('./bin/m001-100.bin','wb')

f.write(xytedata)

f.close()