「全球及び日本域150年連続実験データ」を可視化するその３ー統合プログラム領域20km150年連続実験データの時系列データをQGISでアニメーション表示する。

はじめに

前回は、150年連続実験データの時系列データをgifアニメーションにする方法を紹介しました。

cci-labo.hateblo.jp

ただ、やはり、興味のあるエリア（ROI）を拡大したり、動画を止めたり、ほかの情報と重ね合わせたい！

それなら、時間軸を持ったnetCDF形式で出力し、QGISにメッシュレイヤとして読み込むことで実現できそう。

ということで、やってみましょう。

netCDFの出力にはXarrayを使います。なければインストールしておいてください。

cci-labo.hateblo.jp

時系列netCDFの作成

まずは、各年8月上旬の地上気温の平均値の3次元配列（time,y,x）を作成します。

# ライブラリをインポート
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)

# 空のNumpyArrayを用意する
temp_arr = np.zeros([149,155,191],np.float32)

# 1年ずつ期間平均し、配列に代入
for h_index,h_row in griblist_df.iterrows():
    ds = xr.open_dataset(h_row["grib_path"], engine="cfgrib")
    da = ds["t"]
    da_eary = da.sel(time=slice("{}-08-01 1:00:00".format(h_row["year"]), "{}-08-11 0:00:00".format(h_row["year"])))
    da_eary_mean = da_eary.mean(dim="time")
    temp_arr[h_index,:,:]  = da_eary_mean.data

netCDFファイルの作成

xarray.Datasetを作成し、投影パラメータを設定してnetCDF形式のファイルへ出力します。

# x軸、y軸の設定
xx_vec = np.arange(191) * 20000 + 10000 -1930000.1287495417
yy_vec = np.arange(155) * -20000 - 10000 + 5833346.794147097

# xarray.Datasetを作成
temp_ds = xr.Dataset(
    {
        "temperature": (["time", "y","x" ], temp_arr[:,::-1,:])
    },
    coords={
        "y": ("y", yy_vec),
        "x": ("x", xx_vec),
        "time": pd.date_range("1951-08-01",freq="AS-AUG",periods=149),
        "reference_time": pd.Timestamp("1950-08-01"),
    },
)

# 投影情報を設定
temp_ds = temp_ds.rio.write_crs("+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")
temp_ds = temp_ds.rio.write_transform()

# netCDFへ保存
temp_ds.to_netcdf("NHRCM20_RCP85_tmp_early_Aug.nc")