CesiumのShaderを用いた気象データ3Dアニメーションについて

はじめに

今年の11月にFOSS4G 2024 Japanが開催され、コアデイでは弊社も「CesiumのShaderを用いた気象データ3Dアニメーションについて」というタイトルで発表させていただきました。
本記事は上記の発表内容をまとめたものとなっております。

背景

昨年のFOSS4G Advent Calendarでは、弊社はGraphCastの気象予測データをCesiumで可視化してみたという記事で参加させていただきました。こちらの記事を簡単に説明いたしますと、AIが予測した気象データを3D地球儀ライブラリのCesiumで可視化してみた、という内容になっております。こちらの記事の中では、気圧を画像レイヤーでアニメーションさせたり、地形タイルによって気圧を3D化してみたりといったことを行いました。しかしながら2Dの気圧アニメーションは行ったものの3Dでのアニメーションは行っておりませんでした。
そこで気圧データを3Dアニメーションできないか試してみました。

3Dアニメーションの実装

気圧データのようなGridデータの3DアニメーションをCesium上で実現させるにはどうすればよいのか実装方法を模索していたところ、Cesium公式のSandcastleで点群の3Dモデルをアニメーションさせるデモを発見しました。こちらのデモでは点群の3Dモデルに対してカスタムシェーダーを適用することで高速かつ滑らかなアニメーションを実現しています。

こちらのデモを参考に気圧データを点群でアニメーションさせることを試してみました。
点群の3Dモデルを用意し、カスタムシェーダーにより各点を時系列の気圧データに従って高さを変化させてアニメーションさせます。
また、インタラクティブに任意の日時を指定できるようにするため、Cesiumの時系列スライダーで操作できるようにします。

点群3Dモデル作成

アニメーションさせる点群3Dモデルを作成します。

気圧の時系列データをGeotiffから読み込んで点群3Dモデルに反映させてアニメーションさせます。
この時、点群の頂点とGeotiffの各ピクセルを1対1で対応させて、対応するGeotiffのピクセルの値を点群の頂点の高さなどに反映させる実装を考えました。
これを実現するには、点群の各頂点からGeotiffの対応するピクセルを参照できるようにしなければなりません。そのため点群3DモデルにGeotiffの対応するピクセルの情報を持たせる必要があります。
幸いCesiumに表示可能な3Dモデル形式である「glTF」は属性情報を保持可能な3Dモデルの形式になりますので、Geotiffの属性情報をもつ点群glTFモデルを作成しました。

glTFモデルはGeotiffを元にPythonのスクリプトを組んで生成しました。
まずPythonライブラリのGDALでGeotiffを読み込み、座標情報を3Dモデル用の座標値に変換します。
3次元地理座標変換ライブラリであるpymap3dを使用して以下のようにWGS84座標からECEF座標へ変換を行いました。

import pymap3d as pm

def xyz_to_ecef(lon, lat, alt):
    # WGS84座標をECEF座標に変換
    x, y, z = pm.geodetic2ecef(lat, lon, alt)
    return x, y, z

glTFの読み書きを行うライブラリpygltfを使用してglTFモデルを作成します。3Dモデルの頂点の座標値として先ほど変換した座標値を設定し、さらに対応するGeotiffのピクセルのインデックス情報と高さ更新時の位置情報の計算で使用されるWGS84の座標値を属性情報として設定してglTFモデルを書き出します。

import pygltflib
...

