「日本域気候予測データ‐観測地点データ」を可視化する その２－地点ごとの将来予測のグラフ化
はじめに
「日本域気候予測データ‐観測地点データ」の年統計値データには、それぞれ20 年分（現在：20世紀末（1980～1999年）、将来：21世紀末（2076～2095年）の年統計値データが収録されています。気候予測データセット2022解説書を読むと、
気候モデルによる予測結果には避けられない不確実性が伴います。そのため、その不確実性を理解し、適切に評価したうえで、予測データを気候変動対策に利用することが推奨されています。（気候予測データセット解説書第２章．P2-82より引用）
とあります、また、
気候モデルの予測結果に現れる年々～数十年周期の自然変動の位相等は予測の対象外であり、20 年間といった長い期間にわたる計算結果を解析することで、地球温暖化に伴う長期的な気候の変化を抽出することができる。（気候予測データセット解説書第２章．P2-82,８３より引用）
とされています。
つまり、「予測データを気候変動対策に利用する」ためには「20 年間といった長い期間にわたる計算結果を解析」し、比較する必要があります。
札幌管区気象台の公表している「北海道地方 地球温暖化予測情報」では、各振興局の代表地点の将来変化が、不確実性を表現した”箱ひげ図”で示されています。
前回ダウンロードした年統計値を用いて、他のアメダス地点の20年分データの統計的傾向の変化を比較し、同様の図を作成してみましょう。
夏日の年間日数のデータ集計
ここでは、夏日（日最高気温25℃以上）の日数データを使用します。要素名は"t2mx_ov25"です。
現在気候実験(SPA)のデータフレーム作成
現在気候再現実験（SPA）の年集計値は/mnt/c/Users/hoge/JWP9/dias/data/csv/SPA/annually/に解凍されており、ファイル名は次の規則でつけられています。
adjust_&&&&_annually-nhrcmc-@@@@_YYYY.csv
ここで、 &&&&：要素名
@@@@：メンバー名
YYYY：4 桁の西暦
pythonのpandasライブラリを用いて、20年分の20ファイルを読み込んで1つのデータフレームにします。
# ライブラリのインポートと作業ディレクトリ変更 import os import pathlib import numpy as np import pandas as pd os.chdir("/mnt/c/Users/hoge/JWP9/") # ファイル名に"t2mx_ov25"を含むcsvのリストを作成します spa_files_df = pd.DataFrame({"file_path":list(pathlib.Path("./dias/data/jmagwp/gwp9/data/csv/SPA/annually").glob("*t2mx_ov25*"))}) # 1ファイルずつ読み込み、年情報を追加した上で結合します。 spa_df = pd.DataFrame(index=[]) #入れ物のデータフレームを用意 for i,one_line in spa_files_df.iterrows(): # ファイルリストを１行ずつ処理 tmp_df = pd.read_csv(one_line["file_path"],header=None,names=[0,1,2,3]) # 年ごとのcsvを４列のデータフレームに読み込み date_str = tmp_df.iloc[0,1] # 年の情報を抽出 sub_df = tmp_df.iloc[2:,:] # データ行を切り出し sub_df.columns = tmp_df.iloc[1,:] # 列名を設定 sub_df = sub_df.assign(date_str=date_str) # 年情報の列を追加 spa_df = pd.concat([spa_df,sub_df],axis=0) # 入れ物のデータフレームに追加する # データ型を実数に変更 spa_df["value"] = spa_df["value"].astype("float") # フラッグ列を追加（この後処理する将来予測実験のデータフレームと結合するために必要） spa_df["MEM-FLG"] = 1 # 20年分のデータを、行方向が地点ID、列方向が年のデータフレームに整形する spa_table = spa_df.set_index(["sfc_no","MEM-FLG","date_str"])["value"].unstack(level=2)
将来気候予測実験(SFA)のデータフレーム作成
RCP2.6 シナリオに基づく将来予測計算結果の年集計値は/mnt/c/Users/hoge/JWP9/dias/data/csv/SFA‗rcp26_ens/annually/に、RCP8.5 シナリオに基づく結果の集計値は/mnt/c/Users/hoge/JWP9/dias/data/csv/SFA‗rcp26_ens/annually/に解凍されていますので、現在気候実験(SPA)と同様にデータフレーム化します。
# rcp26シナリオの結果のデータフレーム化 rcp26_files_df = pd.DataFrame({"file_path":list(pathlib.Path("./dias/data/jmagwp/gwp9/data/csv/SFA_rcp26_ens/annually").glob("*t2mx_ov25*"))}) rcp26_df = pd.DataFrame(index=[]) for i,one_line in rcp26_files_df.iterrows(): tmp_df = pd.read_csv(one_line["file_path"],header=None,names=[0,1,2,3,4]) date_str = tmp_df.iloc[0,1] sub_df = tmp_df.iloc[2:,:] sub_df.columns = tmp_df.iloc[1,:] sub_df = sub_df.assign(date_str=date_str) rcp26_df = pd.concat([rcp26_df,sub_df],axis=0) rcp26_df["value_ens"] = rcp26_df["value_ens"].astype("float") rcp26_df["MEM-FLG"] = rcp26_df["MEM-FLG"].astype("int") rcp26_table = rcp26_df.set_index(["sfc_no","MEM-FLG","date_str"])["value_ens"].unstack(level=2) # rcp85シナリオの結果のデータフレーム化 rcp85_files_df = pd.DataFrame({"file_path":list(pathlib.Path("./dias/data/jmagwp/gwp9/data/csv/SFA_rcp85_ens/annually").glob("*t2mx_ov25*"))}) rcp85_df = pd.DataFrame(index=[]) for i,one_line in rcp85_files_df.iterrows(): tmp_df = pd.read_csv(one_line["file_path"],header=None,names=[0,1,2,3,4]) date_str = tmp_df.iloc[0,1] sub_df = tmp_df.iloc[2:,:] sub_df.columns = tmp_df.iloc[1,:] sub_df = sub_df.assign(date_str=date_str) rcp85_df = pd.concat([rcp85_df,sub_df],axis=0) rcp85_df["value_ens"] = rcp85_df["value_ens"].astype("float") rcp85_df["MEM-FLG"] = rcp85_df["MEM-FLG"].astype("int") rcp85_table = rcp85_df.set_index(["sfc_no","MEM-FLG","date_str"])["value_ens"].unstack(level=2)
統合データフレームの作成
３つのデータフレームを統合して保存します。
# データフレームそれぞれのインデックスにメンバー名を追加 spa_tmp = pd.concat([spa_table],keys=["spa"],names=["scenario"]) rcp26_tmp = pd.concat([rcp26_table],keys=["rcp26"],names=["scenario"]) rcp85_tmp = pd.concat([rcp85_table],keys=["rcp85"],names=["scenario"]) # カラム名を統一 spa_tmp.columns=np.arange(0,20) rcp26_tmp.columns = np.arange(0,20) rcp85_tmp.columns = np.arange(0,20) # 統合して一つのデータフレームにする t2mx_ov25_df = pd.concat([spa_tmp,rcp26_tmp,rcp85_tmp]) t2mx_ov25_df = t2mx_ov25_df.reorder_levels([1,0,2]).sort_index() # マルチインデックスの順位を変更して、ソートをかける t2mx_ov25_df.to_pickle("t2mx_ov25_df.pkl") # pickle形式で保存しておく
地点検索用のポイントベクタ作成
「観測地点データ」のcsvには”sfc_no”という地点IDが振られています。このIDがどのアメダス観測地点に対応しているかは、インデックスファイルに載っています。そのまま使ってもよいのですが、直感的に地図からIDを検索できるよう、ポイントベクタを作成します。
# ライブラリのインポート import os import numpy as np import pandas as pd import geopandas as gpd # 作業ディレクトリへ移動 os.chdir("/mnt/c/Users/hoge/JWP9/") #現在再現実験データの読み込み t2mn_df = pd.read_csv("./dias/data/jmagwp/gwp9/data/csv/SPA/annually/adjust_t2mn_ov25_annually-nhrcmc-SPA_1980.csv",skiprows=1) ppsf_df = pd.read_csv("./dias/data/jmagwp/gwp9/data/csv/SPA/annually/adjust_ppsf_max_annually-nhrcmc-SPA_1980.csv",skiprows=1) sndp_df = pd.read_csv("./dias/data/jmagwp/gwp9/data/csv/SPA/annually/adjust_sndpmx_annually-nhrcmc-SPA_1980.csv",skiprows=1) # データフレームの結合と整形（重複行削除とソート） points_df = pd.concat([t2mn_df[["sfc_no","lat","lon"]], ppsf_df[["sfc_no","lat","lon"]], sndp_df[["sfc_no","lat","lon"]]]) points_df = points_df.drop_duplicates() poitns_df = points_df.sort_values("sfc_no") points_df = points_df.reset_index(drop=True) # 要素の有無を判別し、"???_point"列に記載 points_df = points_df.assign(temp_point=points_df["sfc_no"].isin(t2mn_df["sfc_no"])) points_df = points_df.assign(prec_point=points_df["sfc_no"].isin(ppsf_df["sfc_no"])) points_df = points_df.assign(snow_point=points_df["sfc_no"].isin(sndp_df["sfc_no"])) # インデックスファイルの読み込み temp_index = pd.read_csv("./dias/data/jmagwp/gwp9/meta/temp_index.csv") prec_index = pd.read_csv("./dias/data/jmagwp/gwp9/meta/prec_index.csv") snow_index = pd.read_csv("./dias/data/jmagwp/gwp9/meta/snow_index.csv") # インデックスファイルの結合 all_index = pd.concat([temp_index[["sfc_no","sfc_name"]], prec_index[["sfc_no","sfc_name"]], snow_index[["sfc_no","sfc_name"]]]) all_index = all_index.drop_duplicates() all_index = all_index.set_index("sfc_no").sort_index() # ポイントベクタの作成 points_df = points_df.assign(sfc_name=all_index.loc[points_df["sfc_no"],"sfc_name"].values) points_geometry = gpd.points_from_xy(points_df["lon"],points_df["lat"],crs="EPSG:4326") gdf = gpd.GeoDataFrame(points_df[["sfc_no","sfc_name","temp_point","prec_point","snow_point"]],geometry=points_geometry) # 今流行りのfgb形式で保存 gdf.to_file("kansoku_points.fgb",index=False, driver="FlatGeobuf", spatial_index="NO")
QGISに読み込んで凡例をつける
geopandasで作成したポイントベクタをQGISで読み込み、凡例等を設定し、往路ジェクトを作成しておきます。
QGISを立ち上げます
「レイヤ（L）」→「レイヤを追加」→「ベクタレイヤを追加」
データソースマネージャを設定
- ファイル名：
- 観測地点ポイントベクタのfgbファイルを指定する
「追加」をクリック
- ファイル名：
レイヤパネルの「kansoku_point」を右クリックして、レイヤプロパティを立ち上げ、シンボロジタブを選択します
シンボロジタブで「ルールによる定義」を選択し、１行目をダブルクリックします
”Edit Rule”ウィンドウで、気温に関する統計値の存在するポイントの表示を設定します。
- ラベル → ”Temp"
- フィルタ → ”temp_point"
- シンボル ：
- シンボルレイヤ型 → シンプルマーカー
- 色 → ここでは赤系を指定
ルールを追加し、降水に関する統計値、雪に関する統計値についても同様に設定します。
- ラベル：”prec" フィルタ：”prec_point"
- ラベル：”snow" フィルタ：”snow_point"
シンボルを選択することで、それぞれの統計値が存在する地点を表示することができます
「プロジェクト（J)」→「名前を付けて保存」で、プロジェクトを保存しておきます。 （JWP9_kansoku.qgz)
AMeDAS[手稲山口] の夏日の年間日数の将来変化
ここでは、「北海道石狩地方 山口」地点における夏日の年間日数の将来変化の箱ひげ図を作成してみます。
地点IDの検索・確認
観測地点ポイントベクタから、「北海道石狩地方 山口」地点の地点ID（"sfc_no"）を確認します。
QGISを立ち上げ、「プロジェクト」→「開く」で”JWP9_kansoku.qgz”を開き、”temp"を選択します
「北海道石狩地方 山口」地点にズームし、「地物情報を表示」アイコンを選択してから、「北海道石狩地方 山口」地点をクリックすると、右側に地物情報が表示されます。
「北海道石狩地方 山口」のsfc_noは14116であることが確認できます
将来変化グラフ（箱ひげ図）を作成
初めに作成した”t2mx_ov25_df”からsfc_noが14116のデータを取り出し、plotlyで箱ひげ図を描きます。
# ライブラリのインポート import os import numpy as np import pandas as pd import plotly.express as px import plotly.graph_objects as go # 作業ディレクトリの変更とデータベースの読み込み os.chdir("/mnt/c/Users/hoge/JWP9/") t2mx_ov25_df = pd.read_pickle("t2mx_ov25_df.pkl") # sfc_noが14116のデータを切り出し、plotlyに入力できる形に整形 data_df = t2mx_ov25_df.loc[("14116",),:].reset_index(level=1,drop=True).T[["spa","rcp26","rcp85"]] # 箱ひげ図の作成 fig = go.Figure() fig.add_trace(go.Box(y=data_df["spa"],marker_color="black",name="historical")) fig.add_trace(go.Box(y=data_df["rcp26"],marker_color="blue",name="rcp26") ) fig.add_trace(go.Box(y=data_df["rcp85"],marker_color="red",name="rcp85")) fig.update_traces(boxpoints='all', jitter=0) fig.update_yaxes(title="夏日日数（日/年）") fig.update_layout( width=400,height=300, margin=dict(t=10, b=20, l=10, r=40)) fig.show()
表示されます。
ほかの地点、ほかの統計値についても、同じ方法で箱ひげ図を作ることができます。 ぜひお近くのアメダス地点について将来予測を見える化してみてくださいね～。
※本記事では「気候予測データセット2022 ②日本域気候予測データ」を使用しました。同データセットは、気象庁気象研究所が開発した気候モデルを利用して、文部科学省気候変動リスク情報創生プログラム（RCP8.5シナリオ）及び統合的気候モデル高度化研究プログラム（RCP2.6シナリオ）において計算されたデータを元に作成されたものです。
