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

はじめに

今年の夏は北海道でもうだるような暑さが続きました。気象庁の”夏の天候のまとめ”によれば、2023年（令和5年）夏（6～8月）の日本の平均気温は1898年以降で夏として最も高くなったとのこと。

www.jma.go.jp

別紙（概況、統計値等）によれば、北海道の平均気温は平年より3℃も高かったらしい。

https://www.data.jma.go.jp/obd/stats/data/stat/tenko2023jja_besshi.pdf

実際、これまでに経験したことのない暑さだったことは間違いないのですが、気温の変化をよりわかりすく示すヒートマップというものがあるようです

www.sakigake.jp

www.nishinippon.co.jp

これはわかりやすい！ということで、１５０連続実験データでもこの表現手法を使ってヒートマップを作成してみたいと思います。

データの切り出し

ここでは１例として札幌管区気象台に最も近い格子点における８月の平均気温のヒートマップを作成します。

地理座標は北緯43度3.6分、東経141度19.7分ですので、対応する格子点の位置を特定します。

https://www.jma.go.jp/jma/kishou/know/amedas/ame_master.pdf

# ライブラリをインポート
import os
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon

os.chdir("/mnt/c/Users/hoge/ONEFIFTY")

# 札幌管区気象台の緯度経度
src_lat = 43 + 3.6/60
src_lon = 141 + 19.7/60

# 日本域20㎞150年連続実験データの投影法での座標を計算
wgs_SRS = osr.SpatialReference()
wgs_SRS.ImportFromEPSG(4326)

lcc_SRS = osr.SpatialReference()
lcc_SRS.ImportFromProj4(
"+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=135. +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")

trans = osr.CoordinateTransformation(wgs_SRS,lcc_SRS)
out_x,out_y,z = trans.TransformPoint(src_lat,src_lon)

#最も近い格子点の行番号、列番号を特定する
x_vec = (-1930000.1287495417 + 10000) + np.arange(0,191) * 20000
y_vec = ((5833346.794147097 - 10000) - np.arange(0,155) * 20000)[::-1]

x_ind = np.argmin(np.abs(x_vec - out_x))
y_ind = np.argmin(np.abs(y_vec - out_y))

# 確認のためのベクタを作成する
x = x_vec[x_ind]
y = y_vec[y_ind]

df = pd.DataFrame(["Sapporo"])
df.columns=["name"]
gm = [Polygon([(x-10000,y-10000),(x-10000,y+10000),(x+10000,y+10000),(x+10000,y-10000),(x-10000,y-10000)])]

sapporo_gdf = gpd.GeoDataFrame(df,geometry=gm,crs=lcc_SRS.ExportToWkt())
sapporo_gdf.to_crs("EPSG:4326").to_file("sappro_grid.fgb", driver="FlatGeobuf")
print(x_ind,y_ind) # 121 121

格子点を中心とした20㎞メッシュに札幌が含まれていることが確認できます。

150年連続実験データの処理

「全球及び日本域150年連続実験データ」の20km150年連続実験データから、（121,121）に位置するデータを切り出し、地上気温の日平均値を計算します。

# ライブラリをインポート
import os
import pathlib
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray

# 処理するファイルの一覧を作成する
os.chdir("/mnt/c/Users/hoge/ONEFIFTY")

griblist_df = pd.DataFrame({"grib_path":list(pathlib.Path("./").glob("**/*.grib"))})
griblist_df["grib_name"] = griblist_df["grib_path"].map(lambda x: x.name)
griblist_df["year"] = griblist_df["grib_name"].str.extract(".*_(\d{4})\d{2}.grib").astype(int)
griblist_df["month"] = griblist_df["grib_name"].str.extract(".*_\d{4}(\d{2}).grib").astype(int)
griblist_df = griblist_df.sort_values("year").reset_index(drop=True)

# 年ごとの日平均値を収納する配列を作成
TA_daily_mean_arr = np.zeros((griblist_df.shape[0],31))

