「全球及び日本域150年連続実験データ」を可視化する その２ー統合プログラム領域20km150年連続実験データの時系列変化を可視化する。

はじめに

「全球及び日本域150年連続実験データ」は気候変動がいつの段階でどのくらいまで進むか？その時間推移を把握することができます。

その変化を視覚的に捉えるため、マップの時系列アニメーションを作成してみましょう。

時系列ラスタの作成

前回入手した、RCP8.5シナリオの8月の149年分のデータを集計し、8月上旬（1日～10日）の気温の時系列GeoTIFFを作成します。

1951/8/1~1951/8/10の平均マップ

まずは1951年8月上旬（1日～10日）について、地上気温の平均マップを作ってみます。

# ライブラリをインポート
import os
import pathlib
import numpy as np
import pandas as pd
import datetime
import docker
from osgeo import gdal, gdalconst, gdal_array 
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# "wgrib_docker”コンテナを起動する
client = docker.from_env()
wgrib_cont =client.containers.get("wgrib_dock_mntd")
wgrib_cont.start()

# 処理ファイルのパスを設定
ctl_path = pathlib.Path("/mnt/c/Users/hoge/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.ctl")
grib_path = pathlib.Path("/mnt/c/Users/hoge/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.grib")

# ctlファイルのVAR行を特定
substr_df = pd.read_fwf(ctl_path,widths=[10],header=None,skip_blank_lines=False).fillna("#")
n_skip = substr_df.loc[substr_df[0].str.startswith("var")].index.values[0] + 1

# VARの内容を読み込み
parm_df = pd.read_fwf(ctl_path,header=None,skiprows=n_skip,nrows=17,
　　　　　　　　　　colspecs=[(0,8),(11,13),(21,100)],
　　　　　　　　　　names=["varname","parm_no","description"],
　　　　　　　　　　dtype={"varname":str, "parm_no":int, "description":str})

# gribファイルのヘッダ情報を読み込み
res_code,res_out = wgrib_cont.exec_run('wgrib -min {}'.format(str(grib_path)), 
                                        workdir=str(grib_path.parents[2]),stderr=False)

# レコード情報を整理
reco_df =  pd.DataFrame({"v_info":res_out.decode().split("\n")})
reco_df = reco_df.loc[reco_df["v_info"].str.match("^\d")]
reco_df["band_no"] = reco_df["v_info"].str.extract("(^\d*):.*",expand=False).astype(int)
reco_df["date_time"] = pd.to_datetime(
             ("19" + reco_df["v_info"].str.extract(".*:d=(\d*-d*):.*",expand=False)),
             format="%Y%m%d%H-%M")
reco_df["parm_str"] = reco_df["v_info"].str.extract(".*:d=\d*-\d*:([A-Z]*):.*",expand=False)
reco_df["parm_no"] = reco_df["v_info"].str.extract(".*kpds5=(\d*)",expand=False).astype(int)
reco_df = reco_df.reset_index(drop=True)

# parm_dfを結合する
reco_df = reco_df.merge(parm_df)

# 8/1~8/10の地上気温レコードを絞り込む
reco_tmp = reco_df.loc[
    (reco_df["varname"] == "tmp") &
    (reco_df["date_time"] >= datetime.datetime(1951,8,1,0,0)) & 
    (reco_df["date_time"] <= datetime.datetime(1951,8,11,0,0))
    ] 

# gribファイルをGDALデータセットとして読み込む
grib_ds = gdal.Open(str(grib_path),gdalconst.GA_ReadOnly)

# 1レコードずつ読み込み、平均値を求める
tmp_arr_sum = np.zeros((155,191))

for h_ind,h_row in reco_tmp.iterrows():
    tmp_arr_sum += grib_ds.GetRasterBand(h_row["band_no"]).ReadAsArray()

tmp_arr_mean = tmp_arr_sum/reco_tmp.shape[0]

# Cartpyで表示してみる。
fig = plt.figure(figsize=(10, 5))
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=135.0,
                       central_latitude=0,standard_parallels=(30.0,60.0)))
ax.imshow(tmp_arr_mean,origin="upper",
          transform=ccrs.LambertConformal(central_longitude=135.0,
                       central_latitude=0,standard_parallels=(30.0,60.0)),
          extent=[-1930000, 1890000, 2733347, 5833347])
