北極振動 パターン

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