# 年ごとに日平均値を計算
for h_index,h_row in griblist_df.iterrows():
    ds = xr.open_dataset(h_row["grib_path"], engine="cfgrib")
    ds["time"] = ds["time"] - np.timedelta64(1, 'h')
    da = ds["t"][:,121,121].resample(time='1D').mean()
    TA_daily_mean_arr[h_index,:] = da.data

np.save("TA_daily_mean.npy",TA_daily_mean_arr)

ヒートマップの作成

plotlyでヒートマップを描きます。

# ライブラリをインポート
import os
import numpy as np
import plotly.express as px
import matplotlib.cm as cm

os.chdir("/mnt/c/Users/hoge/ONEFIFTY")

#　データをロード
TA_daily_mean_arr = np.load("TA_daily_mean.npy")

#　平年値（統計期間1991～2020年）を計算
TA_normal = TA_daily_mean_arr[40:70,:].mean()

# ヒートマップの描画
fig = px.imshow(TA_daily_mean_arr, aspect="auto",width=550,height=800, 
color_continuous_scale='RdBu_r',
color_continuous_midpoint=TA_normal)

# レイアウト編集
fig.update_layout(
    title=dict(text="札幌周辺の８月の日平均気温シミュレーション（1951~2099）",
                y=0.98),
    margin=dict(t=40, b=45, l=25, r=120)
)   

#脚注１を作成
fig.add_annotation(
    font_size=10,
    font_family="sans-serif",
    x=1.04,
    y=0.83,
    xref="paper",
    yref="paper",
    xanchor="left",
    yanchor="top",
    borderwidth=1,
    bordercolor='black',
    align="left",
    text="「日本域150年連続<br>実験」による、札幌<br>に最も近い格子点に<br>
　　　　おける1951年から<br>2099年までの８月の<br>日平均気温シミュレ<br>
　　　　ーション結果。縦軸<br>が '年' , 横軸が '月日' 。<br>8月の日平均気温の<br>
　　　　平年値（1991~2020<br>年）を白とし、高温<br>を赤、低温を青とし<br>
　　　　た色分けにより気温<br>の違いを表現した。<br>結果は特定の日の「<br>
　　　　気象」の状態（何月<br>何日に気温は何度か）<br>を予測しているもの<br>
　　　　ではないことに注意！"
)

# 脚注２を作成
fig.add_annotation(
    font_size=10,
    font_family="sans-serif",
    x=1.03,
    y=0.12,
    xref="paper",
    yref="paper",
    xanchor="left",
    yanchor="top",
    align="left",
    showarrow=False,
    text="本ヒートマップは、<br>文部科学省「統合<br>的気候モデル高度<br>
　　　　化研究プログラム」<br>において、地球シ<br>ミュレータを用い<br>
　　　　て作成されたデー<br>タを使用させてい<br>ただいています。"
)

# y軸を設定
fig.update_yaxes(
    tickvals=[0,9,19,29,39,49,59,69,79,89,99,109,119,129,139,148],
    ticktext=["1951","1960","1970","1980","1990","2000","2010","2020","2030","2040","2050","2060","2070","2080","2090","2099"]
    )

# x軸を設定
fig.update_xaxes(
    tickvals=[0,4,9,14,19,24,29],
    ticktext=["08/01","08/05","08/10","08/15","08/20","08/25","08/30"],
    tickangle=-50
)

# カラーバーを設定
fig.update_coloraxes(
    colorbar_len=0.4,
    colorbar_y=0.32,
    colorbar_yanchor="middle",
    colorbar_ticklabeloverflow="allow",
    colorbar_tickvals=[284,310],
    colorbar_ticktext=["low","high"],
    colorbar_title="日平均気温",
    colorbar_title_side="top"
    )

# ロゴを描画
fig.add_layout_image(
    sizex=0.2,sizey=0.2,
    source="https://pbs.twimg.com/profile_images/1664190624701222912/8gaaJj-l_400x400.png",
    xref="paper",
    yref="paper",
    x=1.05,
    y=0.995)

fig.show()