*make coords for xr.dataarray
xarrayのdataarrayがあれば，to_netcdf関数でnetcdfデータができるのだけど，適切に座標系を定義してあげないと，後でxarrayなどで読み込んでうまく使えない．また，lat, lon gridであるならGrADsのsdfopenでも読めるほうがお得だ．以下のようにするとxarrayで読めば座標情報が使えるし，GrADsのsdfopenでも読めるようだ．GrADsで読むためにはlon, latのunitsを指定することが必須．これがないとsdfopenで(ctl fileなしでは)読めない．時間の設定はちょっと面倒でもっと良い方法がありそう．import numpy as npimport pandas as pdimport xarray as xrlon1d=np.arange(0.5, 360., 1.0).astype(np.float32)lat1d=np.arange(-89.5, 90., 1.0).astype(np.float32)xsz=len(lon1d)ysz=len(lat1d)tsz=10data=np.random.randn(tsz,ysz,xsz)# 時間の設定time_list=[]for t in range(tsz): strtime='2000-%2.2d-15'%(t+1) time_list.extend([strtime])datetime_array=pd.to_datetime(time_list)xr_data=xr.DataArray(np.float32(data), name='slp', \ # ディスク容量を節約するためにfloat32で coords={'time':datetime_array,\ 'lat':('lat',lat1d,{'units':'degrees_north'}), 'lon':('lon',lon1d,{'units':'degrees_east'})}, \ dims=['time','lat','lon'])
#---------- 次の例では，時間を含む任意次元数かつ任意独立変数名のデータ# をもとにして，時間以外は共通にして時間を新しく用意したdatetime_array# で設定しているin_xc=data_xa.coordsndimsz=len(in_xc.dims)out_xc={'time':datetime_array}for ndim in np.arange(1,ndimsz): dim_name=in_xc.dims[ndim] coord_values=in_xc[dim_name].values out_xc[dim_name]=(dim_name,coord_values,{'units':in_xc[dim_name].units,'axis':in_xc[dim_name].axis})data_out_xa=xr.DataArray(np.float32(data_out), name=var, \ coords=out_xc, dims=in_xc.dims)print(data_out_xa)data_out_xa.to_netcdf('x.nc')
#---------------次の例では入力よりも切り詰めた時間を出力データにセットしているdata_xs=xr.open_dataset(fn_fl,decode_cf=False)time_data=data_xs.coords['time']tsz=len(time_data)# cut the first & last 2 years for coorinateout_coords=data_xs[var].isel(time=slice(2,tsz-2)).coords data_xa=data_xs[var]data_na=np.array(data_xa)data_na[data_na<-1.e+7]=np.nandata_rav5_na=runave_wei.runavnd(data_na,5,5,axis=0)# cut the first & last 2 years for datadata_rav5_out_na=data_rav5_na[2:-2,:,:,:]data_rav5_xa=xr.DataArray(np.float32(data_rav5_out_na), name=var, \ coords=out_coords, \ dims=data_xa.dims)
#-------------- 次の例では４次元（時間，鉛直，南北，東西）のデータから作った相関マップデータに，元の座標から時間を除いた３次元の座標情報（鉛直，南北，東西）を付けている．coords=mvdt_xrds.coordsdel(coords['time'])xr_data=xr.DataArray(np.float32(crr_map), name='crr', \ coords=coords, dims=['plev','lat','lon'])
#------------- 次の例では一か月一つのデータから12カ月の気候値を作り，timeを足すためにcoordsを直接指定している．time_list=[]for mth in range(1,12+1): strtime='2000-%2.2d-15'%mth time_list.extend([strtime])datetime_array=pd.to_datetime(time_list)xr_data=xr.DataArray(clm, name=var_name, \ coords={'time':datetime_array, \ 'depth':('depth',np.array(temp.coords['DEPTH']),{'units':'meters'}), \ 'lat':('lat',np.array(temp.coords['LAT']),{'units':'degrees_north'}), \ 'lon':('lon',np.array(temp.coords['LON']),{'units':'degrees_east'})}, \ dims=['time','depth','lat','lon'])
#-------- 次の例では時間を新しくpandasのdate_rangeで定義して元からある時間と取り替えている。なおfrequencyを'M'にすると月の終わりの日が設定され，これは気候科学では一般的ではない。より一般的な月初めを設定するために，'MS'を指定している。
#---------- 次の例では，時間を含む任意次元数かつ任意独立変数名のデータ# をもとにして，時間以外は共通にして時間を新しく用意したdatetime_array# で設定しているin_xc=data_xa.coordsndimsz=len(in_xc.dims)out_xc={'time':datetime_array}for ndim in np.arange(1,ndimsz): dim_name=in_xc.dims[ndim] coord_values=in_xc[dim_name].values out_xc[dim_name]=(dim_name,coord_values,{'units':in_xc[dim_name].units,'axis':in_xc[dim_name].axis})data_out_xa=xr.DataArray(np.float32(data_out), name=var, \ coords=out_xc, dims=in_xc.dims)print(data_out_xa)data_out_xa.to_netcdf('x.nc')
#---------------次の例では入力よりも切り詰めた時間を出力データにセットしているdata_xs=xr.open_dataset(fn_fl,decode_cf=False)time_data=data_xs.coords['time']tsz=len(time_data)# cut the first & last 2 years for coorinateout_coords=data_xs[var].isel(time=slice(2,tsz-2)).coords data_xa=data_xs[var]data_na=np.array(data_xa)data_na[data_na<-1.e+7]=np.nandata_rav5_na=runave_wei.runavnd(data_na,5,5,axis=0)# cut the first & last 2 years for datadata_rav5_out_na=data_rav5_na[2:-2,:,:,:]data_rav5_xa=xr.DataArray(np.float32(data_rav5_out_na), name=var, \ coords=out_coords, \ dims=data_xa.dims)
#-------------- 次の例では４次元（時間，鉛直，南北，東西）のデータから作った相関マップデータに，元の座標から時間を除いた３次元の座標情報（鉛直，南北，東西）を付けている．coords=mvdt_xrds.coordsdel(coords['time'])xr_data=xr.DataArray(np.float32(crr_map), name='crr', \ coords=coords, dims=['plev','lat','lon'])
#------------- 次の例では一か月一つのデータから12カ月の気候値を作り，timeを足すためにcoordsを直接指定している．time_list=[]for mth in range(1,12+1): strtime='2000-%2.2d-15'%mth time_list.extend([strtime])datetime_array=pd.to_datetime(time_list)xr_data=xr.DataArray(clm, name=var_name, \ coords={'time':datetime_array, \ 'depth':('depth',np.array(temp.coords['DEPTH']),{'units':'meters'}), \ 'lat':('lat',np.array(temp.coords['LAT']),{'units':'degrees_north'}), \ 'lon':('lon',np.array(temp.coords['LON']),{'units':'degrees_east'})}, \ dims=['time','depth','lat','lon'])
#-------- 次の例では時間を新しくpandasのdate_rangeで定義して元からある時間と取り替えている。なおfrequencyを'M'にすると月の終わりの日が設定され，これは気候科学では一般的ではない。より一般的な月初めを設定するために，'MS'を指定している。