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

# はじめに

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

https://www.eorc.jaxa.jp/EARTHCARE/research_product/registration_j.html


今回使用するのはJAXAレベル2(L2)標準プロダクトのうち、CPR単体雲プロダクト(CPR_CLP)です。[JAXA L2標準プロダクトリスト](https://www.eorc.jaxa.jp/EARTHCARE/data/L2prd_list_std_j.html)には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標準プロダクトリスト](https://www.eorc.jaxa.jp/EARTHCARE/data/L2prd_list_std_j.html)のページにあります。

# データの可視化

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

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

```Python
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の利点です。

![スクリーンショット 2025-03-21 22.43.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/c01c74a0-a87b-4b4d-b613-39b8d6d9ac2d.png)

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

```Python
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"})
```

![スクリーンショット 2025-03-21 22.47.19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/b3294c14-5947-4b7e-94b5-e189b87a7c83.png)

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

```Python
ds_Data["cloud_doppler_velocity_1km"].plot()
```

![output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/bffe9e74-97b3-44b6-894f-628af96bd27d.png)

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

```Python
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()
```

![output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/ef28023d-a862-4558-b374-563ae666f50b.png)

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

```Python
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()
```

![output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/e253f056-59c3-49e4-8c78-4162ae4ac533.png)

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

[^1]: 全てのHDF5ファイルがXarrayで読み込めるわけではないです。
[^2]: 変数には`_1km`と`_10km`があります。どちらもグリッド間隔は1kmですが、水平解像度（軌道に沿って平均する範囲）が1kmと10kmで異なります。シミュレーションによる評価では、ドップラー速度の不確実性は500m平均だと約1m/s、5km平均だと0.3-0.4m/s未満となるそうです([Kollias et al. 2023](https://doi.org/10.5194/amt-16-1901-2023))。しかし10kmの平均をすると積乱雲内部の急激な速度変化（上昇流核と下降流など）を潰してしまうため、強い対流の詳細を捉えるには空間分解能とのトレードオフが問題となります​。
[^3]: PDDに書いてあるのでしょうが、調べるのがとても面倒なので確認していません。