3.　全球大気再解析データ（４次元データ）の解析¶

3.0 はじめに¶

3.1 データの読み込み：NetCDFデータを読み込む¶

3.2 データを切り出す¶

3.2 統計解析・時系列解析¶

3.4 絵を描く¶

$\ \ \ \ $ 3.4.1 二次元プロットの方法¶

$\ \ \ \ $ 3.4.2 地図描画¶

3.5 【参考】ベクトル場プロット¶

3.0 はじめに¶

　この章では大気再解析 (Reanalysis)を用いた解析を行う。前章までと大きく異なるのはデータの次元が一気に増えたことである（「時間１次元」->「時間１次元＋空間３次元」）。したがって単純に考えれば、どの軸に沿った図を描くかによって、ラインプロットでは４通り、二次元平面図では$_4\text{C}_2 = 6$通りの選択肢がある。

　再解析データを含むグローバルデータセットの多くはNetCDF形式で提供されている。これは自己記述式と言われるもので、データのヘッダーにデータに関する様々な情報が記載されている。

　ここではNCEP/NCAR再解析を使って、NetCDF形式データの読み込みと情報の縮約、そして二次元データの様々な描画方法について実習を行う。再解析データでは様々な物理量のデータが提供されているが、本チュートリアルでは月平均東西風データ (uwnd.mon.mean.nc)を用いることにする（ベクトル場では、これに加えてvwnd.mon.mean.ncも用いる）。

　この章で新たに用いるライブラリは以下の通り：

