import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import random
m = 2000
#r1 = [random.normalvariate(mu=0,sigma=1.0) for i in range(0,n)]
#r2 = [random.normalvariate(mu=0,sigma=1.0) for i in range(0,n)]
factor = np.sqrt(3)*2.0 #標準偏差が1になるように規格化
r1 = 1.5*factor*np.array([random.random()-0.5 for i in range(0,m)])
r2 = factor*np.array([random.random()-0.5 for i in range(0,m)])
%matplotlib notebook
# (1) AOが真(Eq. 9.1)
xa = -r1 + r2 #x軸
xp = -r1 - r2 #y軸
xn = 2*r1 #z軸
xanom = np.vstack([xa,xp,xn]).T
print(xanom.shape)
# (2) AOが偽(NAO-PNA) (Eq. 9.2)
ya = -r1
yp = -r2
yn = r1 + r2
yanom = np.vstack([ya,yp,yn]).T
# 特異値分解
xU, xs, xVh = np.linalg.svd(xanom)
yU, ys, yVh = np.linalg.svd(yanom)
# 描画
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121,projection='3d')
ax2 = fig.add_subplot(122,projection='3d')
ax.scatter(xa,xp,xn,alpha=0.3,s=5)
ax2.scatter(ya,yp,yn,alpha=0.3,s=5)
for i in range(0,2):
x2 = xs[i]*xVh[i,0]/m**0.5; y2 = xs[i]*xVh[i,1]/m**0.5; z2 = xs[i]*xVh[i,2]/m**0.5
ax.quiver([0],[0],[0],[x2],[y2],[z2],color='red') #主成分ベクトル(AO真)
print("EOF" + str(i+1) +" for AO-NCM: ",[x2,y2,z2])
x2 = ys[i]*yVh[i,0]/m**0.5; y2 = ys[i]*yVh[i,1]/m**0.5; z2 = ys[i]*yVh[i,2]/m**0.5
ax2.quiver([0],[0],[0],[x2],[y2],[z2],color='red') #主成分ベクトル(AO偽)
print("EOF" + str(i+1) +" for NAO-PNA: ",[x2,y2,z2])
# AO, PNAの基底ベクトル
ax.quiver([0],[0],[0],[-1],[-1],[2],color='orange')
ax.quiver([0],[0],[0],[1],[-1],[0],color='orange')
# NAO, PNAの基底ベクトル
ax2.quiver([0],[0],[0],[-1],[0],[1],color='orange')
ax2.quiver([0],[0],[0],[0],[-1],[1],color='orange')
ax.set_xlabel('Atlantic (A)'); ax.set_ylabel('Pacific (P)'); ax.set_zlabel('North Pole (N)')
ax2.set_xlabel('Atlantic (A)'); ax2.set_ylabel('Pacific (P)'); ax2.set_zlabel('North Pole (N)')
ax.set_xlim(-4,4);ax.set_ylim(-4,4);ax.set_zlim(-4,4)
ax2.set_xlim(-4,4);ax2.set_ylim(-4,4);ax2.set_zlim(-4,4)
plt.show()
xanom = np.vstack([xa,xp]).T
yanom = np.vstack([ya,yp]).T
# 特異値分解(Numpyの関数を使う)
xU, xs, xVh = np.linalg.svd(xanom)
yU, ys, yVh = np.linalg.svd(yanom)
# 描画
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax.scatter(xa,xp,s=10)
ax2.scatter(ya,yp,s=10)
for i in range(0,2):
x2 = xs[i]*xVh[i,0]/m**0.5; y2 = xs[i]*xVh[i,1]/m**0.5
ax.quiver([0],[0],[x2],[y2],color='red',label ="EOF"+str(i))
print("EOF" + str(i+1) +" for AO-NCM: ",[x2,y2])
x2 = ys[i]*yVh[i,0]/m**0.5; y2 = ys[i]*yVh[i,1]/m**0.5
ax2.quiver([0],[0],[x2],[y2],color='red',label ="EOF"+str(i))
print("EOF" + str(i+1) +" for NAO-PNA: ",[x2,y2])
# NAO, PNAの基底ベクトル
ax.quiver([0],[0],[-1],[-1],color='orange')
ax.quiver([0],[0],[1],[-1],color='orange')
# NAO, PNAの基底ベクトル
ax2.quiver([0],[0],[-1],[0],color='orange')
ax2.quiver([0],[0],[0],[-1],color='orange')
ax.set_xlabel('Atlantic (A)'); ax.set_ylabel('Pacific (P)')
ax2.set_xlabel('Atlantic (A)'); ax2.set_ylabel('Pacific (P)')
ax.set_xlim(-4,4);ax.set_ylim(-4,4)
ax2.set_xlim(-4,4);ax2.set_ylim(-4,4)
plt.show()
二点EOFの場合は、AOが"真"の場合(左)でも"偽"の場合(右)でも
ことがわかる(符号の不定性は除いて)。
これが、AOの実在性を調べる為に中緯度EOFを施す
根拠となっている。