熱帯赤道上空東西風についてのEOF解析（QBOのindex化）

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import xarray as xr
import datetime

from pandas.plotting import register_matplotlib_converters #これを描かないと警告文がでる
register_matplotlib_converters()
In [2]:
# データの読み込み
# 以下を参考にしました
# https://arg.usask.ca/docs/LOTUS_regression/_modules/LOTUS_regression/predictors/download.html

file = "data/qbo_1953_2021.dat"
#--> "data/dd/qbo_1953_2020.dat"に変更のこと

def date_parser(s):
        s = int(s)
        return datetime.datetime(2000 + s // 100 if (s // 100) < 50 else 1900 + s // 100, s % 100, 1)

dataORG = pd.read_fwf(file,
                       skiprows=92, header=None,
                       colspecs=[(0, 5), (6, 10), (12, 16), (19, 23), (26, 30), (33, 37), (40, 44), (47, 51), (54, 58)],
                         delim_whitespace=True, index_col=1, parse_dates=True, date_parser=date_parser,
                         names=['station', 'month', '70', '50', '40', '30', '20', '15', '10'])

#最初の方（1950年代は歯抜けが多いので、skiprows=92として読み飛ばしています）

#dataORG.index = dataORG.index.to_period(freq='M')

data_sub = dataORG.loc['1960-01':'2018-12',:].copy()
data_sub.drop('station', axis=1, inplace=True)
data_sub /= 10.0 # units: m/s に変換


print(data_sub.head()) #列方向は気圧座標が入っている（高度に換算すると 70 hPa~20 km, 10 hPa~ 30 km, に相当）

             70   50   40    30    20    15    10
month                                            
1960-01-01  6.5  5.0  3.0  -0.2  -1.8  -7.4 -14.0
1960-02-01  2.0  6.7  5.4  -1.8  -5.1  -9.0 -17.0
1960-03-01  5.0  7.0  5.0  -4.6  -8.3 -11.0 -19.0
1960-04-01  6.5  6.0  3.0 -12.4 -15.2 -16.0 -20.0
1960-05-01  8.0  6.9 -0.5 -18.2 -17.0 -18.0 -23.0

その他のヒント・補足

  • Deseasonalization（季節変化の除去）について：
    • 色々やり方はあるが、各月でコンポジットを取って月平均の東西風（の高度分布）を算出し、オリジナルデータから差し引くのが簡単。

  • 行列の積を計算する際に、1次元の縦ベクトルと横ベクトルを区別したいときがある。uを大きさnの配列とする

    • uh = np.reshape(u, (1,n)):横ベクトルにする
    • uv = np.reshape(u, (n,1)):縦ベクトルにする
In [3]:
# Fig 1 (upper)

fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)

data_plot = data_sub.loc['1960':'1990','30'].copy()
#data_plot.plot(ax=ax)
x = data_plot.index
ax.scatter(x,data_plot)
plt.show()
In [ ]:
In [4]:
# (準備): Deseasonalization（季節変化の除去）

# その１(大畠 くん) こちらがシンプルで良いですね

#month = data_sub.index.month
#data_ds = data_sub.copy()
#for im in range(1,13):
#    data_ds[month == im] = data_ds[month == im] - data_sub[month == im].mean(axis=0)


key = data_sub.index.month
monthly = data_sub.groupby(key).mean()

data_ds = data_sub.copy()


print(data_ds.shape)

m,n = data_sub.shape

for i in range(0,m):
    mm = i % 12 
    data_ds.iloc[i,:] = data_sub.iloc[i,:] - monthly.iloc[mm,:]

    
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)

data_sub.iloc[:,3].plot(ax=ax,label='Original')
data_ds.iloc[:,3].plot(ax=ax,label='Deseasonalized')

ax.legend()
plt.show()

(708, 7)
In [5]:
# 値だけ取り出す
data_ds_v = data_ds.values

#
#mean = data_ds.mean(axis=0) #　平均
#anom = data_ds - mean #　平均偏差行列

# 特異値分解(Numpyの関数を使う)
U, s, Vh = np.linalg.svd(data_ds_v)

m,n = data_ds_v.shape
In [6]:
#Fig. 1（EOF1, 2モードの鉛直構造）

EOF1 = s[0]*Vh[0,:]/m**0.5
EOF2 = s[1]*Vh[1,:]/m**0.5

pres_s = data_sub.columns #気圧座標値を取得

pres_v = list(map(float,pres_s))

fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)

