Upgrade to Pro — share decks privately, control downloads, hide ads and more …

[FOSS4G 2024 Japan] CesiumのShaderを用いた気象データ3Dアニメ...

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

★ このスライドは、FOSS4G 2024 Japan の発表資料です ★

気圧の谷って、なんで谷と呼ばれるのか気になりませんか?
気圧線を見ても直感的に谷とかよくわからないので、Cesiumを使って3Dで可視化できるかを試してみました

発表者: 小林 知愛(株式会社ノーザンシステムサービス 研究開発部)

More Decks by ノーザンシステムサービス | Northern System Services

Other Decks in Technology

Transcript

  1. FOSS4G Advent Calender 2023 19日目 GraphCastの気象予測データをCesiumで可視 化してみた1 • AIが予測した気象データをCesiumで 可視化してみた記事

    • 今回の発表はこちらの記事のアペンド 的内容 • 最初に記事の内容を簡単に紹介 [1] https://qiita.com/wayama_ryousuke/items/f1823be92153e9b2928b 3
  2. 背景 • 近年、異常気象による災害の被害が増加 早期に異常気象を発見するために 長期間を正確に気象予測する技術が重要 [1]気候変動監視レポート 2023 | 気象庁 https://www.data.jma.go.jp/cpdinfo/monitor/2023/pdf/ccmr2023_all.pdf

    [2] 2023 年(令和 5 年)の世界の主な異常気象・気象災害(速報) | 気象庁 https://www.jma.go.jp/jma/press/2312/22d/2023matome_besshi2-2.pdf 日本での気温の上昇、降水量の増加1 世界的に異常気象が発生2 4
  3. 背景 • Google DeepMind 気象予測AI モデル「GraphCast」1を公開 • 従来のシステムよりも高精度かつ高速に最大 10 日間先の気象を予測

    • 予測が難しいとされる異常気象の予測にも強い [1] GraphCast: AI model for faster and more accurate global weather forecasting https://deepmind.google/discover/blog/graphcast-ai-model-for-faster-and-more-accurate- global-weather-forecasting/ 5
  4. GraphCast による気象予測 • ヨーロッパ中期予報センター(ECMWF)から公開されている GraphCast での推論ツール1で気象予測 2023/12/7 の気象データを入力して GraphCast で

    予測された 2023/12/9 の気圧のデータ 2023/12/9 に実際に観測された気圧データ 7 [1] ecmwf-lab/ai-models https://github.com/ecmwf-lab/ai-models
  5. CS 立体図を用いた気圧可視化アニメーション • 気圧の可視化 • 一般的に等高線で表されることが多い • 等高線だと専門知識がないとぱっと見ではわ かりにくい •

    「気圧の谷」っていうけど実際谷になってるの? • 地形判読に使用されるCS立体図で気圧を 可視化してアニメーション • 気圧の高低差が分かりやすく表現できるので は? 等高線の例1 12 [1] 日本周辺域 実況天気図の説明 | 気象庁 https://www.jma.go.jp/jma/kishou/know/kurashi/sokuhou_kaisetu.html
  6. CS 立体図を用いた気圧可視化アニメーション • 「FOSS4G japan 2020 online1」で発表したGeotiff 画像を地図タイ ル化し、Cesium 上にWebGL

    で CS 立体図を描画させる手法2 • 実装に当たりtmizu23様3、frogcat様4の記事を参考 [1] FOSS4G japan 2020 online https://www.osgeo.jp/events/foss4g-2020/foss4g-2020-japan-online [2] Geotiff.js で始めるリアルタイム演算 in foss4g japan 2020 online (https://www.slideshare.net/makinux7/geotiffjs-in-foss4g-japan-2020-online) [3] GeoTIFFからUTFGridを作る方法 -自然環境保全のための周辺技術https://tmizu23.hatenablog.com/entry/20121107 [4] 標高PNGタイルと WebGL による地形表現 - Qiita https://qiita.com/frogcat/items/7e91d3070a7a8d3e2c94 GraphCastが推論した気圧のCS立体図 CS立体図+カラーマップ 13
  7. おまけ:気圧を 3D で可視化してみる • 地形タイルにして3Dで立体的に可視化してみた • Cesium terrain builder1を使用してGeotiff 画像から地形タイルを作成

    • 起伏が小さいので高さ方向のスケールを 100 倍程度に拡大 2023 年 8 月の台風 7 号の気圧を地形タイル化 ここまでが去年の記事の内容 15 [1] geo-data/cesium-terrain-builder https://github.com/geo-data/cesium-terrain-builder
  8. 3D点群アニメーション 実装項目 • 点群3Dモデル作成 • Cesium上に表示する点群の3Dモデルの作成 • カスタムシェーダー • カスタムシェーダーによって点群3Dモデルをアニメーション

    • 時系列データ読み込み • アニメーションに反映させる時系列気圧データの読み込み • 時系列スライダーで任意の日時のデータを読み込むようにする 19
  9. 3D点群アニメーション 点群3Dモデル作成 • PythonでGeotiffから点群glTFモデルを生成 • PythonライブラリのGDALでGeotiffを読み込み • モデルの頂点用にGeotiffの座標値を3次元空間の座標値に変換 (WGS84座標→ECEF座標) •

    pymap3d1:3次元地理座標変換ライブラリ 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 22 [1] geospace-code/pymap3d https://github.com/geospace-code/pymap3d
  10. 3D点群アニメーション 点群3Dモデル作成 • PythonでGeotiffから点群glTF モデルを生成 • glTFモデルを生成 • pygltflib1:PythonでglTFの読み書 きを行うライブラリ

    • 対応するピクセルのインデックス とWGS84の座標値を属性情報とし て持たせる gltf = pygltflib.GLTF2( ... accessors=[ pygltflib.Accessor( bufferView=0, componentType=pygltflib.FLOAT, count=len(points), 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), 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(), ), ], ... 3Dモデルの 頂点の座標値 (ECEF座標) WGS84の座標値 ピクセルの インデックス 23 [1] pygltflib https://gitlab.com/dodgyville/pygltflib
  11. 3D点群アニメーション カスタムシェーダー • カスタムシェーダーで気圧の値が点の高度に反映されるように 実装 • vertexShaderで点群の頂点の位置を操作 void vertexMain(VertexInput vsInput,

    inout czm_modelVertexOutput vsOutput) { vec2 lonlat = vsInput.attributes.longlat; 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); 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; } モデルの属性情報から座標値と Geotiffのインデックスを取得 Geotiffデータからインデックスが参照する気圧の値を取得 気圧の値から高さを計算し、 頂点の新しい座標値を生成 頂点に新しい座標値を設定 25
  12. 3D点群アニメーション カスタムシェーダー • CS立体図やカラーマップ等が点群の色に反映されるように実装 • CS立体図や等高線等の計算にはfrogcat様1の記事を参考 • fragmentShader 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; vec2 uv_pos = vec2(float(geotiff_idx.x) / 1439.0, float(geotiff_idx.y) / 719.0); material.diffuse = get_cs_map(uv_pos); } モデルの属性情報からGeotiffのインデックスを取得 Geotiffデータから気圧値を取得、CS立体図等を計算して点群の色に反映 [1] 標高PNGタイルと WebGL による地形表現 - Qiita https://qiita.com/frogcat/items/7e91d3070a7a8d3e2c94 26
  13. 3D点群アニメーション 時系列データ読み込み • Geotiffから気圧データを時系列で読み込む • Geotiffの形式 • 画像サイズ:1440x721(1ピクセル0.25度の全球データ) • バンド数:40(10日間6時間ごとの時点)

    • 浮動小数点型の気圧データを格納 • Geotiff.js1を使用してブラウザ上から直接Geotiffを読み込む • GraphCastの気象予測データを格納したGeotiffは40バンドの時系列 データを持つ • 40バンドのGeotiffをいっぺんに読み込むとメモリエラーになる 27 [1] geotiff.js https://geotiffjs.github.io/
  14. 3D点群アニメーション 時系列データ読み込み • 1バンドのGeotiffを40枚読み込むとメモリエラーにはならない • 40バンドのGeotiffを1バンドのGeotiff x 40枚に分離して1枚ずつ順番 に読み込み •

    最初に全部の時系列の気圧データを読み込む形式で実装 • 時系列スライダーが動いた場合、読み込まれた時系列気圧データから 現在の日時のデータを参照 28
  15. メッシュアニメーション メッシュ3Dモデル作成 • Geotiffからメッシュモデルを作成 • 点群モデルを作成するPythonスクリプトを改修 • Pythonライブラリopen3d1で点群をメッシュ化して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の設定に追加する 32 [1] Open3D https://www.open3d.org/