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

データの読み込み部分のみ以下に示す

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr

from pandas.plotting import register_matplotlib_converters #これを描かないと警告文がでる（オマジナイ）
register_matplotlib_converters()

file = "data/radiosonde_tropical_eastward_wind_195301-202509.nc"
# --> "data/dd/radiosonde_tropical_eastward_wind_195301-202509.nc"　に変更

data_org = xr.open_dataset(file)["u"].sel(time=slice("1960","2018"),pressure=[70,50,40,30,20,15,10])
print(data_org)

/usr/local/lib/python3.9/dist-packages/gribapi/__init__.py:23: UserWarning: ecCodes 2.31.0 or higher is recommended. You are running version 2.20.0
  warnings.warn(

<xarray.DataArray 'u' (time: 708, pressure: 7)>
[4956 values with dtype=float64]
Coordinates:
  * pressure  (pressure) int32 70 50 40 30 20 15 10
  * time      (time) datetime64[ns] 1960-01-01 1960-02-01 ... 2018-12-01
Attributes:
    standard_name:  eastward_wind
    long_name:      Eastward Wind
    units:          m/s
    cell_methods:   time: mean (interval: 1 month)
    u_qc:           0

time = data_org["time"]
pressure = data_org["pressure"]

time_v = time.values
pressure_v = pressure.values
data_org_v = data_org.values 
# 後々のために、データ行列（数値配列）にしておく。


# サンプル描画　その１：ラインプロット
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)
ax.plot(time_v,data_org_v[:,0],label='Original') # "0"なので、pressure levelの最初の層である 70 hPa面
ax.legend()
plt.show()

# サンプル描画　その２：時間ー高度断面
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)
clb = ax.contourf(time_v,pressure,data_org_v.T) 
# 横軸が時刻、縦軸が気圧になるように、データは転置している（.T）

ax.set_yscale('log'); ax.invert_yaxis()
ax.set_xlabel('Date'); ax.set_ylabel("Pressure [hPa]")
plt.colorbar(clb)
plt.show()

Deseasonalization（季節変化の除去）について：

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

• 「月」の情報はmonth = data["time"].dt.monthなどとすれば得られます。

---

解析にあたってのヒント・補足

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

  • uh = np.reshape(u, (1,n)):横ベクトルにする
  • uv = np.reshape(u, (n,1)):縦ベクトルにする

---

描画方法のヒント・補足

• (1): EOF1 & EOF2 を 気圧座標（縦軸）の関数として描画

• (2): PC1 & PC2 を 時刻（横軸）の関数として描画

• (3): PC1 (横軸）➖ PC2（縦軸）の散布図

  • 散布図を描くメソッド: ax.scatter(x,y,....)

• (4): 横軸：時刻(x)、縦軸(y):気圧座標として、東西風速(z)をコンター表示

  • コンター（塗りつぶし）メソッド: ax.contourf(x,y,z) (z = z(y, x))

• (5:発展): 横軸: 緯度(x)、縦軸(y):気圧座標として、回帰係数をコンター表示

---

発展課題（5）について

• 月平均東西風データはdata/ncep.reanalysis.derived/pressure/uwnd.mon.mean.ncを使用。

• 以下、データの読み込み例（の一部）を示す:

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

data = xr.open_dataset(file)

print(data)

<xarray.Dataset>
Dimensions:  (level: 17, lat: 73, lon: 144, time: 855)
Coordinates:
  * level    (level) float32 1e+03 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

uwnd = data['uwnd'].sel(time=slice('1960-01','2018-12')).mean(dim='lon')

plevel = uwnd['level']
lat = uwnd['lat']
time = uwnd['time']

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