GCS上のERA5データをColabで解析する

# はじめに
Google Cloud Storage (GCS) 上にある大気再解析データERA5をGoogle Colaboratory (Colab)で解析する方法についての記事です。

一般にERA5データを利用する場合は、[Climate Data Store (CDS)](https://cds.climate.copernicus.eu/)からサーバー等にデータをダウンロードします。あらかじめ研究室サーバー等でデータを整備していない場合[^1]は、各個人で大量のデータをダウンロードする必要があり時間とディスク容量が必要になります。

別の方法として、GCSのパブリックデータセットである[Analysis-Ready, Cloud Optimized (ARCO) ERA5 Data](https://console.cloud.google.com/storage/browser/gcp-public-data-arco-era5)を利用することができます。このデータセットはZarr形式で整備されているため、ColabからアクセスすることでダウンロードすることなくXarrayからすぐに開くことができます[^2]。

# Google Colab
Googleアカウントがあれば、Colabは無料で利用できます。ただし、本格的に利用する際には、バックグラウンド実行が可能で多くのメモリが利用できるColab Pro (月額1,179円、2024/1現在)をお勧めします。

ARCO ERA5 Data は GCSのus-central1 (アイオワ) リージョンにあります。なるべくデータと同じリージョンでColabを利用したいところですが、幸いなことにus-central1が割り当てられることが多いようです。割り当てられたリージョンは以下の記事の方法で調べることができます。

https://qiita.com/kasagon/items/2dd5475ce00f16333a6a

# 実行例
ここでは特定の緯度経度の可降水量の日別の時系列データをnetCDF形式で保存する例を紹介します。主に公式で用意されている[ノートブック](https://github.com/google-research/arco-era5/blob/main/docs/0-Surface-Reanalysis-Walkthrough.ipynb)を参考にしています。

## 準備
まず、GCSに接続するためにユーザー認証をします
```python
from google.colab import auth 
auth.authenticate_user()
```

また、データを保存するため、Google Driveをマウントします
```python
from google.colab import drive
drive.mount('/content/drive')
```

`/content/drive`以下のパスについては、左側の「ファイル」パネルから確認することができます。

XarrayでZarrを読み込むためには、以下を実行しておく必要があります。
```
!pip install xarray[io]
```

## データセットの選択
データセットの説明は[公式ページ](https://cloud.google.com/storage/docs/public-datasets/era5?hl=ja
)に書かれています。
- `raw/` ECMWF から取り込まれたソースデータ(nc)
- `co/` Cloud Optimized、ガウス格子、Zarr形式
- `ar/` Analysis-Ready、緯度経度格子、Zarr形式

ここではデータ解析に便利な`ar`を利用することにします。`ar`に含まれるデータセットの一覧は`fsspec`を用いて確認できます
```Python
import fsspec

fs = fsspec.filesystem('gs')
fs.ls('gs://gcp-public-data-arco-era5/ar/')
```

```
['gcp-public-data-arco-era5/ar/1959-2022-1h-240x121_equiangular_with_poles_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-1h-360x181_equiangular_with_poles_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-128x64_equiangular_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-128x64_equiangular_with_poles_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-1440x721.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-240x121_equiangular_with_poles_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-512x256_equiangular_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-64x32_equiangular_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-6h-64x32_equiangular_with_poles_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-full_37-1h-0p25deg-chunk-1.zarr-v2',
 'gcp-public-data-arco-era5/ar/1959-2022-full_37-6h-0p25deg-chunk-1.zarr-v2',
 'gcp-public-data-arco-era5/ar/1959-2022-full_37-6h-0p25deg_derived.zarr',
 'gcp-public-data-arco-era5/ar/1959-2022-wb13-6h-0p25deg-chunk-1.zarr-v2',
 'gcp-public-data-arco-era5/ar/1959-2023_01_10-6h-64x32_equiangular_conservative.zarr',
 'gcp-public-data-arco-era5/ar/1959-2023_01_10-full_37-1h-0p25deg-chunk-1.zarr',
 'gcp-public-data-arco-era5/ar/1959-2023_01_10-full_37-1h-1440x721.zarr',
 'gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3']
 ```
オリジナルの水平0.25°格子・hourlyデータの他に、解像度を荒くしたデータセットがあります[^3]。
ここでは水平1°格子データを利用することにします。

## ファイルの読み込みとGoogle Driveへの保存
データセットを開きます。
```Python
import xarray as xr

ds = xr.open_zarr(
    'gs://gcp-public-data-arco-era5/ar/1959-2022-1h-360x181_equiangular_with_poles_conservative.zarr',
    chunks={'time': 24},
    consolidated=True,
)
ds
```
`chunks`や`consolidated`オプションをつけているため、この段階ではデータはメモリ上に読み込まれていません(lazy load)。

データセットから1年分の可降水量を読み込みます。このデータセットは空間方向にチャンク分割されていないので、1地点を`.sel()`してから`.compute()`しても`.compute()`してから1地点を抽出しても同じように全球のデータが読み込まれます。

```python
ds_t = ds["total_column_water_vapour"] \
    .sel(time=slice("2000-01-01", "2000-12-31")) \
    .compute() 
```    
試したところ、この処理は20秒程度で完了しました。水平1°格子データとはいえ、hourlyデータをこの時間で読み込めるのは速いと思います。

特定の地点の選択したのち日平均し、ファイルに保存します。出力先ディレクトリ(`/content/drive/MyDrive/ERA5`)はあらかじめ作成しておきます。
```python
ds_t.sel(latitude=5, longitude=180, method="nearest") \
    .resample(time='1D').mean() \
    .to_netcdf("/content/drive/MyDrive/ERA5/output.nc")
```
すでにデータはメモリ上にあるので、1地点だけであればこの処理は数秒で終わりました。

[^1]: 本来は大学や研究所単位でデータを整備するのが効率が良いのでしょうが、どこでもそのような取り組みが行われているわけではありません。また、オンプレの解析用サーバー自体を持たない気象学の研究室もあるようです。日本だと軽視されがちですが、環境の整備は重要だと私は思います。
[^2]: データをロードする際には時間がかかりますが、CDSのように待ち時間は発生しません。
[^3]: 元々はAI数値予報モデルの学習に用意されたものだと思われます。