3. 北極振動(Arctic Oscillation: AO) の解析

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.util as cutil
import cartopy.feature as cfeature
import datetime
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

◆ データの読み込み・描画サンプル

以下のサンプルプログラムでは、データの読み込みをxarrayで行い、ndarray型に変換した上で、描画にはmatplotlib (+cartopy)を用いている。

（注意）

前期の演習DBでは、データの読み込みから描画まで xarrayのメソッドを用いて行った。 これは、幸い凝った解析をしなかったために一連の操作を「xarrayのメソッドだけ」で 対応できたからである（ndarrayに変換する必要がなかった）。

参考）xarrayについては演習DB(3.全球大気再解析データの解析参照。

In [2]:
file = "../DD/Fortran/data/slp.mon.mean.nc"
#-->"data/ncep.reanalysis.derived/surface/slp.mon.mean.nc"

dataset = xr.open_dataset(file) # データセットの読み込み（データ変数＋軸）
data = dataset['slp'] # 海面気圧("slp")読み込み
lat = data['lat']
lon = data['lon']

#描画サンプル、適当な年月データを切り出して描く
data1 = data.sel(time = '2000-01').squeeze()

# ここで軸を除いてデータ行列(ndarray)にする:
data1_v = data1.values
lat_v = lat.values
lon_v= lon.values

# 描画ここから-------

# 経度方向のデータをサイクリックにする（本演習ではここまでしなくても良いが）
data_cyclic_v, lon_cyclic_v = cutil.add_cyclic_point(data1_v,coord=lon_v)

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111,projection=ccrs.NearsidePerspective(central_latitude=90,central_longitude=180))
# オプションでcaropyの地図描画方法を指定

levels = np.linspace(970,1050,17) #コンターレベル設定 (例: 970から1050まで17レベル)
cf = ax.contourf(lon_cyclic_v,lat_v,data_cyclic_v,transform=ccrs.PlateCarree(),levels=levels,cmap='jet') 

# カラーバー
cbar = plt.colorbar(cf) #,orientation='horizontal')
cbar.set_label('SLP [hPa]')

ax.gridlines()  #罫線
ax.coastlines() #海岸線 
ax.set_global() #オマジナイ

plt.show()

◆ 月平均の計算について(課題演習DBの復習）

時間軸の"月"情報は以下のようにして抜き出せる。後は前回(QBO解析)と同様に、、、

In [3]:
key = data['time'].dt.month #"dt"をつけることで、時刻型として認識される
print(key.head())

<xarray.DataArray 'month' (time: 5)>
array([1, 2, 3, 4, 5])
Coordinates:
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 1948-05-01

◆ ブロードキャストについて

（"operands could not be broadcast together with shapes"などというエラーが出たら...）

　形状の異なる２つの配列間での演算（ブロードキャストという）においては注意が必要である。 ここでは例えば、データ行列($X=X(x,y,t)$) に緯度($y$)に応じた面積重み($w=w(y)$)を乗じる演算がそうである。

　xarray 型で演算を行えば、共通の軸（上の例では緯度）を認識してくれるので気にする必要は無い（気にしなくてもウマく演算してくれる）。一方で、ndarray型での演算においては、少々工夫が必要となる。以下の簡単な例を参考に。

（※さらに詳しいことを知りたい人は、こちらの「例題２」（2019年度課題演習DB資料）　などを参考にしてください。）

In [4]:
A = np.zeros((2,3,4))
B = np.zeros(3) + 1

print(A + B) #型が違うのでエラーがでる

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-01f8fa6438c4> in <module>
      2 B = np.zeros(3) + 1
      3 
----> 4 print(A + B) #型が違うのでエラーがでる

ValueError: operands could not be broadcast together with shapes (2,3,4) (3,)
In [5]:
B2 = B[np.newaxis,:,np.newaxis] # 形が合うように、新しく軸を挿入する
print(B2.shape)
print(A + B2)

(1, 3, 1)
[[[1. 1. 1. 1.]
  [1. 1. 1. 1.]
  [1. 1. 1. 1.]]

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

となって、計算できていることがわかります。

◆ 特異値解析ルーチンについて補足

　これまで実践してきたように、特異値分解（ルーチン）は二次元データ配列についての処理なので、 ここでは元の三次元データ行列(nt, nlat, nlon)を、予め二次元データ行列(nt, nlat$\times$nlat)に変形して おく必要がある。これには、np.reshape()を用いるのが便利（以下の例を参考に）。

　もちろん、結果を描画にする際には再び元の形状に戻しておくこと。

In [6]:
A = np.zeros((2,3,4))
print("original size = ", A.shape)

Ar = A.reshape(2,3*4)
print("New size = ", Ar.shape)

original size =  (2, 3, 4)
New size =  (2, 12)
In [ ]: