串本・浦神の潮位差と黒潮流路緯度 Scipyでregression

regressionはscipyでも行える。stats.linregress関数（実際にはmasked array用のstats.mstat.linregress)を使う。
参照　Cookbook/Linear Regression

#read module
from scipy import stats

from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()   

#calculate regression
(a_s,b_s,r,prob,stderr,c)=stats.mstats.linregress(kushimoto_uragami,kuroshio_lat)
print "slope=",a_s
print "intersect=",b_s
print "correlation=",r
print "p value=",prob
print "std. of residuals",stderr
print "(Var. Latitude)/(Var. Sea level)=",c

最後のcはマニュアルには無いが、ソースコードを見てみたところ(ipythonではstats.mstats.linregress??で見ることができる。scipy version 0.7.0）、それぞれの変数の分散の比になっていた。結果は

slope= 0.0666657
intersect= 31.9559306074
correlation= 0.618468
p value= 4.38452479663e-62
std. of residuals 0.62129592706
(Var. Latitude)/(Var. Sea level)= 0.011619

regressionの線も引いてみよう。

 

プログラムはmatplotlibを使って以下のようになる。

#plot
import matplotlib.pyplot as plt
import numpy as np

x=np.arange(-7.,35.,1.) # x data from -7 to 35 step 1
y=a_s*x+b_s         # regression line

plt.figure()
plt.plot(kushimoto_uragami,kuroshio_lat,'.') # plot data
plt.plot(x,y,'r-') # plot regression line
plt.xlabel("Kushimoto-Uragami")
plt.ylabel("Kuroshio Latitude")
plt.xlim((-7.,35.))
plt.ylim((30.,33.5))
plt.savefig("kurosho_reg_scipy.png")
plt.show()

今回はつかわなかったが、以下の情報もある。
Cookbook/Fit Statistics 

Scikits.statsmodels

multiple regression

串本・浦神の潮位差と黒潮流路緯度 Rで線形モデル （続き2）

前回のさらに続き。推定値・予測値の信頼区間も描くことができる。

前々回もとめた係数をα、切片をβとすると、
(1) -7から35を1刻みで分割した串本・浦神の各点x0
(2) αx0+βの推定値、およびその99%信頼区間
(3)  αx0+βの推定値、および値αx0+β+ε0（ε0はx0で混入する誤差）の予測値の99%信頼区間
(4)  αx0+βの推定値および予測値の信頼区間を黒、推定値の信頼区間を青でプロット　（計6本の線。直線ではない）。
(5) グラフを重ねがきする。
(6) データ点をプロットする。

rscript="""
new<-data.frame(rlevel=seq(-7,35,1))  # (1)
png(\'kuroshio_reg_R3.png\')
dc<-predict(fm,new,interval="confidence",level=0.99) #(2)   
dp<-predict(fm,new,interval="prediction",level=0.99) #(3)
matplot(new$rlevel,cbind(dp,dc),lty=1,col=c("black","black","black","black","blue","blue"),type="l",xlab="",ylab="",xlim=c(-7, 35),ylim=c(30,33.5))　 #(4)
par(new=T) #(5)
plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\",xlim=c(-7, 35),ylim=c(30,33.5))  #(6)
dev.off()     
"""
r(rscript)

図

 

参考:　やさしい医学統計手法　6.二標本の関連性を見る。

http://www3.ocn.ne.jp/~stat/medical/med_023.htm

串本・浦神の潮位差と黒潮流路緯度 Rで線形モデル （続き）

前回の続き。
以下のようにして残差(residuals)に関するプロットを作ることができる。
(1) 図を四分割し、上下左右の間隔を指定。現在の作図パラメータはoparに退避。
(2）4種類の診断図を一度に描く。
(3) 作図パラメータをもとにもどす。

r('png(\'kuroshio_reg_R2.png\')')
r('opar<-par(mfrow=c(2,2),oma=c(0,0,1.1,0))')　　#(1)
r('plot(fm)')                                       #(2)
r('par(opar)')                                      #(3)  
r('dev.off()')

図

 

左上：　予測値対残差。残差が1.96σを超える（正規分布で95%内に入っていない）データにはデータ番号が示されている。

右上：　標準化された残差に対する累積比率に対する正規Q-Qプロット。残差が正規分布しているか。

左下：　残差の標準化の絶対値の平方根を予測値に対してプロットしている。

右下：　leverage(梃子比: 回帰直線への影響度）と標準化された残差。

例えば、右上は黒潮の緯度が高い（日本に近い）場合には、回帰直線は過大評価になっている。これは陸地があるので、それ以上は高い位置にいけないからであろう。緯度が低い（沿岸から離れる）場合に値がばらついているのは、黒潮が十分に岸から離れている場合には串本・浦神の潮位差は0に近く、黒潮の離岸距離との関係がうすまっていると考えられる。

串本・浦神の潮位差と黒潮流路緯度 Rで線形モデル

前回の続き。今回はrpy2を通してRを使い、串本・浦神の潮位差と黒潮流路緯度を線形モデルとして評価してみる。
(1) データをRに渡し、scatter plotを描く。
(2) 潮位差(rlevel)と黒潮緯度(rlat)の線形モデルを考える。
(3) 線形モデルをもとに、regression line の直線を引く。
(4) 線形モデルの結果を打ち出す。
missing valueもmasked arrayから自動的に選別して除いてくれているようだ。

#read module
#for data
from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

#for R
import rpy2.robjects as robjects

#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()   


#for R
rlat=robjects.FloatVector(kuroshio_lat)
rlevel=robjects.FloatVector(kushimoto_uragami)
robjects.globalEnv["rlat"] = rlat
robjects.globalEnv["rlevel"] = rlevel
r=robjects.r
r('png(\'kuroshio_reg_R.png\')')    
r('plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\")') #scatter plot 
r('fm<-lm(rlat ~ rlevel)')  # linear model
r('abline(fm)')             # plot regression line
r('dev.off()')            
summary=r('sm<-summary(fm)')
print summary

図

 

summaryの出力

Call:
lm(formula = rlat ~ rlevel)

Residuals:
     Min       1Q   Median       3Q      Max
-3.17592 -0.32259  0.08407  0.40907  1.30406

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 31.955937   0.038044  839.98   <2e-16 ***
rlevel       0.066666   0.003672   18.15   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6225 on 532 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared: 0.3825,    Adjusted R-squared: 0.3813
F-statistic: 329.5 on 1 and 532 DF,  p-value: < 2.2e-16

線形モデルの係数は0.066666、切片は31.955937。
ともにp値は小さく(2e-16)、0であるという帰無仮説は棄却できる。
線形モデルのR^2値は約38%で、前回求めた相関係数0.62の2乗になっている。

係数などは以下のようにとりだせる。

print summary.names

 [1] "call"          "terms"         "residuals"     "coefficients"
 [5] "aliased"       "sigma"         "df"            "r.squared"  
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"  

cf=r('cf<-sm$coefficients')
print cf.names

[[1]]
[1] "(Intercept)" "rlevel"   

[[2]]
[1] "Estimate"   "Std. Error" "t value"    "Pr(>|t|)"

for i in range(8):
    print cf[i]

31.9559366577
0.066665747972
0.0380437732903
0.00367237036394
839.978106637
18.1533291486
0.0
1.13483941203e-57

参考：単回帰http://www.geocities.jp/kashiwabara_fact/statistics-tankaiki.htm

関連ポスト
北極振動　Rで統計
北極振動　Rで統計グラフ編

串本と浦神の潮位差と黒潮流路位置の相関

串本と浦神の潮位差が黒潮の流路のよい指標となることは経験的に知られている（参考:気象庁）。
前々回作った黒潮の流路の緯度を読むプログラムと、前回作った串本と浦神の潮位差を読むプログラムをモジュールとして用いて、潮位差を横軸に、流路緯度を縦軸にプロットしてみた。

 

#read module
from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

import matplotlib.pyplot as plt
import numpy as np

#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()  

#plot
plt.figure()
plt.plot(kushimoto_uragami,kuroshio_lat,'.')
plt.xlabel("Kushimoto-Uragami")
plt.ylabel("Kuroshio Latitude")
plt.savefig("kurosho_corr.png")
plt.show()

潮位差が大きいと、黒潮の緯度は高くなる（岸に近い流路をとる）傾向があることがわかる。
相関をとってみよう。

c1=np.corrcoef([kushimoto_uragami,kuroshio_lat])
print c1 
[[ 1.          0.11969213]
 [ 0.11969213  1.        ]]

相関係数は 0.11969213　!?   小さ！！！

numpyのmasked arrayでマスクした潮位差の欠損値のところでうまく計算できていないらしい。ちゃんと欠損値を抜いて計算してみよう。

mask=(kushimoto_uragami.mask==False)
c2=np.corrcoef([kushimoto_uragami[mask],kuroshio_lat[mask]])
print c2

[[ 1.          0.61846868]
 [ 0.61846868  1.        ]]

今度はそれなりの相関係数 0.61846868 が出た。

実は初めからmasked array用の相関係数の関数ma.corrcoefを使えば、そのままでも計算できる。

c3=np.ma.corrcoef([kushimoto_uragami,kuroshio_lat])
print c3

[[1.0 0.61846868176]
 [0.61846868176 1.0]]

串本と浦神の潮位差

串本と浦神の潮位差は黒潮の流路と密接な関係があることが知られている（参考:気象庁）。串本と浦神の潮位差の月平均値のデータは気象庁から手に入れることができる。http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/ksur.txt

以下はそれを読むプログラム。基本的には前回と同じだが、欠損値があるので一手間が必要。

 

read_kushimoto_uragami.py

import numpy as np
import scikits.timeseries as ts

def read_kushimoto_uragami():
    day=[]
    value=[]
    ds = np.DataSource(None)
    f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/ksur.txt')
    for line in f:
         fields = line.strip().split()
         if len(fields)>=1: day.append(fields[0])
         if len(fields)==1:
             value.append(-999.)
         elif len(fields)==2:
             value.append(fields[1])
    f.close()
    valuea=np.asarray(value,dtype='float32')
    kushimoto_uragami=ts.time_series(valuea,start_date=day[0],freq='M',mask=(valuea<-990.))
    return kushimoto_uragami

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import scikits.timeseries.lib.plotlib as tpl

    kushimoto_uragami=read_kushimoto_uragami()

    fig = tpl.tsfigure()
    fsp = fig.add_tsplot(111)
    fsp.tsplot(kushimoto_uragami, '-')
    plt.xlabel("YEAR")
    plt.ylabel("Kushimoto-Uragami (cm)")
    plt.savefig("kushimoto_uragami.png")
    plt.show()

黒潮の流路変化

*
*

初ガツオ：水揚げ　勝浦漁港、活気づく　／千葉

勝浦市の勝浦漁港に２５日、初ガツオが水揚げされた。江戸中期の俳人・山口素堂が「目には青葉山ほととぎす初鰹」と詠む通り、黒潮に乗って日本列島に近づく初ガツオの本格シーズンは５月ごろ。冬の終わりを予感させる水揚げに、漁港は活気づいた。

（毎日新聞 2010/1/25）

黒潮の流れは、日本南岸を岸に近づいたり、離れたりして流れる。
参考： 黒潮の数か月から十年規模の変動（流路）　気象庁

黒潮の流軸の緯度のデータは上の気象庁のサイトから入手できる。

以下はそのデータを読取りプロットするpythonプログラム。

import numpy as np
import scikits.timeseries as ts

def read_kuroshio_lat():
    ds = np.DataSource(None)
    f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/kuro_slat.txt')
    data=np.genfromtxt(f,dtype=[('date','S7'),('lat','f4')])
    f.close()
#    ds.__del__
    kuroshio_lat=ts.time_series(data['lat'],start_date=data['date'][0],freq='M')
    return kuroshio_lat

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import scikits.timeseries.lib.plotlib as tpl

    kuroshio_lat=read_kuroshio_lat()

    fig = tpl.tsfigure()
    fsp = fig.add_tsplot(111)
    fsp.tsplot(kuroshio_lat, '-')
    plt.xlabel("YEAR")
    plt.ylabel("Kurosho latitude")
    plt.savefig("kurosho_lat.png")
    plt.show()

北極振動時の日本付近の冬の気圧パターン

*
*
12月の1000hpaのgeopotentail heightに前にもとめた北極振動のパターンを重ねてみる。北極振動指数が正だと西高東低が弱まり、負だと強まる。

北極振動指数0の時

 

北極振動+2の時

 北極振動-2の時

  

以下、プログラム

from mpl_toolkits.basemap import Basemap  # for Basemap toolkit

import matplotlib.pyplot as plt
import numpy as np

lat=np.load('lat.npy')  # latitude
lon=np.load('lon.npy')  # longitude
lat[-1]=lat[-1]-1.E-4  # lat=90 does not work

#average field
mhgt=np.load('mhgt.npy') # seasonal climatology
dec_average=mhgt[11,:,:] # december

#AO pattern
factor = np.sqrt(np.cos(np.pi*lat/180.))
eof1=np.load('eof1.npy')
AO_pattern=eof1/factor[:,np.newaxis]


#------------- plot geopotential height (AO index 0, 2,-2)
for index in [0.,2.,-2.]:
    data=dec_average+AO_pattern*index
    fig=plt.figure()
    ax = fig.add_axes([0.1,0.1,0.8,0.8])
    map = Basemap(width=15.e6,height=12.e6,\
            projection='gnom',lat_0=60.,lon_0=135.)
    map.drawcoastlines()

    lons,lats=np.meshgrid(lon[:],lat[:])
    x,y=map(lons,lats)
    cy=map.contourf(x,y,data,np.arange(-100,300,20),extend='both')
    cz=map.contour(x,y,data,np.arange(-100,300,20),colors='k')
    plt.clabel(cz,fmt='%3.0f',)
    plt.savefig('AO_patter_'+str(index)+'.png')
    plt.show()

北極振動指数の時系列 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()

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

*
*
前回の続き。前回もとめたfirst EOFの時系列が、NOAA CPCの発表している北極振動の値と一致しているか確かめる。黒がNOAAのCPCの値で、赤が自分でもとめたもの。
ちょっとだけ値がずれているが、ほとんどあっているのでいいんじゃないかな。前回よくわからなかった北極のとりあつかいとかがずれる原因かもしれない。

以下がプログラム。以前作ったreadAOモジュールでNOAA CPCの北極振動を読む。時系列オブジェクトを作るのと、プロットにはscikits.timeseriesを用いた。

# plot pc1
import scikits.timeseries as ts
import scikits.timeseries.lib.plotlib as tpl

dbeginmonth=ts.Date('M', '1979-1')
pc1_series=ts.time_series(pc1,start_date=dbeginmonth,freq='M')

from readAO import readAO
AOindex_series=readAO()
   
fig = tpl.tsfigure()
fsp = fig.add_tsplot(111)
l1=fsp.tsplot(AOindex_series, 'k-')
l2=fsp.tsplot(pc1_series, 'r-')
istart=ts.Date('M','1979-1').value
iend=ts.Date('M','2000-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.png")
plt.show()

北極振動 パターン

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