EarthCARE/CPR L2aプロダクトをXarrayで眺めてみる

はじめに

2025年03月17日にEarthCARE衛星のレベル2プロダクトが公開されました。搭載されているCloud Profiling Radar (CPR)のドップラー速度を、Pythonを使って描画して眺めてみます。以下のページからデータ利用登録をすることで、誰でもデータを利用することができます。

今回使用するのはJAXAレベル2(L2)標準プロダクトのうち、CPR単体雲プロダクト(CPR_CLP)です。JAXA L2標準プロダクトリストにはCPR_CLPの主要パラメータは「雲マスク・雲粒子タイプ・雲水量・雲氷量・水雲有効粒径・氷雲有効粒径・光学的厚さ」と書かれていますが、そのほかにレーダー反射因子やドップラー速度、再解析由来の気温・湿度など様々な変数が含まれています。

CPR単体エコープロダクト(CPR_ECO)にもドップラー速度は含まれていますが、こちらは雲によるマスクがされていません。また、水平方向のrayはCPR_CLPとCPR_ECOで一致していますが、鉛直方向のbinのheightは一致していません。雲を見たい場合は、基本的にはCPR_CLPを見れば十分なようです。

アルゴリズムの詳細についてはAlgorithm Theoretical Basis Document (ATBD)、データの仕様の詳細についてはProduct Definition Document (PDD)を参照してください。リンクはJAXA L2標準プロダクトリストのページにあります。

データの可視化

EarthCAREのプロダクトはHDF5形式で提供されています。HDF5を扱うPythonのパッケージにはh5pyがありますが、これは比較的低レベルな操作を行うもので、データ解析に利用するにはやや面倒です。今回はXarrayで開くことができるため、これを利用します¹。

HDF5はnetCDFと異なり内部にgroupと呼ばれる階層構造をもつことができるため、Xarrayでopen_datasetをする場合には読み出すgroupを指定する必要があります。今回読み出すCPR_CLPプロダクトでは、経度・緯度・高度・時刻といったいわゆる軸情報は/ScienceData/Geoに、ドップラー速度などの変数は/ScienceData/Dataにあります。

import xarray as xr

h5file = "CPR/2A.CPR_CLP/vBa/2025/03/13/ECA_J_CPR_CLP_2AS_20250313T0213_20250313T0225_04483E_vBa.h5"
ds_Geo = xr.open_dataset(h5file, group="/ScienceData/Geo")
ds_Data = xr.open_dataset(h5file, group="/ScienceData/Data")

Jupyterなどの対話的な環境を利用している場合は、以下のようにDatasetの中身を一覧することができます。これはh5pyに対するXarrayの利点です。

必須ではないですが、軸の名前がphony_dim_*となっているのは分かりづらいので、rayやbinといった中身に変更しておきましょう。

ds_Geo = ds_Geo.rename({"phony_dim_7": "ray", "phony_dim_8": "bin"})
ds_Data = ds_Data.rename({"phony_dim_3": "ray", "phony_dim_4": "bin"})

雲のドップラー速度cloud_doppler_velocity_1kmを描画してみましょう²。Xarrayではとても簡単に描画することができます。

ds_Data["cloud_doppler_velocity_1km"].plot()

どうやらFillValueが-9999のようです³。これをマスクして描画してみましょう。

cloud_doppler_velocity = ds_Data["cloud_doppler_velocity_1km"]
cloud_doppler_velocity = cloud_doppler_velocity_1km.where(cloud_doppler_velocity_1km != -9999.0)

cloud_doppler_velocity.plot()

うまく描画できていそうです。ただ、縦軸と横軸を逆にした方が自然ですね。また、鉛直座標は高度にしたいです。高度をds_Geoから読み出して描画してみましょう。

import matplotlib.pyplot as plt

# pcolormeshに渡す都合上、x軸の座標変数も2次元にする必要がある
ray_2d = np.tile(ds_Geo["ray"].values[:,np.newaxis], (1, ds_Geo["bin"].size))
height_2d = ds_Geo["height"].values

plt.pcolormesh(ray_2d, height_2d, cloud_doppler_velocity, cmap="RdBu_r", vmin=-10, vmax=10)
plt.colorbar()

欲しい図が書けました。右側の下層の雨が強く雨滴の落下速度の大きいところ(5000m以下)では、下向きのドップラー速度が大きいです。上層の氷雲の領域では、ドップラー速度は正と負の両方が見られます。

1. 全てのHDF5ファイルがXarrayで読み込めるわけではないです。 ↩

2. 変数には_1kmと_10kmがあります。どちらもグリッド間隔は1kmですが、水平解像度（軌道に沿って平均する範囲）が1kmと10kmで異なります。シミュレーションによる評価では、ドップラー速度の不確実性は500m平均だと約1m/s、5km平均だと0.3-0.4m/s未満となるそうです(Kollias et al. 2023)。しかし10kmの平均をすると積乱雲内部の急激な速度変化（上昇流核と下降流など）を潰してしまうため、強い対流の詳細を捉えるには空間分解能とのトレードオフが問題となります​。 ↩

3. PDDに書いてあるのでしょうが、調べるのがとても面倒なので確認していません。 ↩