cartopy: 地図描画ライブラリ (https://scitools.org.uk/cartopy/docs/latest/)
xarray: 多次元配列解析ライブラリ(http://xarray.pydata.org/en/stable/index.html)

XArrayについて¶

　　xarrayはpandasをベースにした多次元配列解析ライブラリである。pandasでは一次元(Series型）あるいは二次元(DataFrame型）までしか扱えなかった（注）が、これを多次元に拡張したものと考えることができる。

（注）正確にはDataFrameのインデックスに階層構造を持たせることで扱う次元を増やすこともできる（階層型インデックス)が、少し面倒。

　　XArrayを用いる最大の利点は、データ配列に軸（インデックス）の情報が結合したオブジェクトのまま解析できる点である。ここで軸とは経度・緯度・高度・時間の情報である。軸情報を保持しているので、２章で紹介したようにメソッドを用いて時系列解析や描画を行うことができる。

...さらに凝った解析や描画を行うには...¶

　一方で、xarrayのメソッドで行える解析・描画には限界がある。アドバンストな解析・描画をするには、xarrayで読み込んだデータを数値配列(ndarray)に変換して、(ndarray型に定義された）関数を用いて解析や描画を行えばよい。例えば、XArrayのオブジェクトについては.valuesの属性を参照することで数値データ（ndarray型）を取り出せる（軸情報は落とされる）。3.5節のベクトル場プロットはその一例である。

以下、3.1-3.4において

【練習課題】 2010-2018年平均の日本上空(500hPa)の東西風の分布（緯度-経度断面）を図示する¶

ことを目標に順を追って解説する。

3.1 データの読み込み：NetCDFデータを読み込む¶

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import cartopy.util as cutil
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
#from pandas.plotting import register_matplotlib_converters
#register_matplotlib_converters()

file = "../../DB_2020/data/ncep.reanalysis.derived/pressure/uwnd.mon.mean.nc" # --> 【パスを適宜変更すること！！！！】 
#月平均東西風データ

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

csvファイルと異なり(1節）、データの読み込みが圧倒的に楽であることがわかる。

データ配列だけでなく、軸情報やその他の属性も含むオブジェクト(Dataset型）になっていることが分かる。

これから以下の情報が読み取れる：

次元(Dimensions): （緯度:73, レベル: 17, 経度: 144, 時間: 855)
軸(Coordinates)の値
- レベル(level)： 1000, 925, 850, 700, ...
- 緯度(lat): 90, 87.5, 85.0, ...
- 経度(lon): 0, 2.5, 5.0, ...
- 時間(time): 1948-01-01, 1948-02-01,...
データ変数 (Data variables): uwnd
データ属性（Attributes): ...

この情報を確認した上で、まずは必要な変数を切り出しておく（uwnd = data["uwnd"])

uwnd = data["uwnd"] #変数 "uwnd"の切り出し

3.2 データを切り出す¶

3.2.1 ラベルで範囲指定¶

(1) 軸の名前と値を指定した切り出し: .sel(), .isel()　（お勧め！）¶

.sel(軸1の名前=軸1の値,軸2の名前=軸2の値,method=,..) (.isel()の場合は、（値でなく）整数インデックスで指定)
- method: None（一致するもののみ；デフォルト）、nearest（指定した点に一番近い点を取ってくる）
- 範囲で指定する際にはslice(start,stop,step)を用いる。ただし、データの並び（降順か昇順か）には注意。この例では、緯度と気圧レベルは降順になっているので（上の"Coordinates"を見ればわかる）、例えば30N-50Nの領域を抜き出す際には、lat=slice(50,30)とする必要がある。

(1)と違って軸の並びを気にしなくて良いのでオススメ。

(2) 各軸のラベル・整数のインデックス値を用いた切り出し：.loc[ ]、[ ]¶

　pandasにおける切り取りと同じ。各軸の範囲をラベル（日付、経度など）もしくはインデックス位置を用いて指定する。

　軸の並びに注意すること:ここでは(time, level, lat, lon)の順に並んでいる。

3.2.2 条件指定して切り出し¶

pandasと同様に、条件を指定してデータを切り出すこともできる。

◆ 時刻データを元に条件指定する際の注意点¶

　時刻軸のデータは（上の軸情報を参照して）

index = uwnd['time']

などとして取得できるが、ここから月・日などの情報抜き出すにはさらに.dtをつける：

month = index.dt.month　（月の情報を抜き出す）

(参考）DateimeIndex属性の一覧はこちら

【補足】 .squeeze() メソッド: 大きさ１の次元を削除¶

　ある断面でデータを切り出したり、軸に沿った統計操作を行ったりした際、大きさ１の次元が残ってしまう場合がある(e.g., A(1,1,36,36))。このままの型の配列ではうまく二次元描画できないので、A_sub = A.squeeze()として、大きさ１の次元を削除する必要がある。

# データの切り出し（2010-2018年, 500 hPaレベル）
uwnd_sub = uwnd.sel(time=slice('2010','2018'),level=500)

3.3 統計解析・時系列解析¶

　時系列データの解析では、時間方向の処理を仮定していた。４次元データの解析では、どの軸に沿った処理を行うかに常に注意して解析を行う必要がある（例えば「平均」するにしても、時間平均？ or 東西平均？）。

◆ 統計解析¶

　xarrayの型にも、pandasと同様に統計量を求めるメソッドが定義されている。

.mean(), .std()などのメソッドがある。
ただし、処理を行う軸をaxis=(次元の番号）、もしくは、dim=(次元の名前)としてオプションで指定する（指定しなければ全平均の値が出力される）。例えば、data.mean(dim="lon")（経度平均）やdata.mean(dim=("lon","lat"))（経度・緯度平均）など。

◆ 時系列データ解析¶

　pandasと同様に、reshapeやgroupbyなどのメソッドを使うことで時系列操作も簡単に行える。ただし、処理を行う軸を明示する必要がある:

・resampleの例（年平均データの作成）：¶

uwnd.resample(time="A-DEC").mean(dim="time")

・groupbyの例（各月のコンポジット）：¶

3.2.2を参考に、key = uwnd["time"].dt.monthとしておいて、

uwnd.groupby(key).mean(dim="time")

あるいは、

uwnd.groupby("time.month").mean(dim="time")

と簡単に書くこともできる。

◆ さらに凝った解析を行うには？？¶

　xarrayのメソッドでは事足らない場合には、.valuesとしてndarray型に変換して適宜解析を行えばよい。

#時間平均をとる
uwnd_mean_t = uwnd_sub.mean(dim=('time'))

3.4 データの二次元描画¶

　最後にデータの二次元描画を演習する。最初に述べたようにどの断面で切り出すかを常に意識しながら描画すること。

3.4.1 二次元描画：.plot.contourf（塗りつぶし）, .plot.contour（コンター）¶

.plot.contour,.plot.contourf(ax=,levels=,cmap=,colors=,...):
- ax=: 用いるサブプロット領域
- levels: レベル。(1) 整数で指定した場合：レベル数。(2) 配列で指定した場合：レベルの値を陽に指定。
- cmap: カラーマップ(e.g., "gray"（白黒）、"bwr"（青-赤)、"rainbow"など。詳細はこちら
- colors: コンター、もしくは、塗りつぶしの色。リストで渡すと各レベルの色を陽に指定する。
- cbar_kwargs={}: カラーバーの調整オプション。例：'orientation': 'horizontal'（横向きにする）; 'label': 'hogehoge'（カラーバーのラベルをhogehogeにする）。なお、カラーバーはデフォルトで描いてくれる（素晴らしい！）。
- 詳細はこちら
- ...

【コンターラベルの作成】¶

.plot.contour() メソッドの返り値 を ax.clabel() に渡すことで等高線のラベルを作成できる。

ax.clabel(cntr, fmt=, fontsize=): コンターラベルを作成
- cntr: contour()の返り値
- fmt=: フォーマット(e.g., "%f3.1", "%3i")
- fontsize=: フォントサイズ
- ...

3.4.2 地図描画¶

Cartopyというライブラリを使って、地図投影を行う。

最低限、サブプロット領域の作成(ax = fig.add_subplot(...)、および、図の描画(.contourf(...))の段階で、地図投影オプションを指定することで描画できる（cartopyというライブラリを呼び出している）。

◆ 作成の手順（通常の描画と違う点を中心に）¶

ax（サブプロット領域）のcartopy版を作る: ax = fig.add_subplot(projection=)
- projection=のオプションを入れることで地図描画モードになる（入れなければ通常の絵）。
- projectionの選択肢:
  - ccrs.PlateCarree(central_longitude=): 正距円筒図法（緯度・経度をそれぞれ地図の縦・横に読み替えた円筒図法）
  - ccrs.NearsidePerspective(central_latitude=,central_longitude=): 宇宙から見た鳥瞰図
  - ...
  - 地図投影の一覧はこちら
図を描く: ax.contourf(...,transform = ccrs.PlateCarree()) etc..
- transform=: データ自身の座標系を明示する。再解析データは正距円筒図法でデータが入っているのでccrs.PlateCarree()とする。詳しい（英語の）説明はこちら
- その他はmatplotlibのときと同じ
[オプション] 目盛りの設定: ax.set_xticks(xticks,crs=), ax.set_yticks(yticks,crs=)
- xticks,yticks: 目盛りの位置
- crs=：　目盛り自身の座標系（緯度経度で指定するときはccrs.PlateCarree()とすればよい; 4.も参照）
- ax.xaxis.set_major_formatter(lon_formatter): 目盛りのラベルを描いている。lon_formatterは地図描画用のの書式($^\circ$Eとか）。
[オプション] 描画領域の設定: ax.set_extent([x0,x1,y0,y2],crs=)
- [x0, x1, y1, y2]: データを描く部分領域を指定する
- crs=: 上記のリストの値の座標系（緯度経度で指定するときはccrs.PlateCarree()とすればよい; 4.も参照）
[オプション] その他の情報を付加
- 海岸線: ax.coastlines(resolution=)
- 罫線: ax.set_gridlines()
- 陸域追加: ax.add_feature(cartopy.feature.LAND)
- 海洋域を追加: ax.add_feature(cartopy.feature.OCEAN)
- 陸地の形状を追加: ax.stock_img()
- ...
[全球データを使う場合には] オマジナイ: ax.set_global()

（注意） 3 (目盛りの設定）と4（描画領域の設定）の順序を変えてはいけない！変えると、4（範囲）が3によってにリセットされてしまうようである（私はこれを理解するのに大分苦しみました）。

### 最後に地図上に描画しましょう。

clevels = np.arange(-20,50,5) # 塗りつぶし用

# 描画
fig = plt.figure(figsize=(8,5))

# (1) 領域を用意（地図投影オプションをつけて）

# PlateCaree: 正距円筒図法, central_longitude: 図の中心の経度
ax = fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=180))

# (2) 図を描く
uwnd_mean_t.plot.contourf(ax=ax,transform=ccrs.PlateCarree(),levels=clevels,cmap='rainbow',cbar_kwargs={'label':'U-velocity (m s$^{-1}$)'})
# transform = にはデータ自身の座標系を指定する; 大抵の場合ccrs.PlateCarree()としておけば問題ない

# (3) 目盛り
xticks = np.arange(13)*30
yticks = -90 + np.arange(7)*30
ax.set_xticks(xticks,crs=ccrs.PlateCarree())
ax.set_yticks(yticks,crs=ccrs.PlateCarree())
##　地図投影用の書式を設定（演習では無視（あるいはコピー）してもらって良い）
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)


# (4) 描画領域の設定
ax.set_extent([90,180,20,70],crs=ccrs.PlateCarree())


# (5) オプション
ax.coastlines() #海岸線
ax.gridlines(draw_labels=False) #罫線：ラベルはすでに上で描いたので"False"

plt.show()

別の地図投影も試してみよう（宇宙から見た図）¶

# 描画
fig = plt.figure(figsize=(5,5))

# (1) 領域を用意（地図投影オプションをつけて）

# NearsidePerspective: 宇宙から見た絵
ax = fig.add_subplot(111,projection=ccrs.NearsidePerspective(central_latitude=40,central_longitude=140))

# (2) 図を描く
uwnd_mean_t.plot.contourf(ax=ax,transform=ccrs.PlateCarree(),levels=clevels,cmap='rainbow',\
                          cbar_kwargs={'label':'U-velocity (m s$^{-1}$)'})

### 目盛りはかけない。また全球データを使うので、範囲指定もできない###

# (5) オプション
ax.coastlines() #海岸線
ax.gridlines(draw_labels=False) #罫線：ラベルはすでに上で描いたので"False"

ax.set_global() # 全球データをつかうときのオマジナイ

ax.set_title('U-wind at 500 hPa')

plt.show()

【練習課題（1）】2010-2018年データを用いて、時間・東西平均した東西風を緯度-高度（気圧）断面に描いてみよ（一行目の処理を埋めよ）¶

【練習課題（２）】各月（例えば、１月と７月）の平均値を用いて（１）と同様の図を描いてみよ。（「月」の情報の取り出し方は、3.2.2を参照）¶

uwnd_mean_xt = # ここを各自埋めよ



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

# コンターレベルを設定
clevels = np.arange(-20,80,10) # 塗りつぶし用
cnt_levels = np.arange(-20,80,10) #コンター用

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

# ベタ塗り
uwnd_mean_xt.plot.contourf(ax=ax,cmap='rainbow',levels=clevels)


# コンター
ctr = uwnd_mean_xt.plot.contour(ax=ax,levels=cnt_levels,colors='black')
ctr.clabel(fmt='%3i',fontsize=12,colors='black')

# Y軸を上下反転させる(気圧は大きいほど)
ax.invert_yaxis()
# 気圧座標をlogスケールに
ax.set_yscale('log')

ax.tick_params(labelsize=15)

# 軸のタイトル
ax.set_ylabel('Pressure (hPa)',size=15)
ax.set_xlabel('Latitude ($^\circ$N)',size=15) # $...$はTexでの数式の表現

plt.show()

このような図が出れば正解！（対流圏のジェット気流がきれいに見えますね）

課題¶

【課題 3-1】6時間間隔のジオポテンシャルハイトデータ（2019年01月、北緯40度）を用いてホフメラー図を描いてみよ（＋4 点）。高度によって(例：500 hPaと10 hPa）でどのように特徴が変わるだろうか。¶

（Tips) 東西偏差（東西平均からの偏差）を描画すると見やすい。

【課題 3-2】熱帯上空で経度-高度（気圧）断面図を作成し、ウォーカー循環の構造を確認せよ。東西風と鉛直風(鉛直p-速度）を重ねて描いてみると良い（どちらかがコンター、どちらかがカラー）。（注）鉛直風は12層しかデータがないので注意（＋5 点）¶

（ヒント）　熱帯は例えば、10S-10Nの緯度平均とする。

【課題 3-3】30 hPa面（～高度25 km）の気温（東西平均）データを用いて時間-緯度断面図を作成し、極域の気温の季節変化やその南北半球の対称性について考察せよ。但し、1980年以前（衛生観測以前）のデータは（南半球は特に）信頼性が低いので最近のデータを中心に眺めると良い（＋3 点）¶

（7/14 文言修正しました:「低温領域（極渦）の季節変化」->「極域の気温の季節変化」）

【課題 3-4】赤道上空東西風を用いて時間-高度（気圧）断面図を作成し、東西風準二年周期変動が見られることを確かめよ。（＋2 点）¶

2.5【発展】ベクトル場¶

ベクトル場をプロットするメソッドは無いので、matplotlibを使う。

fileu = "../../DB_2020/data/ncep.reanalysis.derived/pressure/uwnd.mon.mean.nc" #要変更
filev = "../../DB_2020/data/ncep.reanalysis.derived/pressure/vwnd.mon.mean.nc" #要変更

# ファイルの読み込み
U = xr.open_dataset(fileu)
V = xr.open_dataset(filev)

p_cut = 500; time_cut = '2017-01'
Usub = U['uwnd'].sel(level=p_cut,time=time_cut).squeeze()
Vsub = V['vwnd'].sel(level=p_cut,time=time_cut).squeeze()

# low-resolution data
# (経度-緯度方向は適当に間引いて読み込む; そうしないと、描いたベクトルがギチギチになるので）
dx = 4; dy = 2
U_lr = Usub[::dy,::dx].values #[::c]は "c毎にデータを取得する"]の意味, ".values"はndarrayへの変換（軸情報を落として単なる数値配列に）
V_lr = Vsub[::dy,::dx].values

# 軸もそれに合わせて間引く
lon_lr = Usub['lon'][::dx].values; lat_lr = Vsub['lat'][::dy].values


# ここから描画
fig = plt.figure(figsize=(10,7))
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
#ax = plt.axes(projection=ccrs.NearsidePerspective(central_latitude=30,central_longitude=140))

# ベクトル場を描く
Q = ax.quiver(lon_lr,lat_lr,U_lr,V_lr,color='blue',transform=ccrs.PlateCarree())

# 単位ベクトルを描く
ax.quiverkey(Q,0.9,0.8,20,label='20 m s$^{-1}$',labelpos='E',coordinates="figure")

# 目盛り
xticks = np.arange(13)*30
yticks = -90 + np.arange(7)*30
ax.set_xticks(xticks,crs=ccrs.PlateCarree())
ax.set_yticks(yticks,crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)


ax.coastlines()
#ax.stock_img() # 地形を描く
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)

#plt.savefig('test_vector.png') #画面出力だとユニットベクトルが見えないが、保存データにはちゃんと描かれている

plt.show()

C:\Users\takas\Anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:1752: UserWarning: Some vectors at source domain corners may not have been transformed correctly
  u, v = self.projection.transform_vectors(t, x, y, u, v)

3. 全球大気再解析データ（４次元データ）の解析¶

3.0 はじめに¶

3.1 データの読み込み：NetCDFデータを読み込む¶

3.2 データを切り出す¶

3.2 統計解析・時系列解析¶

3.4 絵を描く¶

$\ \ \ \ $ 3.4.1 二次元プロットの方法¶

$\ \ \ \ $ 3.4.2 地図描画¶

3.5 【参考】ベクトル場プロット¶

3.0 はじめに¶

XArrayについて¶

...さらに凝った解析や描画を行うには...¶

【練習課題】 2010-2018年平均の日本上空(500hPa)の東西風の分布（緯度-経度断面）を図示する¶

3.1 データの読み込み：NetCDFデータを読み込む¶

3.2 データを切り出す¶

3.2.1 ラベルで範囲指定¶

(1) 軸の名前と値を指定した切り出し: .sel(), .isel() （お勧め！）¶

(2) 各軸のラベル・整数のインデックス値を用いた切り出し：.loc[ ]、[ ]¶

3.2.2 条件指定して切り出し¶

◆ 時刻データを元に条件指定する際の注意点¶

【補足】 .squeeze() メソッド: 大きさ１の次元を削除¶

3.3 統計解析・時系列解析¶

◆ 統計解析¶

◆ 時系列データ解析¶

・resampleの例（年平均データの作成）：¶

・groupbyの例（各月のコンポジット）：¶

◆ さらに凝った解析を行うには？？¶

3.4 データの二次元描画¶

3.4.1 二次元描画：.plot.contourf（塗りつぶし）, .plot.contour（コンター）¶

【コンターラベルの作成】¶

3.4.2 地図描画¶

◆ 作成の手順（通常の描画と違う点を中心に）¶

別の地図投影も試してみよう（宇宙から見た図）¶

【練習課題（1）】2010-2018年データを用いて、時間・東西平均した東西風を緯度-高度（気圧）断面に描いてみよ（一行目の処理を埋めよ）¶

【練習課題（２）】各月（例えば、１月と７月）の平均値を用いて（１）と同様の図を描いてみよ。（「月」の情報の取り出し方は、3.2.2を参照）¶

課題¶

【課題 3-1】6時間間隔のジオポテンシャルハイトデータ（2019年01月、北緯40度）を用いてホフメラー図を描いてみよ（＋4 点）。高度によって(例：500 hPaと10 hPa）でどのように特徴が変わるだろうか。¶

【課題 3-4】赤道上空東西風を用いて時間-高度（気圧）断面図を作成し、東西風準二年周期変動が見られることを確かめよ。（＋2 点）¶

2.5【発展】ベクトル場¶

3.　全球大気再解析データ（４次元データ）の解析¶

(1) 軸の名前と値を指定した切り出し: .sel(), .isel()　（お勧め！）¶