*
*
NOAA Climate Prediction Centerの記述
http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/history/method.shtml を参考にEOF解析を用いて北極振動のパターンをもとめる。
用いる”python interface for grads”には
EOF関数があるのだが、計算がうまくいかなかったので(緯度が90度の点があるとうまくいかないようだ）、自分でつぎのようなモジュールを作っておく。
”python interface for grads”から
from grads.ganum import GaNum
ga = GaNum(Bin='grads')
ga.eof??
で関数の中身がわかるのでそれを参考に作る。
eof.py
def eof (data):
"""
calculates Empirical
Orthogonal Functions (EOFS) using Singular Value Decomposition (SVD).
V, d, c = eof(data)
where
V --- eigenvectors
d --- NumPy array with eigenvalues
c --- NumPy array with principal components
"""
import numpy as np
# At least 2 time steps
# ---------------------
if data.shape[0] < 2:
raise GrADSError, \
'need at least 2 time steps for EOFS but got nt=%d'% data.shap
e[0]
nt, ny, nx = (data.shape[0], data.shape[1], data.shape[2])
# Export N-dimensional array
# --------------------------
u = data
# Reshape as 3D
# -------------
u = u.reshape(nt,ny,nx)
# ------------------------------------------
pc, d, u = np.linalg.svd(u.reshape(nt,ny*nx),
full_matrices=0)
# Eigenvalues/coefficients
# ------------------------
d = d * d
pc = pc.transpose()
# Reshape eigenvectors as original array
# --------------------------------------
u = u.reshape(nt,ny,nx)
return (u, d, pc)
モジュールが準備できたら、EOF解析を実行する。
使用するモジュールをインポートする。最後に、上のeofモジュールをインポートしている。
from grads.ganum import GaNum # for python interface for grads
from mpl_toolkits.basemap import Basemap # for Basemap toolkit
import matplotlib.pyplot as plt
import numpy as np
from eof import eof # import EOF
python interface for gradsを用いて、1971年か2000までのデータを読み込む。データはOPeNDAPを用いてNCEP/NCARのデータを
NOAA ESRLのサイトから読み込んでいる。あるいはローカルにnetcdfを落としておいてそこから読んでも良い。1000hpaのGeopotential Heightを用いる。
Bin GuanさんのGradsスクリプトを用いて季節変動(mhgt)を引いて、アノマリ（ahgt)だけを取り出す。
#------------ data import using python interface for grads
ga = GaNum(Bin='grads')
fh=ga.open("http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.mean.nc") # data open
#fh=ga.open("hgt.mon.mean.nc") # data open
ga("set lat 20 90")
ga("set time JAN1979 DEC2000") # use data from 1979 to 2000
ga("set lev 1000") # 1000 hpa
ga("deseason hgt ahgt mhgt JAN1979 DEC2000") #use "deseason" at http://www.atmos.umd.edu/~bguan/grads/GrADS_Scripts.html agt is anomaly
ahgt = ga.exp("ahgt")
mhgt = ga.exp("mhgt")
del ga
EOF解析をおこなう。cos（緯度)^0.5のファクターをかけるのであるが、緯度90度の部分があると、svdがもとまらなかった。どうするのが正しいかよく分からないが、ちょっとだけ緯度をずらしてみた。北極振動指数は、1979年から2000年までの時系列の分散が1になるように設定されているので、そのように規格化する。また通常の定義にあうように符号も逆転させておく。
#------------- EOF ----------------------------
lat=ahgt.grid.lat # latitude
lon=ahgt.grid.lon
lat[-1]=lat[-1]-1.E-4 # lat=90 does not work
factor = np.sqrt(np.cos(np.pi*lat/180.))
datain = ahgt[:,:,:] * factor[np.newaxis,:,np.newaxis]
v,d,c=eof(datain)
pc1=c[0,:]*np.sqrt(d[0])
pfactor1=pc1.std(axis=0)
pc1=-pc1/pfactor1
eof1=-v[0,:]*pfactor1
あとはプロット。cos（緯度)^0.5ももとに戻しておく。後で使うかもしれないのでデータも保存しておく。
#------------- plot loading pattern
fig=plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
map = Basemap(projection='nplaea',boundinglat=20,lon_0=0.,resolution='i')
map.drawcoastlines()
lons,lats=np.meshgrid(lon[:],lat[:]+1.E-4)
x,y=map(lons,lats)
cz=map.contourf(x,y,eof1/factor[:,np.newaxis],[-40,-35,-30,-25,-20,-15,-10,-5,5,10,15,20],extend='both',cmap=plt.cm.RdBu_r)
# get axes position, add colorbar axes to right of this.
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.05, b, 0.03, h]) # setup colorbar axes.
plt.colorbar(cax=cax) # draw colorbar
plt.savefig("AO_loading.png")
plt.show()
# save data for later use
np.save('eof1.npy', eof1.data)
np.save('pc1.npy', pc1.data)
np.save('mhgt.npy',mhgt[0:12,:,:].data)
np.save('pfactor1.npy',pfactor1)
np.save('lon.npy',lon)
np.save('lat.npy',lat)
できたプロット。
NOAAのサイトにあるパターンと似たパターンを得ることができた。
このパターン（EOF第一モード）が説明する割合は19%。
d[0]/d.sum(axis=0)
0.19079784
2010/2/7 追記
緯度90度でうまく計算できないのは、計算精度により、cos90度 が負になり、その平方根をとっているからのようである。
2010/2/7 訂正
pfactor1はmasked arrayではないので、
np.save('pfactor1.npy',pfactor1.data)からnp.save('pfactor1.npy',pfactor1)に修正。
2010/2/8 追記
緯度(lat)、経度(lon)も保存するようにした。
np.save('lon.npy',lon)
np.save('lat.npy',lat)
