アメダス気象データ分析チャレンジ！（Python版）動作確認プログラム

Copyright 2023 気象ビジネス推進コンソーシアム、岐阜大学　吉野純
(C) 2023 WXBC、岐阜大学　吉野純
＜利用条件＞
本書は、本書に記載した要件・技術・方式に関する内容が変更されないこと、および出典 を明示いただくことを前提に、無償でその全部または⼀部を複製、翻案、翻訳、転記、引 用、公衆送信等して利用できます。なお、全体を複製、翻案、翻訳された場合は、本書にあ る著作権表示および利用条件を明示してください。
＜免責事項＞
本書の著作権者は、本書の記載内容に関して、その正確性、商品性、利用目的への適合性 等に関して保証するものではなく、特許権、著作権、その他の権利を侵害していないことを 保証するものでもありません。本書の利用により生じた損害について、本書の著作権者は、 法律上のいかなる責任も負いません。

Anaconda3のインストールが完了したら，Pythonのサンプルプログラム（python_test.ipynb）が動くかどうか確かめてみましょう．

プログラムは以下の手順で実行します．

1. 下のセル（プログラムが書かれている場所）をマウスで選択します．
   
2. 上のツールバーの「▶実行」をクリックします．
   
3. 下の方に５つの図と１つの表が出ることを確認します．
   

上記の1～3が上手くいかない場合（図が出ずにエラーメッセージのみが出る場合），Anaconda3のインストールに失敗している可能性があります．再度，以下のホームページから

https://www.anaconda.com/

最新版（2023.07）の Anaconda3-2023.07 をインストールしてください．それでもうまく行かない場合には，以下のお問い合わせ先にご連絡ください．

上記の1～3が上手くいった場合（エラーメッセージは出ずに図が表示される場合），Anaconda3の環境構築はこれで終了です．お疲れさまでした．アメダス気象データ分析チャレンジ！（Python版）の当日にお目にかかれますことを楽しみにしています．

応募者多数により惜しくも抽選で外れてしまった場合には，大変申し訳ございません．当日の資料は後日ホームページ上で公開しますので，この環境を活用して是非，アメダス気象データ分析に挑戦してみてください．

# numpyをnpという別名でインポートします．
import numpy as np
# matplotlibをpltという別名でインポートします．
import matplotlib.pyplot as plt
# Seabornをインポートします．
import seaborn as sns
# pandasをpdという別名でインポートします．
import pandas as pd
# sklearn(scikit-learn)は機械学習関連のライブラリーです．インポートします．
from sklearn import linear_model
%precision 3
%matplotlib inline
"""
4次のルンゲ-クッタ法によるローレンツモデル計算
dx/dt = -s*x + s*y
dy/dt = -y + r*x - x*z
dz/dt = -b*z + x*y
"""

# １つ目の方程式dx/dt = -s*x + s*yを関数として設定
def f1(t,x,y,z,s):
    return -s*x + s*y

# ２つ目の方程式dy/dt = -y + r*x - x*zを関数として設定
def f2(t,x,y,z,r):
    return -y + r*x - x*z

# ３つ目の方程式dz/dt = -b*z + x*yを関数として設定
def f3(t,x,y,z,b):
    return -b*z + x*y

# パラメータs,r,bの設定
s = 10.0
r = 28.0
b = 8.0/3.0

# データx,y,zの初期値の設定
x = 0.1
y = 0.1
z = 0.1

# 時間刻み
h = 0.01
# t=0からt=100までh刻みで計算します．時間tをnumpy配列で設定．
tvalues = np.arange(0.0,100.0,h)
# データxの保存先をnumpy配列とする
xvalues = []
# データyの保存先をnumpy配列とする
yvalues = []
# データzの保存先をnumpy配列とする
zvalues = []

