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()