gltf = pygltflib.GLTF2(

...

accessors=[
  pygltflib.Accessor(
    bufferView=0,
    componentType=pygltflib.FLOAT,
    count=len(points),                       # 3Dモデルの頂点の座標値(ECEF座標)
    type=pygltflib.VEC3,
    max=points.max(axis=0).tolist(),
    min=points.min(axis=0).tolist(),
  ),
  pygltflib.Accessor(
    bufferView=1,
    componentType=pygltflib.FLOAT,
    count=len(lonlat_list),                  # WGS84の座標値
    type=pygltflib.VEC2,
    max=lonlat_list.max(axis=0).tolist(),
    min=lonlat_list.min(axis=0).tolist(),
  ),
  pygltflib.Accessor(
    bufferView=2,
    componentType=pygltflib.UNSIGNED_INT,
    count=len(idx_list),                     # ピクセルのインデックス情報
    type=pygltflib.VEC2,
    max=idx_list.max(axis=0).tolist(),
    min=idx_list.min(axis=0).tolist(),
  ),
],

...

こうして点群glTFモデルを作成することができました。
GraphCastのGeotiffを元に作成した3Dモデルは全球を覆うような形状となっています。3Dモデルをアップで見ると等間隔で点が並んでいることが分かりますが、これらの点はGeotiffの各ピクセルと同じ座標値に設置されています。
また、属性情報として生成元のGeotiffの対応するピクセルのインデックスを保持させることができました。

カスタムシェーダー

カスタムシェーダーで気圧の値によって点群モデルの頂点の高度が変わるように実装します。
Cesium Sandcastleの点群モデルアニメーションのコードを元に実装しました。

点群の頂点の位置の操作についてはvertexShaderで行います。
以下が実際のコードで、点群モデルの属性情報からGeotiffのインデックス情報を取得して気圧の値を参照し、点群の頂点に新しい高さを設定しています。

　　　void vertexMain(VertexInput vsInput, inout czm_modelVertexOutput vsOutput)
      {
        vec2 lonlat = vsInput.attributes.longlat; /* モデルの属性情報から座標値と
Geotiffのインデックスを取得 */
        vec2 geotiff_idx = vsInput.attributes.idx;
        vec2 uv_pos = vec2(float(geotiff_idx.x) / 1439.0, float(geotiff_idx.y) / 719.0);
        float zval = get_tex_value(uv_pos); /* Geotiffデータからインデックスが参照する気圧の値を取得 */
        float newz = 1000.0 + ((zval - zmin) * heightBuf);
        vec3 newXyz = geodeticToECEF(lonlat.x, lonlat.y, newz); /* 気圧の値から高さを計算し、頂点の新しい座標値を生成 */
        vsOutput.positionMC.x = newXyz.x; /* 頂点に新しい座標値を設定 */
        vsOutput.positionMC.y = newXyz.y;
        vsOutput.positionMC.z = newXyz.z;
        vsOutput.pointSize = 3.0;
      }

点群の高さに加えて色も反映されるようにfragmentShaderを作成します。
取得した気圧値によってCS立体図やカラーマップなどの色を計算して点群の各点に反映させます。
CS立体図等の計算はfrogcat様の記事を参考にさせていただきました。

vec3 get_cs_map(vec2 uv_pos) {
  mat3 h;
  h[0][0] = get_tex_value(uv_pos + (vec2(-1,-1) * unit));
　...
  return cs_map * dem_color * contour;
}

void fragmentMain(FragmentInput fsInput, inout czm_modelMaterial material)
{
  vec2 geotiff_idx = fsInput.attributes.idx; /* モデルの属性情報からGeotiffのインデックスを取得 */
  vec2 uv_pos = vec2(float(geotiff_idx.x) / 1439.0, float(geotiff_idx.y) / 719.0);
  material.diffuse = get_cs_map(uv_pos); /* Geotiffデータから気圧値を取得、CS立体図等を計算して点群の色に反映 */
}

時系列データ読み込み

時系列データは以下の形式のGeotiffから時系列で読み込もうとしました。

• 画像サイズ：1440x721（1ピクセル0.25度の全球データ）
• バンド数：40（10日間6時間ごとの時点）
• 浮動小数点型の気圧データを格納

JavascriptライブラリのGeotiff.jsを使用してブラウザから40バンドの時系列データのGeotiffをいっぺんに読み込もうとしたところ、メモリエラーとなってしまいました。

そこから試行錯誤した結果、40バンドのGeotiffを1バンドのGeotiff40枚に分離して、1枚ずつ順番に40枚読み込む場合はメモリエラーにはならないことが分かりました。
そのため最終的に、最初に40枚すべての時系列データを読み込む形式で実装しました。
時系列スライダーを動かすと、読み込み済みの時系列データから現在の日時のデータを参照します。

3D点群アニメーション

こうして時系列気圧データの点群アニメーションができました。
時系列によって変化する気圧にしたがって点群の各頂点の高さが変化し立体的なアニメーションになっています。

こちらの問題点として、点群であるためアップでみた場合に点の間が空いてしまっているため、色や起伏が分かりにくい点があります。

そこで、点群のモデルに対して頂点を動かしていましたが、頂点間が面で埋められているメッシュモデルに対しても適用することが可能なのか試してみました。

メッシュ3Dモデル作成

アニメーションさせるメッシュ3Dモデルを作成していきます。
点群モデルを作成するPythonスクリプトを改修し、open3dというライブラリで点群をメッシュ化してglTFモデルを作成します。先ほどの点群glTFモデルにメッシュの情報を追加します。

# メッシュ作成
# 点の法線推定
ptCloud = o3d.geometry.PointCloud()
ptCloud.points = o3d.utility.Vector3dVector(points_ecef_mesh)
ptCloud.estimate_normals()
ptCloud.orient_normals_consistent_tangent_plane(100)
# メッシュ再構成
distances = ptCloud.compute_nearest_neighbor_distance()
avg_dist = np.mean(distances)
radius = 2*avg_dist   
radii = [radius, radius * 2]
recMeshBPA = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(
        ptCloud, o3d.utility.DoubleVector(radii))
triangles = np.array(recMeshBPA.triangles, dtype=np.uint32) # 三角ポリゴンの情報をglTFの設定に追加する

メッシュの生成がうまくいかず穴が空いている箇所もありますが、全球を覆うようなメッシュ3Dモデルが作成できました。

3Dメッシュアニメーション

先ほどの点群モデルアニメーションをメッシュモデルに差し替えてアニメーションさせてみます。カスタムシェーダーは変更せず、点群モデルに適用したものと同じものを使用しました。

点群モデルを使用した場合と比較してアップで見た場合でも起伏の変化が分かりやすくなったかと思います。

去年のアドベントカレンダーで行った複数枚の画像レイヤーを使用したアニメーションとも比較すると、何十枚と画像レイヤーを重ねていたため動作が重かったのですが、3Dメッシュモデルはそれよりも挙動が軽い結果となりました。
また、リアルタイムでパラメータを変更して、色味の調節や等高線の描画を反映することも可能です。

他のデータも可視化してみる

FOSS4Gでの発表ではGraphCastで予測した全球の気圧データを可視化した結果をお見せしていましたが、他にも2023年8月の台風7号の時系列気圧、雨量データも可視化を試していました。（元データは京都大学の生存圏データベースで公開している気象庁データから取得）

2023年8月 台風7号の気圧データ

台風の気圧データでは、上陸前までは台風の部分が穴のようになっており非常に気圧が低いことが分かりやすく可視化できています。

横から見ると台風の気圧の変化がより分かりやすいです。

2023年8月 台風7号の雨量データ

一部分が棘のように変化していくアニメーションとなっており、局所的に雨量が増加していることが分かります。

まとめ

気圧データを時系列3Dアニメーションさせてみました。
これにより専門知識が無い人でも気圧の起伏の変化をわかりやすく表現することができました。
今回は気圧という気象データを取り上げましたが、こちらの手法は気象以外でもグリッド形式のデータであれば適用できるため、気象に限らず専門家ではない人にもわかりやすいアニメーション可視化が期待できます。

また、今までは全球クラスの時系列データでは2Dで表示することが多く3D表示では動作が重くなってしまう印象がありました。
しかしCesiumのカスタムシェーダーを使用することで従来よりも軽い動作で3D可視化を実現することができました。
実際に今回のアニメーションはGPUがないスマートフォン上でも動作することを確認しております。

今回作成した時系列データの3DアニメーションのソースコードはGithubで公開していますので、興味がある方はぜひ試してみてください。

最後に今回の実装に当たり、各種オープンソースのGISツールを使用することで比較的簡単に実装することができました。
OSSツールの開発者の皆様ありがとうございました。