# ルンゲ-クッタ法で時間積分
for t in tvalues:
    xvalues.append(x)
    yvalues.append(y)
    zvalues.append(z)
    k1 = h*f1(t,x,y,z,s)
    l1 = h*f2(t,x,y,z,r)
    m1 = h*f3(t,x,y,z,b)
    k2 = h*f1(t+0.5*h,x+0.5*k1,y+0.5*l1,z+0.5*m1,s)
    l2 = h*f2(t+0.5*h,x+0.5*k1,y+0.5*l1,z+0.5*m1,r)
    m2 = h*f3(t+0.5*h,x+0.5*k1,y+0.5*l1,z+0.5*m1,b)
    k3 = h*f1(t+0.5*h,x+0.5*k2,y+0.5*l2,z+0.5*m2,s)
    l3 = h*f2(t+0.5*h,x+0.5*k2,y+0.5*l2,z+0.5*m2,r)
    m3 = h*f3(t+0.5*h,x+0.5*k2,y+0.5*l2,z+0.5*m2,b)
    k4 = h*f1(t+h,x+k3,y+l3,z+m3,s)
    l4 = h*f2(t+h,x+k3,y+l3,z+m3,r)
    m4 = h*f3(t+h,x+k3,y+l3,z+m3,b)
    x += (k1+2*k2+2*k3+k4)/6
    y += (l1+2*l2+2*l3+l4)/6
    z += (m1+2*m2+2*m3+m4)/6

# 数値積分おわり
# 後処理
# xの時系列図
print("xの時系列図です")
plt.plot (tvalues, xvalues)
plt.xlabel("t")
plt.ylabel("x(t)")
plt.show()

# yの時系列図
print("yの時系列図です")
plt.plot (tvalues, yvalues)
plt.xlabel("t")
plt.ylabel("y(t)")
plt.show()

# zの時系列図
print("zの時系列図です")
plt.plot (tvalues, zvalues)
plt.xlabel("t")
plt.ylabel("z(t)")
plt.show()

# x-z断面図
print("x-zの断面図です")
plt.plot (xvalues, zvalues)
plt.xlabel("x(t)")
plt.ylabel("z(t)")
plt.show()

# Pandasのデータフレーム形式に変換
print("データフレームに割り当てて，相関係数を計算します")
df = pd.DataFrame({'t':tvalues,'x':xvalues,'y':yvalues,'z':zvalues})
corr=df['x'].corr(df['y'])

# 線形回帰直線
print("xとyの散布図と回帰直線を引きます")
clf = linear_model.LinearRegression()
X=df.loc[:,['x']].values
Y=df.loc[:,['y']].values
clf.fit(X,Y)
plt.plot(X, Y, 'o')
plt.plot(X, clf.predict(X))
plt.xlabel("x(t)")
plt.ylabel("y(t)")
plt.show()

print('xとyの相関係数は',corr)
df

xの時系列図です

yの時系列図です

zの時系列図です

x-zの断面図です

データフレームに割り当てて，相関係数を計算します
xとyの散布図と回帰直線を引きます

xとyの相関係数は 0.8867624304628634

	t	x	y	z
0	0.00	0.100000	0.100000	0.100000
1	0.01	0.101300	0.126889	0.097481
2	0.02	0.105051	0.154219	0.095059
3	0.03	0.111094	0.182642	0.092737
4	0.04	0.119344	0.212771	0.090521
...	...	...	...	...
9995	99.95	4.805565	2.329803	26.451100
9996	99.96	4.572819	2.392494	25.864153
9997	99.97	4.369234	2.476375	25.290816
9998	99.98	4.193928	2.578975	24.731966
9999	99.99	4.045893	2.698347	24.188313

10000 rows × 4 columns

ちなみに

上のサンプルのプログラムは，気象学における熱対流（積乱雲も対流の１つ）の問題を表す微分方程式であり，「ローレンツモデル」と呼ばれています（Lorenz, 1963）．方程式は以下の通りです．

$\frac{dx}{dt} = - \sigma x + \sigma y $

$\frac{dy}{dt} = - y + r x - xz $

$\frac{dz}{dt} = - b z + x y $

$x$は水平方向の風速の変動，$y$は水平方向の温度の変動，$z$は鉛直方向の温度の変動（のようなもの）を表しています．ルンゲクッタ法という方法で計算しています．余裕がある人はプログラムを解読してみてください．

$x$と$y$と$z$が連動しながら不規則に変化するカオス(chaos) の振る舞いを示します．$x$と$z$の断面図はチョウのような模様（ローレンツバタフライといいます）になっています．初期値（33～35行目）を変えると結果も大きく変わります．このように大気はカオス的性質を持つために，天気予報で用いられる数値予報モデルでは初期値を少しずつ変えたアンサンブル予報が行われています．

話がそれてしまいました．．．

難しそうな式が出てきてびっくりしたかもしれませんが，アメダス気象データ分析チャレンジ！を受講するにあたって上記を理解している必要は全くありません のでご安心ください．Pythonのプログラムのサンプルとして用意しました．