ax.coastlines()
ax.gridlines()

cartpyの使い方は以下を参考にさせていただいています。

yyousuke.github.io

処理の関数化

この処理を149年分繰り返し、149バンドのGeoTIFFを作成します。 ループを回すために、まずはそれぞれの処理を関数化しておきます。

# ライブラリをインポート
import os
import pathlib
import numpy as np
import pandas as pd
import datetime
import docker
from osgeo import gdal, gdalconst, gdal_array
from osgeo import osr

# ctlファイルを読み込みparameter_dfを返す関数
def res_parameter_df(ctl_path):
    # 行を特定
    substr_df = pd.read_fwf(ctl_path,widths=[10],header=None,skip_blank_lines=False).fillna("#")
    n_skip = substr_df.loc[substr_df[0].str.startswith("var")].index.values[0] + 1

    # 読み込み
    parm_df = pd.read_fwf(ctl_path,header=None,skiprows=n_skip,nrows=17,
                                         colspecs=[(0,8),(11,13),(21,100)],
                                         names=["varname","parm_no","description"],
                                         dtype={"varname":str, "parm_no":int, "description":str})

    # parm_dfを戻す
    return parm_df

# gribファイルのヘッダ情報を読み込み、records_dfを返す関数
def res_records_df(grib_path,year,parm_df):
    # ヘッダ情報を読み込み
    res_code,res_out = wgrib_cont.exec_run('wgrib -min {}'.format(str(grib_path)), 
                                    workdir=str(grib_path.parents[2]),stderr=False)

    # レコード情報を整理
    reco_df =  pd.DataFrame({"v_info":res_out.decode().split("\n")})
    reco_df = reco_df.loc[reco_df["v_info"].str.match("^\d")]
    reco_df["band_no"] = reco_df["v_info"].str.extract("(^\d*):.*",expand=False).astype(int)
    reco_df["date_time"] = pd.to_datetime((str(h_row["year"])[:2] + 
             reco_df["v_info"].str.extract(".*:d=(\d*-\d*):.*",expand=False)),
             format="%Y%m%d%H-%M")
    reco_df["parm_str"] = reco_df["v_info"].str.extract(".*:d=\d*-\d*:([A-Z]*):.*",expand=False)
    reco_df["parm_no"] = reco_df["v_info"].str.extract(".*kpds5=(\d*)",expand=False).astype(int)
    reco_df = reco_df.reset_index(drop=True)

    #parm_dfを結合してから戻す
    return reco_df.merge(parm_df)

# gribファイルから特定要素、期間のデータを切り出し、平均値の配列を返す関数
def res_mean_arr(grib_ds,reco_df,varname,datetime_start,datetime_end):

    reco_sub = reco_df.loc[
        (reco_df["varname"] == varname) & 
        (reco_df["date_time"] >= datetime_start) & 
        (reco_df["date_time"] <= datetime_end)
        ] 

    # 1レコードずつ読み込み、平均値を求める
    arr_sum = np.zeros((155,191))

    for h_ind,h_row in reco_sub.iterrows():
        arr_sum += grib_ds.GetRasterBand(h_row["band_no"]).ReadAsArray()

    arr_mean = arr_sum/reco_sub.shape[0]

    # 平均値の配列を戻す
    return arr_mean

時系列GeoTIFFの作成

149年分のデータについて、8月上旬（1日～10日）の地上気温の平均マップを作成し、1つのGeoTIFFにします。

# 前のコードブロックのライブラリの読み込みと関数を設定しておく。

# "wgrib_docker”コンテナを起動する
client = docker.from_env()
wgrib_cont =client.containers.get("wgrib_dock_mntd")
wgrib_cont.start()

# 地理座標系のパラメータを設定
ul_x = -1930000.1287495417
ul_y = 5833346.794147097
x_res = 20000
y_res = -20000
x_size = 191
y_size = 155

outRasterSRS = osr.SpatialReference()
outRasterSRS.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")

# 処理するファイルの一覧を作成する
os.chdir("/mnt/d/Users/kazuh/ONEFIFTY")
griblist_df = pd.DataFrame({"grib_path":list(pathlib.Path("./").glob("**/*.grib"))})
ctllist_df = pd.DataFrame({"ctl_path":list(pathlib.Path("./").glob("**/*.ctl"))})

