In [1]:
import numpy as np
import matplotlib.pyplot as plt
import netCDF4
import cartopy.crs as ccrs
import cartopy.util as cutil
import datetime
import pandas as pd
import copy
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
## 図を簡単に描くために関数として定義しておく
def fig_y_z(lat,z,data,cmap='rainbow',min=-20,max=60,int=10):
# min = -20 などとなっているのは、関数を呼び出す際にに指定しなければこの値を使うという意味（デフォルト値）
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
clevels = np.arange(min,max+int,int)
cnt_levels = np.arange(min,max+int,int)
cf = ax.contourf(lat,z,data,levels=clevels,cmap=cmap)
# Y軸を反転させる
ax.invert_yaxis()
# 気圧座標をlogスケールに
ax.set_yscale('log')
# 軸のタイトル
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Latitude ($^\circ$N)')
###カラーバー
cbar = plt.colorbar(cf,orientation='vertical')
#ax.contour(lon,lat,zm,cnt_levels,colors = 'white')
cbar.set_label('Zonal wind [m/s]')
plt.show()
def fig_x_y(lon,lat,data,cmap='rainbow',min=-20,max=60,int=10,title='dummy',var_name ='var_name',var_unit='var_unit'):
# min = -20 などとなっているのは、関数を呼び出す際にに指定しなければこの値を使うという意味（デフォルト値）
# 0..357.5（経度方向）のデータを、0..360の範囲に拡張しておく。
data_cyclic, lon_cyclic = cutil.add_cyclic_point(data,coord=lon)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=180))
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.gridlines(draw_labels=False) #ラベルはすでに上で描いたので"False"
##ax.set_extent([-180,180,-90,90],crs=ccrs.PlateCarree(central_longitude=180))
clevels = np.arange(min,max+int,int)
cnt_levels = np.arange(min,max+int,int)
cf = ax.contourf(lon_cyclic,lat,data_cyclic,transform=ccrs.PlateCarree(),levels=clevels,cmap=cmap)
# transform = にはデータ自身の座標系を指定する
ax.coastlines()
cbar = plt.colorbar(cf,orientation='horizontal')
cbar.set_label(var_name + var_unit)
plt.show()
In [2]:
file = "../../DD/Fortran/data/uwnd.mon.mean.nc"
# ファイルの読み込み
data = netCDF4.Dataset(file,'r')
# 軸・変数の取得
lon = data.variables['lon'][:] # "[:]" は全要素の値を取得するという意味, ".getValue"としても同じ
print(type(lon))
lat = data.variables['lat'][:]
level = data.variables['level'][:]
uwnd = data.variables['uwnd'][:,:,:,:]
var_name = data.variables['uwnd'].long_name
var_unit = data.variables['uwnd'].units
print(var_name,var_unit)
# 日付軸
time = data.variables['time']
time_date = netCDF4.num2date(time[:],units=time.units)
time_date_pd = pd.to_datetime(time_date)
In [3]:
zm = uwnd.mean(axis=(0,-1))
fig_y_z(lat,level,zm,min=-20,max=50,int=10)
（例題 2）時間➖東西平均からの偏差場をx-y平面に描け。¶
形状の異なる配列同士の計算：ブロードキャスト¶
元の時系列と時間➖東西平均場、および、偏差場は：
\begin{align} U(x,y,z,t) = [\overline U(y,z)] + U'(x,y,z,t) \end{align}という関係になっている。ここで、
- 元の時系列($U$): 大きさ(843,17,73,144)の配列
- 時間➖東西平均($[\overline U]$)： 大きさ(17,73)の配列
- 偏差場($U'$): 大きさ(843,17,73,144)の配列
である。
偏差場を求める際には$U$から$[\overline U]$を差し引く操作が必要であるが、
形状の異なる２つの配列間での演算が必要にとり、注意が必要である。
NumPyではこういった状況においても、以下の条件を満たしている場合には
"うまいこと"演算してくれる（ブロードキャストという）。
ブロードキャストの規則（参考文献[]のp.500より引用; 一部改変）¶
２つの配列に対してブロードキャストが適用できる条件は、両配列の最も高い次元から低い方向へ、対応する次元をそれぞれ調べていった 場合に、
(a) 両者の配列の長さが一致しているか
(b) 一方の配列の長さが１であることです。この場合、ブロードキャストは、
(I) 一方に存在しない次元や
(II) 長さが1の次元に適用されます。
これだけでは良くわからないので、実例を用いて説明する。まずは単純に$U-[\overline U]$として計算してみよう。
In [4]:
anom = uwnd - zm
、、とエラーが出てしまった。「ブロードキャストがうまくできない」と怒られている。
上の規則と照らし合わせて、uwnd(元データ: $U$)とzm(時間-東西平均: $[\overline U]$)の次元の大きさを比較してみると:
| 次元 | 3 | 2 | 1 | 0 |
|---|---|---|---|---|
| uwnd | 144 | 73 | 17 | 843 |
| zm | 73 | 17 | - | - |
となっていて、最も高い次元(3)のところで上記規則の(a)(b)のいずれも満たしていないことが分かる。
これを解決するには、以下のように新たにzmに大きさ１の次元を追加すれば良い（np.newaxisを使う；下記参照）：
| 次元 | 3 | 2 | 1 | 0 |
|---|---|---|---|---|
| uwnd | 144 | 73 | 17 | 843 |
| zm | 1 | 73 | 17 | - |
実際このときには、高い次元から順に比べていくと
- 次元3: 規則(b)を満たす -> 適用規則(II)によりブロードキャスト
- 次元2: 規則(a)を満たす
- 次元1: 規則(a)を満たす
- 次元0: zmには存在しないが -> 適用規則(I)によりブロードキャスト
となってブロードキャストが適用され、うまく演算を行える。
大きさ１の新たな軸を追加する： np.newaxis¶
np.newaxisという特殊な属性を用いると、任意の場所に大きさ１の軸を挿入することができる（次元は増えるが、大きさ１なので実態はなにも変わらない）。
In [5]:
# 新しい次元を挿入
zm1 = zm[:,:,np.newaxis]
# 確認
print(zm1.shape)
# ブロードキャストで計算できる
anom = uwnd - zm1
iz_mask = (level == 500)
it_mask = (time_date_pd == datetime.datetime(2017,1,1))
fig_x_y(lon,lat,anom[it_mask,iz_mask,:,:].squeeze(),min=-30,max=30,int=5)
（課題）：12-2月、6-8月それぞれの東西平均と標準偏差の場を図示せよ（配布資料参考）。¶
In [6]:
mask = (time_date_pd.month == 6) | (time_date_pd.month == 7) | (time_date_pd.month == 8)
zm = uwnd[mask].mean(axis=(0,-1))
#
fig_y_z(lat,level,zm,min=-30,max=80,int=10)
In [7]:
zm = uwnd[mask].mean(axis=-1)
std = zm.std(axis=0)
fig_y_z(lat,level,std,min=0,max=20,int=2)
In [ ]: