これまでの気候の変化（観測結果）を解析する

はじめに

先日、9th GEWEX-OSC 市民講座　「みんなで考えるこれからの気候変動」にオンライン参加しました。

sites.google.com

札幌管区気象台からの話題提供で、北海道全体での気温の上昇傾向について紹介されていたのですが、その中で(恣意的に分割した場合という前提で）温暖化のスピードが１９９０年代に加速しているようにみえるのでは？とのこと。

www.data.jma.go.jp

ここにデータが公開されているので入手してグラフ化してみると、

%cd cd /mnt/c/Users/hoge/AMeDAS/SapporoKanku/

import numpy as np
import pandas as pd
import plotly.express as px
from sklearn import linear_model

hokkaido_tmean_df = pd.read_csv("hokkaido_tmean_ann_obs.csv")
hokkaido_tmean_df.columns = ["year","tmean"]

fig = px.scatter(hokkaido_tmean_df,x="year",y="tmean",width=800,height=500,trendline='ols',trendline_color_override='orange',)
fig.update_traces(selector=0,marker_size=8,mode="lines+markers")
fig.update_traces(selector=1,line_width=6)
fig.update_layout(margin=dict(l=60, r=20, t=20, b=20))
fig.update_xaxes(title="年")
fig.update_yaxes(title="1991-2020年平均からの差(℃)")

たしかに1990年代以前と以降で別々の線を書きたくなります。

本当に1990年代にトレンドが変化したのか？長期的なトレンド＋年々変動のばらつきの範囲内なのか？変化点解析を行ってみました。

データの入手と整形

札幌管区気象台のホームページで公開されている「北海道地方のこれまでの気候の変化（観測結果）」から、長期間観測を継続している6地点（旭川、網走、札幌、帯広、根室、函館）の年平均気温データを入手しました。

%cd cd /mnt/c/Users/hoge/AMeDAS/SapporoKanku/

import numpy as np
import pandas as pd
import plotly.express as px

data_df = pd.DataFrame(index=[])

for h_name in ["abashiri","asahikawa","hakodate","nemuro","obihiro","sapporo" ]:
    tdf = pd.read_csv(f"{h_name}_tmean_ann_obs.csv")
    tdf.columns = ["year","tmean"]
    tdf["sta_name"] = h_name
    data_df = pd.concat([data_df,tdf])

data_df = data_df.set_index(["sta_name","year"]).unstack(level=0)
data_df = data_df.dropna()
data_df.columns = data_df.columns.get_level_values(1)

fig = px.line(data_df,width=800,height=500)
fig.update_xaxes(title="年")
fig.update_yaxes(title="年平均気温(℃)")

1991-2020年平均からの差を計算する
data_std_df = data_df - data_df.loc[1991:2020].mean()

fig = px.line(data_std_df,width=800,height=500)
fig.update_xaxes(title="年")
fig.update_yaxes(title="1991-2020年平均からの差(℃)")

モデル構築

前提として、次の条件を設定しました。

• 変曲点は１つ
• 変曲点の前後は直線のトレンド
• すべての地点でトレンドの変化は同時に起こっている
• 年々変動のばらつきはすべての地点で同じ

変化点の前と後でトレンドが変化するモデルを以下の式で定義しました

ここで、はトレンドが変化する年、は年々変動のばらつきで、年に傾きがからに変化するようにしています。

pymcでモデルを組み立て、サンプリングします。

import pymc
import

with pymc.Model() as model:
    #次元・インデックスを定義
    model.add_coord("year",data_std_df.index,mutable=True)
    model.add_coord("sta",data_std_df.columns,mutable=True)
    model.add_coord("period",[0,1],mutable=True)
    n_year = data_std_df.shape[0]
    n_sta = data_std_df.shape[1]
    #データ読み込み
    X=pymc.Data("X",np.arange(n_year),dims=("year",),mutable=True)
    y=pymc.Data("y",data_std_df,dims=("year","sta"),mutable=True)
    #事前分布
    sigma = pymc.HalfCauchy("sigma", beta=10)
    intercept = pymc.Normal("intercept", 0, sigma=20,dims=("sta",))
    slope = pymc.Normal("slope", 0, sigma=20,dims=("sta","period"))
    cp = pymc.DiscreteUniform('cp', lower=0, upper=n_year-2)
    #muの計算
    X_bcast = pymc.math.broadcast_to(X,[n_sta,n_year])
    mu1 = pymc.Deterministic("mu1",intercept[:]+ slope[:,0]X_bcast.T,dims=("year","sta",))
    mu2 = pymc.Deterministic("mu2",(intercept[:]+(slope[:,0] - slope[:,1])cp) + slope[:,1]*X_bcast.T,dims=("year","sta",))
    #変化点での切り替え
    idx = pymc.math.broadcast_to(np.arange(n_year),[n_sta,n_year]).T
    mu = pymc.math.switch(idx <= cp, mu1, mu2)
    #尤度
    likelihood = pymc.Normal("lklh", mu=mu, sigma=sigma, observed=y)
    #サンプリング
    idata = pymc.sample(3000,idata_kwargs={'log_likelihood':True})

模式図を表示

pymc.model_to_graphviz(model)

サンプリング結果をプロット

az.plot_trace(idata)

うまく混交＆収束していそうです。

結果１　トレンドが変化した年

cpのサンプリング結果から、いつトレンドが変化していそうか？を見ることができます。

import plotly.express as px
fig = px.histogram(x=data_std_df.index[idata.posterior["cp"].data.flatten()],width=500,height=400)
fig.update_xaxes(title="年")
fig.update_layout(margin=dict(l=40, r=20, t=20, b=20))

1975~1991年まで分布しており、グラフの見た目と一致しています。そのうち1985年が最も多く、この年にトレンドが変化していると考えて良さそうです。

結果2　傾きの変化

slopeのサンプリング結果から、それぞれの地点で傾きがどの程度変化しているかを見ることができます

from plotly.subplots import make_subplots
import plotly.graph_objects as go
fig = make_subplots(rows=6, cols=1,subplot_titles=("旭川", "札幌", "網走","根室","帯広","函館"),shared_xaxes=True,vertical_spacing=0.02,x_title="傾き")

fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="asahikawa",period=0).data.flatten(),name="変曲前",marker=go.histogram.Marker(color="lightblue"),legendgroup = '1'),row=1, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="asahikawa",period=1).data.flatten(),name="変曲後",marker=go.histogram.Marker(color="lightsalmon"),legendgroup = '1'),row=1, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="sapporo",period=0).data.flatten(),name="変曲前",marker=go.histogram.Marker(color="lightblue"),showlegend=False),row=2, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="sapporo",period=1).data.flatten(),name="変曲後",marker=go.histogram.Marker(color="lightsalmon"),showlegend=False),row=2, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="abashiri",period=0).data.flatten(),name="変曲前",marker=go.histogram.Marker(color="lightblue"),showlegend=False),row=3, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="abashiri",period=1).data.flatten(),name="変曲後",marker=go.histogram.Marker(color="lightsalmon"),showlegend=False),row=3, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="nemuro",period=0).data.flatten(),name="変曲前",marker=go.histogram.Marker(color="lightblue"),showlegend=False),row=4, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="nemuro",period=1).data.flatten(),name="変曲後",marker=go.histogram.Marker(color="lightsalmon"),showlegend=False),row=4, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="obihiro",period=0).data.flatten(),name="変曲前",marker=go.histogram.Marker(color="lightblue"),showlegend=False),row=5, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="obihiro",period=1).data.flatten(),name="変曲後",marker=go.histogram.Marker(color="lightsalmon"),showlegend=False),row=5, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="hakodate",period=0).data.flatten(),name="変曲前",marker=go.histogram.Marker(color="lightblue"),showlegend=False),row=6, col=1)
fig.add_trace(go.Histogram(x=idata.posterior["slope"].sel(sta="hakodate",period=1).data.flatten(),name="変曲後",marker=go.histogram.Marker(color="lightsalmon"),showlegend=False),row=6, col=1)

