GraphCastの気象予測データをCesiumで可視化してみた

# はじめに

近年、異常気象による災害の被害が増加していることから、早期に異常気象を発見できるように長期間先を正確に気象予測する技術が重要になってきています。現在の気象予測の精度は当日や明日など直近の予測であれば高精度で予測できるようになってきましたが、週間予報のような中長期予報では最近は外れることが多いように感じます。

そんな中で 1 か月ほど前、Google DeepMind が従来のシステムよりも高精度かつ高速に最大 10 日間先の気象の予測を行う AI モデル「[GraphCast](https://deepmind.google/discover/blog/graphcast-ai-model-for-faster-and-more-accurate-global-weather-forecasting/)」を公開しました。この GraphCast は予測が難しいとされる異常気象の予測にも強く、 2023 年 9 月に北米に上陸した大型ハリケーン「Lee」について、従来の予測システムでは上陸の約 6 日前に上陸を予測したのに対し、GraphCast は約 9 日前に上陸を予測し的中したと発表されました。

この GraphCast の発表を見て予想された気象データを分かりやすく伝えるにはどのように可視化するのがよいのだろうかという点に興味を持ちました。そこで、弊社でよく使用している 3D 地球儀ライブラリ Cesium を使用して GraphCast で予測された気象データの可視化を試してみました。

# GraphCast による気象予測

GraphCast は DeepMind 公式から実行用のソースコードと学習済みモデルが公開されています。

- https://github.com/google-deepmind/graphcast

またヨーロッパ中期予報センター（ECMWF）から、入力する気象データを API により自動で取得して GraphCast で推論を実行してくれるツールも公開されています。

- https://github.com/ecmwf-lab/ai-models

今回はこちらの「ecmwf-lab/ai-models」を使用して GraphCast による気象予測を行ってみます。

事前準備として入力する気象データの取得先である [Copernicus Climate Data Store (CDS)](https://cds.climate.copernicus.eu/)にアカウント登録を行います。

動作環境は Colab Pro の GPU A100 の環境で実行しました。
Colab ノートブック上で必要なライブラリをインストールします。

    !pip install ai-models
    !pip install ai-models-graphcast  # Install details at https://github.com/ecmwf-lab/    ai-models-graphcast
    !git clone https://github.com/ecmwf-lab/ai-models-graphcast.git
    !cd ai-models-graphcast && pip install -r requirements-gpu.txt -f https://storage.googleapis.   com/jax-releases/jax_cuda_releases.html

入力する気象データの日時を指定して GraphCast による気象予測を実行します。

    !ai-models --input cds --date 20230813 --time 0000 --download-assets graphcast

途中で CDS からデータを取得するために API キーの入力が求められるので、表示された URL を開いて API キー等をコピペして入力します。

しばらく待つと推論が完了して Grib 形式のファイルが出力されます。デフォルトでは 6 時間間隔で 10 日間の気象予測が行われます。

以下は 2023/12/7 の気象データを入力して GraphCast で予測された 2023/12/9 の気圧のデータを QGIS でカラーマップとして表示した画像です。
![GraphCastで予測された2023/12/9の気圧](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/da8ee4fe-c702-d82b-df9f-e076e4309cc3.png)

以下は 2023/12/9 に実際に観測された気圧データです。
![実際に観測された2023/12/9の気圧](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/a1c7e2b6-b483-9ad0-d766-8e07d7d9047d.png)

比較してみるとかなり高精度で予測出来ていることが分かります。

気象データは上記のカラーマップのように静止画で表されることも多いですが、風のように動きがあるものや時系列での変化を表現する場合にはアニメーションによる可視化が効果的な場合があります。そこで、今回は GraphCast で予測した気象データについてアニメーションによる可視化を試してみました。

# 多次元風速アニメーション

Cesium での風速の可視化については、WebGL により風速をパーティクルアニメーションとして描画することで可視化する手法が発表されており、ソースコードも公開されています。

- https://cesium.com/blog/2019/04/29/gpu-powered-wind/

そこでこちらのソースコードを使用して、GraphCast で予測された気象データから Cesium での風速の可視化を試してみました。
実際に GraphCast で予測された風速データを表示してみると以下のようになります。

![GraphCast予測データの風速パーティクルアニメーションの例](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/df7c2a28-e02b-cf9d-a191-6ea78f31ac2a.png)

こちらの実装では多次元データの読み込みと表示にも対応しているので、次元ごとに高さを変えて立体的に表示することができないか試してみました。次元ごとに高さを変えることで次元による風の動きの違いが見やすくなることが期待できます。

元のソースコードでは地形との衝突判定のためにすべてのパーティクルに対して一律で高さを設定する機能があります。この高さを設定している箇所で次元ごとに高さを変えるように実装すれば実現できそうだと考えていました。以下が変更した箇所です。

- 3D-Wind-Field/Cesium-3D-Wind/glsl/segmentDraw.vert

      vec3 convertCoordinate(vec3 lonLatLev) {
      ...

      float N_Phi = a / sqrt(1.0 - e2 * sinLat * sinLat);
      float h = particleHeight + ((lonLatLev.z - 1.0) * 3000.0); // 次元ごと（lonLatLev.z）に高さを変える

      ...

高さ変更前と次元ごとに高さを変更して表示させたスクリーンショットを以下に示します。分かりやすいように次元ごとに色分けしています。

- 高さの設定変更前
  ![高さ変更前の風速パーティクルアニメーションの例](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/410a4f94-79a7-3706-2a0e-05e0d96e98a6.png)
  
- 次元ごとに高さを変更
  ![次元ごとに高さを変えて表示させた風速パーティクルアニメーションの例](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/86b58804-4fbc-357c-0007-ee1d40025cf0.png)

低い順から 白(約 1400m) 赤(約 4400m) 黄(約 7400m) の高さになっているはずなのですが、いまいち立体的に見えません。高さの設定を変更する前のスクリーンショットと比較してもそこまで変わっていないように思えます。

一応高さはパーティクルの表示位置に影響しているようなのですが、高さから地形に隠れるかどうかの判定は行うものの、パーティクル間の重なりについては計算されていないようでした。そのため本来赤色のパーティクルの下にあるはずの白色のパーティクルが上に重なっているように見えるなどして高さ方向の重なりが分からなくなってしまっています。また、パーティクルは地球の表面上のみに表示されていることや、カメラよりも高い位置にあるパーティクルは表示されない仕様になっていることも立体的に見えない要因となっています。

- 横向きの視点からみた場合のスクリーンショットですが、地球の表面上からはみ出たパーティクルは表示されなくなっています。
  ![横向きの視点からみた場合のスクリーンショット](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/11e5882f-8412-46dd-4c4f-23954834c7b5.png)

これらの問題を修正するとなると時間も労力が大きくかかりそうであったため、今回は立体的な風速データアニメーション表示はあきらめました。

# CS 立体図を用いた気圧可視化アニメーション

風速以外に可視化したら面白そうな気象データがないか探した結果、気圧を可視化してみることにしました。一般的な天気図では気圧は等高線として表されますが、専門知識がないとぱっと見ただけではよくわからないことがあります。そこで、地形判読に用いられる CS 立体図で気圧を可視化することで、よく言われる「気圧の谷」などの気圧の高低差が分かりやすく表現できるのではと考え試してみました。

CS 立体図の表示には「[FOSS4G japan 2020 online](https://www.osgeo.jp/events/foss4g-2020/foss4g-2020-japan-online)」で発表した、Geotiff 画像を地図タイル化して Cesium での表示時に WebGL で CS 立体図を描画させる手法を使用しました。（[Geotiff.js で始めるリアルタイム演算 in foss4g japan 2020 online](https://www.slideshare.net/makinux7/geotiffjs-in-foss4g-japan-2020-online)）

まず、GraphCast で予測したデータを Geotiff に変換してタイル化します。
Geotiff 画像 への変換には QGIS を使用しました。この Geotiff 画像を元の値を保持したままタイル化するのに [tmizu23 様が公開しているタイル生成プログラム](https://tmizu23.hatenablog.com/entry/20121107)を使用させていただきました。

作成した Geotiff タイルを Cesium 上で読み込みリアルタイムレンダリングして CS 立体図を画像レイヤーで表示させます。Geotiff タイルの値の読み込みには、Javascript で Geotiff 画像を読み込むライブラリ Geotiff.js を使用しています。WebGL での CS 立体図のリアルタイムレンダリングは [frogcat 様の標高 PNG タイルから CS 立体図を表示する実装](https://qiita.com/frogcat/items/7e91d3070a7a8d3e2c94)を参考にさせていただきました。

Geotiff タイルをパラメータを調節して CS 立体図として表示したスクリーンショットが以下の画像です。

![Graphcastの気圧予測データをCS立体図として表示](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/c350cca4-ccf1-09ed-d0d9-4ec36432813c.png)

上記の CS 立体図に気圧の値でカラーマップを描画したものを重ねて表示させます。紫色～青色が気圧が低く、黄色～オレンジ色が気圧が高いところになります。

![Graphcastの気圧予測データのCS立体図に気圧カラーマップを重ねて表示](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/da381858-a856-661c-e986-c5d4dcf2b46f.png)


元の気象予測データの解像度が低いので多少表示が変なところもありますが、気圧の起伏が確認できる画像になっています。

さらに、この気圧を表示した CS 立体図を時系列でアニメーションさせてみます。Geotiff タイルを予測データの時間ごとに作成して、Cesium 上で表示レイヤーを切り替えていくことでアニメーションさせます。レイヤー切り替えの際にレイヤーの表示処理を実行すると、タイルの読み込みとレンダリングに時間がかかってしまいアニメーションになりませんでした。そのため、アニメーションに使用する全てのレイヤーを表示順に重ねて表示させた状態にしてから上から順番に消していく、古典的なパラパラ漫画のような手法で実装しました。

以下は 2023/1/10 から 10 日間の気圧予測データをアニメーションさせたものです。

![cs_animation.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/73160ca1-1f7d-e2ae-d97d-84f5a7e8e0c7.gif)

気圧の変化が分かりやすく表示できているのではないでしょうか。

Graphcast の気象予測データが解像度が低めだったので既存のもう少し解像度が高いデータでも試してみました。[京都大学の生存圏データベースで公開している気象庁データ](http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-netcdf.html)から 2023 年 8 月の台風 7 号の気圧データを可視化しました。

![2023年8月の台風7号の気圧データのCS立体図](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/6af901fa-6441-bd90-1390-6fef9258b371.png)

アニメーションさせると以下のようになります。

![2023年8月の台風7号の気圧データのCS立体図アニメーション](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/300835e1-7214-25c2-ec55-eaf3d18507a8.gif)

[frogcat 様が WebGL で等高線を表示させるコード](https://github.com/frogcat/leaflet-tilelayer-glue)を公開されていたのでそちらを参考に等高線も表示させてみました。

![2023年8月の台風7号の気圧データのCS立体図＋等高線のアニメーション](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/9e62eed2-936b-40b9-af48-7705c68050f5.gif)

進んでいくにつれて台風が徐々に勢力を落としていくのが分かります。また、勢力の強いうちは台風の目のようなものも確認できます。

## おまけ：気圧を 3D で可視化してみる

この CS 立体図で可視化した気圧を見ているうちに、気圧を 3D で立体的に可視化したらどうなるのだろうと思いついたのでやってみました。

[Cesium terrain builder](https://github.com/geo-data/cesium-terrain-builder)を使用して台風 7 号の気圧の Geotiff 画像から地形タイルを作成して Cesium に表示してみます。通常の地形よりも起伏が小さいので高さ方向のスケールを 100 倍程度に拡大して表示させました。

![地形タイルで立体的に可視化した気圧 1](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/988ee60f-e510-dd75-f825-7ec6fae36f09.png)
![地形タイルで立体的に可視化した気圧 2](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/690b9709-16d2-8b20-957c-1a6336291004.png)
![地形タイルで立体的に可視化した気圧 3](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58671/dbc480f9-d56a-191e-6858-2a85469c7673.png)

台風の部分だけぽっかり穴が空いているような状態になりました。今回はできませんでしたが 3D 表示をアニメーションさせても面白そうです。

# おわりに

GraphCast で予測した気象データを Cesium により可視化させてみました。普段天気予報などで何気なく見ていた天気図を CS 立体図などの別の表現手法で可視化してみることで、今までとは異なる視点で見ることができ興味深かったです。

今回、OSS の GIS ツールを使うことで思い付きで始めた気象の可視化も比較的簡単に実装することができました。OSS ツールの開発者の皆様ありがとうございました！