griblist_df["dir_path"] = griblist_df["grib_path"].map(lambda x: x.parent.absolute())
ctllist_df["dir_path"] = ctllist_df["ctl_path"].map(lambda x: x.parent.absolute())
griblist_df = griblist_df.merge(ctllist_df)

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)

# 空の出力ファイルを用意する
name = "NHRCM20_RCP85_tmp_early_Aug.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name,x_size,y_size,149,gdal.GDT_Float64)
output.SetProjection(outRasterSRS.ExportToWkt())
output.SetGeoTransform((ul_x, x_res, 0, ul_y, 0, y_res))

# ファイル一覧を1行ずつ処理
for h_index,h_row in griblist_df.iterrows():
    #records_dfを作成
    parameter_df = res_parameter_df(h_row["ctl_path"].absolute())
    records_df = res_records_df(h_row["grib_path"].absolute(),h_row["year"],parameter_df)

    # gribファイルをGDALデータセットとして読み込む
    grib_ds = gdal.Open(str(h_row["grib_path"].absolute()),gdalconst.GA_ReadOnly)

    # 8/１ 1:00～8/11 0:00の地表気温の平均値マップを作成
    tmp_arr = res_mean_arr(grib_ds,records_df,"tmp",
                            datetime.datetime(h_row["year"],h_row["month"],1,0),
                            datetime.datetime(h_row["year"],h_row["month"],11,0))

    # GeiTIFFのバンドに保存
    outband = output.GetRasterBand(h_index+1)
    outband.WriteArray(tmp_arr.astype("float64"))

    # バンド名を設定
    outband.SetDescription("{}-{:02d}".format(h_row["year"],h_row["month"]))
    outband.FlushCache()

# 出力ファイルを閉じる
output = None

アニメーション化

matplotlibを用いてgifアニメーションを作成します。

# ライブラリをインポート
import os
import pathlib
from osgeo import gdal, gdalconst, gdal_array
import matplotlib.pyplot as plt
from matplotlib import animation, rc
import cartopy.crs as ccrs

%matplotlib widget
plt.rcParams["font.family"] = 'TakaoGothic'

# geotiffを読み込む
os.chdir("/mnt/d/Users/kazuh/ONEFIFTY")
geotiff_ds = gdal.Open("NHRCM20_RCP85_tmp_early_Aug.gtiff",gdalconst.GA_ReadOnly)
outband = geotiff_ds.GetRasterBand(1)

# マップを更新する関数
def update(i):
    if i != 0:
        plt.cla()   
    outband = geotiff_ds.GetRasterBand(i+1)
    ax.set_extent([110, 160, 19, 50]) 
    ax.imshow(outband.ReadAsArray(),origin="upper",vmin=5.,vmax=40.,
              transform=ccrs.LambertConformal(central_longitude=135.0,central_latitude=0,
                                              standard_parallels=(30.0,60.0)),
              extent=[-1930000, 1890000, 2733347, 5833347],cmap='afmhot')
    ax.coastlines(resolution="10m")
    ax.gridlines()
    ax.set_title('地上気温 {}'.format(outband.GetDescription()))
    return fig

# アニメーションを作成
fig = plt.figure(figsize=(10, 5))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_extent([110, 160, 19, 50])
im = ax.imshow(outband.ReadAsArray(),origin="upper",vmin=5.,vmax=40.,
                  transform=ccrs.LambertConformal(central_longitude=135.0,central_latitude=0,
                                                    standard_parallels=(30.0,60.0)),
                  extent=[-1930000, 1890000, 2733347, 5833347],cmap='afmhot')
ax.coastlines(resolution="10m")
ax.gridlines()
cb = plt.colorbar(im)
cb.ax.set_title('[℃]')
ax.set_title('地上気温 {}'.format(outband.GetDescription()))

ani = animation.FuncAnimation(fig, update, 148, interval=600)
ani.save("NHRCM20_RCP85_tmp_early_Aug.gif", writer = 'imagemagick')

アニメーションができました!

他の月、気象要素についても同様に作成することができます。

「全球及び日本域150年連続実験データ」の20km150年連続実験データのアニメーションを作る方法を紹介しました。試してみてくださいね～。

※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。