for anot in fig.layout.annotations[:-1]:
    anot.update(x=0.06,yshift=-25)
fig.update_layout(barmode='overlay',height=800,width=500,legend_x=0.8,legend_y=0.99,margin=dict(l=30, r=10, t=20, b=50))

旭川や札幌で変化量が小さく、道東、道南では大きく変化しているようです。

結果3 :北海道地方の平均気温トレンドの解析

先程の地点毎の解析で得られた年々変動の分散と変化点を使って、北海道地方の年平均気温のトレンドを解析してみます。

with pymc.Model() as model:
    # coords
    model.add_coord("year",hokkaido_df["year"],mutable=True)
    model.add_coord("period",[0,1],mutable=True)
    n_year = hokkaido_df.shape[0]
    # data
    X=pymc.Data("X",np.arange(n_year),dims=("year",),mutable=True)
    y=pymc.Data("y",hokkaido_df["tmean"],dims=("year",),mutable=True)
    # Set parameter
    sigma = 0.529
    cp = np.argwhere(hokkaido_df["year"] == data_std_df.index[92])[0][0]
    # Define priors
    intercept = pymc.Normal("intercept", 0, sigma=20)
    slope = pymc.Normal("slope", 0, sigma=20,dims=("period"))
    # mu1,mu2 →mu
    mu1 = pymc.Deterministic("mu1",intercept+ slope[0]X,dims=("year",))
    mu2 = pymc.Deterministic("mu2",(intercept+(slope[0] - slope[1])cp) + slope[1]*X,dims=("year",))
    idx = np.arange(n_year)
    mu = pymc.math.switch(idx <= cp, mu1, mu2)
    # likelihood
    likelihood = pymc.Normal("lklh", mu=mu, sigma=sigma, observed=y)
    # sampling
    idata = pymc.sample(3000,idata_kwargs={'log_likelihood':True})