ax.plot(EOF1,pres_v,label="EOF1")
ax.plot(EOF2,pres_v,label="EOF2")
ax.legend()


ax.set_xlabel('Zonal wind (m/s)')
ax.set_yscale('log')
ax.invert_yaxis()

plt.show()
In [7]:
#Fig. 2

pc1 = U[:,0]*m**0.5
pc2 = U[:,1]*m**0.5

x_date = data_sub.index #日付座標を取得
#print(x_date)

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111)

ax.plot(x_date,pc1,label="PC1")
ax.plot(x_date,pc2,label="PC2")

ax.legend()

plt.show()
In [8]:
#Fig.3

fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)

ax.scatter(pc1,pc2)

ax.set_xlim(-2,2); ax.set_ylim(-2,2)

ax.set_xlabel('PC1'); ax.set_ylabel('PC2')

plt.show()
In [9]:
#Fig.4

reconst = np.zeros((m,n))
for idx in range(0,2):
    vh_s = np.reshape(Vh[idx,:], (1,n))
    u_s = np.reshape(U[:,idx], (m,1))

    reconst = reconst + s[idx]*(u_s@vh_s)

fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.contourf(x_date,pres_v,data_ds.T)
ax2.contourf(x_date,pres_v,reconst.T)

ax1.set_title('Original')
ax2.set_title('Reconstructed with EOF1+2')

ax1.set_yscale('log'); ax1.invert_yaxis()
ax2.set_yscale('log'); ax2.invert_yaxis()

fig.tight_layout()

plt.show()

↑２つのモードで全体の特徴を良く表せている

ここから発展課題（東西平均東西風データとの回帰係数を計算する）

In [10]:
file = "Fortran/data/uwnd.mon.mean.nc" 
#月平均東西風データを指定
#"data/ncep.reanalysis.derived/pressure/"以下に変更

data = xr.open_dataset(file)

print(data)

<xarray.Dataset>
Dimensions:  (lat: 73, level: 17, lon: 144, time: 855)
Coordinates:
  * level    (level) float32 1000.0 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2019-03-01
Data variables:
    uwnd     (time, level, lat, lon) float32 ...
Attributes:
    title:          monthly mean u wind from the NCEP Reanalysis
    description:    Data from NCEP initialized reanalysis (4x/day).  These ar...
    platform:       Model
    Conventions:    COARDS
    NCO:            20121012
    history:        Mon Jul  5 22:36:33 1999: ncrcat uwnd.mon.mean.nc /Datase...
    References:     http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reana...
    dataset_title:  NCEP-NCAR Reanalysis 1
In [11]:
uwnd = data['uwnd'].sel(time=slice('1960-01','2018-12')).mean(dim='lon')

plevel = uwnd['level']
lat = uwnd['lat']
time = uwnd['time']
In [12]:
uwnd_v = uwnd.values

# 12月データだけ抜き出す
mask = (time.dt.month == 12)  
uwnd_1 = uwnd_v[mask]
uwnd_1_anom = uwnd_1 - uwnd_1.mean(axis=0)

pc1_1 = pc1[mask]
pc2_1 = pc2[mask]


nt,nla,nlo = uwnd_1_anom.shape

#回帰係数を計算
a1 = (uwnd_1_anom*pc1_1[:,np.newaxis,np.newaxis]).sum(axis=0)/(pc1_1**2).sum()
a2 = (uwnd_1_anom*pc2_1[:,np.newaxis,np.newaxis]).sum(axis=0)/(pc2_1**2).sum()
In [13]:
#Fig.5 PC2との回帰
fig = plt.figure(figsize=(6,6))

ax = fig.add_subplot(111)

cf = ax.contourf(lat,plevel,a2)

cbar = plt.colorbar(cf,orientation='horizontal')
cbar.set_label('Zonal wind [m/s]')

ax.set_yscale('log')
ax.set_ylabel('Pressure (hPa)')
ax.invert_yaxis()

plt.show()
In [15]:
# QBOでどれだけ説明できるか？

var_tot = uwnd_1.var(axis=0)

ratio = (a1**2+a2**2)/var_tot

fig = plt.figure(figsize=(6,6))

ax = fig.add_subplot(111)
cf = ax.contourf(lat,plevel,ratio)

cbar = plt.colorbar(cf,orientation='horizontal')
cbar.set_label('Ratio')

ax.set_yscale('log')
ax.set_ylabel('Pressure (hPa)')
ax.invert_yaxis()

plt.show()

1を超えるのは、おそらくインデックスが直交しているとは限らないから（１２月で指定しているため）

In [ ]: