はじめに
Google Cloud Storage (GCS) 上にある大気再解析データERA5をGoogle Colaboratory (Colab)で解析する方法についての記事です。
一般にERA5データを利用する場合は、Climate Data Store (CDS)からサーバー等にデータをダウンロードします。あらかじめ研究室サーバー等でデータを整備していない場合1は、各個人で大量のデータをダウンロードする必要があり時間とディスク容量が必要になります。
別の方法として、GCSのパブリックデータセットであるAnalysis-Ready, Cloud Optimized (ARCO) ERA5 Dataを利用することができます。このデータセットはZarr形式で整備されているため、ColabからアクセスすることでダウンロードすることなくXarrayからすぐに開くことができます2。
Google Colab
Googleアカウントがあれば、Colabは無料で利用できます。ただし、本格的に利用する際には、バックグラウンド実行が可能で多くのメモリが利用できるColab Pro (月額1,179円、2024/1現在)をお勧めします。
ARCO ERA5 Data は GCSのus-central1 (アイオワ) リージョンにあります。なるべくデータと同じリージョンでColabを利用したいところですが、幸いなことにus-central1が割り当てられることが多いようです。割り当てられたリージョンは以下の記事の方法で調べることができます。
実行例
ここでは特定の緯度経度の可降水量の日別の時系列データをnetCDF形式で保存する例を紹介します。主に公式で用意されているノートブックを参考にしています。
準備
まず、GCSに接続するためにユーザー認証をします
from google.colab import auth
auth.authenticate_user()
また、データを保存するため、Google Driveをマウントします
from google.colab import drive
drive.mount('/content/drive')
/content/drive以下のパスについては、左側の「ファイル」パネルから確認することができます。
XarrayでZarrを読み込むためには、以下を実行しておく必要があります。
!pip install xarray[io]
データセットの選択
データセットの説明は公式ページに書かれています。
-
raw/ECMWF から取り込まれたソースデータ(nc) -
co/Cloud Optimized、ガウス格子、Zarr形式 -
ar/Analysis-Ready、緯度経度格子、Zarr形式
ここではデータ解析に便利なarを利用することにします。arに含まれるデータセットの一覧はfsspecを用いて確認できます
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への保存
データセットを開きます。
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地点を抽出しても同じように全球のデータが読み込まれます。
ds_t = ds["total_column_water_vapour"] \
.sel(time=slice("2000-01-01", "2000-12-31")) \
.compute()
試したところ、この処理は20秒程度で完了しました。水平1°格子データとはいえ、hourlyデータをこの時間で読み込めるのは速いと思います。
特定の地点の選択したのち日平均し、ファイルに保存します。出力先ディレクトリ(/content/drive/MyDrive/ERA5)はあらかじめ作成しておきます。
ds_t.sel(latitude=5, longitude=180, method="nearest") \
.resample(time='1D').mean() \
.to_netcdf("/content/drive/MyDrive/ERA5/output.nc")
すでにデータはメモリ上にあるので、1地点だけであればこの処理は数秒で終わりました。