トレンドのグラフを作成

summary_df = az.summary(idata)

cp_x = data_std_df.index[92]
cp_y = summary_df.loc["slope[0]","mean"]92+summary_df.loc["intercept","mean"]
x_0 = hokkaido_df["year"].min()
y_0 = summary_df.loc["intercept","mean"]
x_1 = hokkaido_df["year"].max()
y_1 = summary_df.loc["slope[1]","mean"](x_1 - cp_x) + cp_y

fig = px.scatter(hokkaido_df,x="year",y="tmean",width=800,height=500)
fig.update_traces(selector=0,marker_size=8,mode="lines+markers")
fig.update_layout(margin=dict(l=60, r=20, t=20, b=20))
fig.update_xaxes(title="年")
fig.update_yaxes(title="1991-2020年平均からの差(℃)")
fig.add_shape(type="line",x0=x_0, y0=y_0, x1=cp_x, y1=cp_y,line=dict(width=6,color='orange'))
fig.add_shape(type="line",x0=cp_x, y0=cp_y, x1=x_1, y1=y_1,line=dict(width=6,color='orange'))

1984年にトレンドが変化するグラフになりました。

終わりに

観測データのトレンドが1984年に変化しているからといって、その年になにかイベントがあったということではありません。（相関と因果を混同してはいけない）

ただ近年、温暖化が加速していることについてはいろいろ議論されているようです。

www.carbonbrief.org

一昔前の「ハイエイタス」の議論の逆パターンで、今回解析した短期の傾向がそのまま「人為的地球温暖化」が加速している証拠にはならないことに注意が必要です。

ただ、地球沸騰化の実感とは一致しているかも。

ccilabo_director 2024-07-11 18:35 読者になる

広告を非表示にする

関連記事

• 

  2024-09-20

  気候変動が進んだ将来、「かんしょ」の生産はどのように変化するか？ ―その２ 10a当たり収量の変化予…

  はじめに 前回は、「かんしょ」の10a当たり収量について平均気…
• 

  2024-09-12

  気候変動が進んだ将来、「かんしょ」の生産はどのように変化するか？ ―その１ 10a当たり収量予測モデ…

  はじめに 私事ですが、先月無事に某育成講座を終了し、晴れて「…
• 

  2023-10-12

  気候予測データの解析環境を構築する。その7―MCMCサンプラーの実行環境を構築する。

  はじめに 気候予測データは、将来の平均的な「気候」状態を予測…
• 

  2023-09-29

  「全球及び日本域150年連続実験データ」を可視化する その7－札幌の平均気温のシミュレーション結果を…

  はじめに 今年の夏は北海道でもうだるような暑さが続きました。…
• 

  2023-09-15

  「全球及び日本域150年連続実験データ」を可視化する その6－回る地球上に地上気温変化を描く。

  はじめに 前回は、QGIS上で全球の地上気温変化アニメーションを…

• もっと読む

コメントを書く