北極振動指数の時系列 2/2

*
*
前回はEOFと同時に得られる時系列をプロットした。今回はEOFの期間外の北極振動指数を求めてプロットする。例えば2008年から2009年までの値を求めてプロットすると次のとおり。赤が自分で求めた値で、黒がNOAAの発表している値。ほとんどあっている。

 

以下がプログラム。

(1) python interface for grads を使い、2008年から2009年の1000hpaのgeopotential heightを読み込む。データソースはNCEP/NCAR reanalysis。

(2) EOFを求めるときにセーブしておいた、12ヶ月分の1979-2000年気候値(mhgt)を読み込んで、差し引く。mhgtはnumpyのtile関数を用いて時間方向に2回(2年分)の繰り返しにする。

(3)　緯度のファクターを掛ける。緯度90度の取扱いはEOFを求める時とおなじようにしておく。

(4)EOFを求めるときにセーブしておいたEOF1のパターンにプロジェクションをおこなう。また、北極振動を規格化するための定数pfactor1で、EOFを規格化するのと、時系列を規格化するので2回割る。

 

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

import scikits.timeseries as ts        #scikits timeseries module
import scikits.timeseries.lib.plotlib as tpl

#------------  data import using python interface for grads　(1)

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("/home/tmiyama/DATA/Reanalysis/NCEP1/hgt.mon.mean.nc") # data open
ga("set lat 20 90")
ga("set time  JAN2008 DEC2009") # use data from 2008 to 2009
ga("set lev 1000") # 1000 hpa
hgt = ga.exp("hgt")
del ga

#------------- Get AO index  ----------------------------

#remove seasonal climatogical (2)
mhgt=np.load('mhgt.npy')   
hgt=hgt-np.tile(mhgt,(2,1,1))

# multiply latitude factor　(3)　
lat=hgt.grid.lat
lat[-1]=lat[-1]+1.E4
factor = np.sqrt(np.cos(np.pi*(lat/180.)))
datain = hgt[:,:,:] * factor[np.newaxis,:,np.newaxis]


# projection to EOF1 (4)
nt,ny,nx=datain.shape
eof1=np.load('eof1.npy')
pc1=np.dot(np.reshape(datain.data,(nt,ny*nx)),np.reshape(eof1,(ny*nx,1)))
pfactor1=np.load('pfactor1.npy')
pc1=pc1/pfactor1/pfactor1

dbeginmonth=ts.Date('M', '2008-1')
AO_2008_2009_series=ts.time_series(pc1,start_date=dbeginmonth,freq='M') # make timeseries object

# plot
from readAO import readAO
AOindex_series=readAO() # import NOAA data

fig = tpl.tsfigure()
fsp = fig.add_tsplot(111)
l1=fsp.tsplot(AOindex_series, 'k-')
l2=fsp.tsplot(AO_2008_2009_series, 'r-')
istart=ts.Date('M','2008-1').value
iend=ts.Date('M','2009-12').value
fsp.set_xlim(istart,iend)
plt.xlabel("YEAR")
plt.ylabel("AO index")
plt.legend((l1,l2),("NOAA","What I got"))
plt.savefig("check_AOindex-2008-2009.png")
plt.show()