2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@haystacker(大塚 隼斗)
in

オフライン GIS 高さ計算システムの実装

2
Posted at

GIS や測量、3D 地理空間処理では、「高さ」をどの基準面で扱うかが非常に重要です。
特に、正高 $H$（標高） と 楕円体高 $h$ を混同すると、地形解析、3D モデル配置、センサ位置推定、経路計算などで無視できない誤差につながります。

本記事では、DTM（GeoTIFF） と ジオイドモデル（EGM2008 / JPGEO2024） を統合し、正高 $H$ ⇔ 楕円体高 $h$ をオフラインで相互変換する方法について解説します。
さらに、グローバル版の EGM2008 と 日本域向けの JPGEO2024 を比較し、日本国内でどちらを採用すべきかも精度評価を通じて確認します。

この記事で分かること

  • 正高 $H$・楕円体高 $h$・ジオイド高 $N$ の関係
  • GeoTIFF + ジオイド格子 を使ったオフライン高さ変換の実装方法
  • EGM2008 と JPGEO2024 の違い と、日本域での実務上の選び方
  • 高さ変換システムにおける I/O 最適化・補間・例外処理 の考え方

対象読者

本記事は、以下のような読者を想定しています。

  • Python で地理空間データ処理を行いたい人
  • GeoTIFF、CRS、座標変換の基礎をある程度理解している人
  • 標高、楕円体高、ジオイド高の違いを実装レベルで整理したい人
  • API に依存せず、現場・災害時・オフライン環境でも動作する高さ変換を実装したい人

目次

  1. 背景
  2. 事前準備
  3. システムアーキテクチャ
  4. 主要な実装の解説
  5. 実行結果と精度評価
  6. 実装設計の考え方
  7. まとめ
  8. ソースコード
  9. 参考文献・関連資料

1. 背景

1.1. 高さの基準

GIS や測量で「高さ」を扱う場合、まず重要になるのは、どこを高さの基準とするかです。
現実の地球表面には凹凸があり、さらに地球は完全な球体でもありません。そのため、高さは必ず何らかの基準面からの距離として定義されます。

実務では主に次の 2 種類の高さが使われます。

  • 正高 (Orthometric Height) $H$
    平均海水面を基準にした高さ。一般に「標高」と呼ばれる
  • 楕円体高 (Ellipsoidal Height) $h$
    地球を数学的な回転楕円体で近似した地球楕円体を基準にした高さ

これら 2 つの高さを結び付ける量が、ジオイド高 (Geoid Height) $N$ です。

ここでいう ジオイド (Geoid) とは、地球の重力ポテンシャルが等しい面 (等ポテンシャル面)　のうち、平均海水面に最もよく一致する起伏を持つ基準面を指します。
地球内部の質量分布は場所ごとに異なるため、ジオイドは完全な球面や楕円体ではなく、場所ごとに上下する起伏を持つ面になります。

そのため、ジオイド高 $N$ も地域によって変化します。
全球ジオイドモデル（例：EGM2008）や地域ジオイドモデル（例：JPGEO2024）で与えられる $N$ の精度が、そのまま高さ変換の精度に影響します。

これらの関係は次式で表されます。

h = H + N

1.2. なぜ「オフライン」でやるのか

高さ計算をオンライン API に依存すると、次のような制約が発生します。

  • 通信断・帯域制約
    （現場、災害時、海上、山間部など）
  • レイテンシ
    （リアルタイム処理・大量点処理でのレスポンス低下）
  • 利用制限
    （アクセス制限、仕様変更、サービス停止リスク）

そのため、DTM とジオイドモデルをローカルに保持し、オフラインで完結する高さ計算は、実務・研究の両面で非常に強い構成になります。

特に、次のような用途ではオフライン化の価値が高くなります。

  • 現地運用を前提とする GIS / モバイル端末
  • 災害対応やインフラ点検
  • ドローン航法や 3D 地理空間処理
  • API 非依存が求められる長期運用システム

1.3. 実装上の課題

オフラインで精度と速度を両立するには、主に次の点が課題になります。

  1. 大容量データの効率的読み込み
  2. 格子点からの補間
  3. 座標系変換（CRS）
  4. nodata、範囲外、端ピクセルなどの例外処理
  5. DTM 側の高さ基準とジオイド側の基準面の整合

2. 事前準備

2.1. DTM ファイルの準備

  1. 国土地理院の DEM（JPGIS(GML) 型式）をダウンロードする

  2. エコリス社の「基盤地図情報 標高DEMデータ変換ツール」などを用いて GeoTIFF に変換する

    • QGIS の QuickDEM4JP プラグイン でも同様の処理は可能だが、、DEM データサイズが大きい場合にエラーが出るようなので、今回はエコリス社のツールを使用した

2.2. ジオイドモデルの準備

  1. グローバル版のジオイドモデル（egm2008-1.pgm）をダウンロードする

  2. 日本版のジオイドモデル（JPGEO2024.isg）をダウンロードする

2.3. 今回の前提

本記事では比較検証をしやすくするため、次の前提で実装します。

  • 入力座標は WGS84（EPSG:4326）の緯度経度

  • DTM は GeoTIFF

  • ジオイドモデルは以下の 2 系統を使用

    • EGM2008：全球モデル、PGM 形式
    • JPGEO2024：日本域モデル、ISG 形式

  • DTM 側の高さ基準は、明示的に vertical_datum として扱う

この「基準をコード上で明示する」ことが、実務ではかなり重要です。

3. システムアーキテクチャ

今回は比較検証のため、グローバル版 と 日本版 の 2 系統で実装します。

  • グローバル版（main_global.py）
    EGM2008（PGM）を使用
  • 日本版（main_japan.py）
    JPGEO2024（ISG）を使用

両者の違いは、基本的にはジオイド読み込み部のみです。
それ以外の流れは共通化できます。

処理フローは次の通りです。

  1. DTM（GeoTIFF）処理
    座標変換 → 近傍 2×2 の取得 → バイリニア補間
  2. ジオイド格子処理
    格子読み込み → バイリニア補間で $N$ を算出
  3. 統合
    vertical_datum に応じて $H \leftrightarrow h$ を相互変換

image.png

3.1. この構成の重要ポイント

本構成の本質は、「地形の高さ」と「基準面の差」を分離して扱うことです。

  • DTM は地表面の高さ情報を与える
  • ジオイドモデルは基準面差 $N$ を与える
  • 最後に両者を統合して $H$ または $h$ を得る

この分離により、誤差要因を次のように切り分けやすくなります。

  • DTM 由来の誤差
  • ジオイドモデル由来の誤差
  • 補間や例外処理の実装由来の誤差

つまり、実装の正しさだけでなく、改善すべき箇所の特定にも有利です。

4. 主要な実装の解説

4.1. DTM（GeoTIFF）からの高さ取得

(1) WGS84 → DTM の CRS へ変換

入力は緯度経度（EPSG:4326）であることが多いため、pyproj.Transformer を用いて DTM の CRS に変換します。
always_xy=True を指定することで、引数順を lon, lat に固定します。

transformer = Transformer.from_crs(CRS.from_epsg(4326), ds.crs, always_xy=True)
x, y = transformer.transform(lon, lat)

ここは小さく見えて、実務では非常に重要です。
lat, lon と lon, lat の取り違えは、高さ取得失敗の典型原因のひとつです。

(2) (x, y) を浮動小数のピクセル座標へ変換し、2×2 を読む

(~ds.transform) * (x, y) によって、浮動小数のピクセル座標 (col_f, row_f) を得ます。

col_f, row_f = (~ds.transform) * (x, y)

その後、

  • c0 = floor(col_f)
  • r0 = floor(row_f)

を左上セルとして、2×2 の近傍を読み込みます。
端部では範囲外に出ないよう、必要に応じて内側にクランプします。

window=((r0, r0+2), (c0, c0+2)) の +2 は、Python / rasterio の半開区間 [start, stop) により、2 要素を取得するためです。

(3) 補間係数 fx, fy は小数部から作る

補間係数は、浮動小数ピクセル座標の小数部を使います。

fx = col_f - c0
fy = row_f - r0

必要に応じて 0〜1 にクリップし、2×2 の格子点に対してバイリニア補間を行います。

4.2. ジオイドの読み込みと N の算出

4.2.1. EGM2008（PGM：memmap）

egm2008-1.pgm は、EGM2008 のジオイド高 $N$ を格子状に格納したファイルです。
見た目は PGM（Portable GrayMap）ですが、ここで扱うのは一般的な画像用途の PGM というより、ジオイド格子を格納するための PGM ベースの専用フォーマットと考えた方がよいです。

P5
# Geoid file in PGM format for the GeographicLib::Geoid class
# Description WGS84 EGM2008, 1-minute grid
# URL http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
# DateTime 2009-08-31 06:54:00
# MaxBilinearError 0.025
# RMSBilinearError 0.001
# MaxCubicError 0.003
# RMSCubicError 0.001
# Offset -108
# Scale 0.003
# Origin 90N 0E
# AREA_OR_POINT Point
# Vertical_Datum WGS84
21600  10801
65535

EGM2008（PGM）の読み取りポイント

項目 値の例 説明 備考
解像度 1-minute grid 1分格子 1 minute = 1/60°
グリッド数 21600 10801 横 21600 × 縦 10801 width, height
経度範囲 360° 分 全球 21600 = 360 × 60
緯度範囲 北極〜南極 全球 両極を含むため 10800 + 1
データ型 65535 16bit unsigned integer ビッグエンディアン
実値変換 Offset, Scale 格納値から実際のジオイド高へ変換 重要

この形式では、格納値そのものがメートル値ではないため、Offset と Scale を用いて実際のジオイド高へ復元する必要があります。

4.2.2. JPGEO2024（ISG：テキスト格子 + LRU キャッシュ）

JPGEO2024.isg は、日本域向けジオイドモデルを格子形式で保持したファイルです。
ISG 形式は、ジオイド格子交換用のフォーマットです。

begin_of_head ================================================
model name     : JPGEO2024
model year     : 2024
model type     : gravimetric
data type      : geoid
data units     : meters
data format    : grid
data ordering  : N-to-S, W-to-E
ref ellipsoid  : GRS80
ref frame      : ---
height datum   : ---
tide system    : ---
coord type     : geodetic
coord units    : dms
map projection : ---
EPSG code      : ---
lat min        =  15°00'00"
lat max        =  50°00'00"
lon min        = 120°00'00"
lon max        = 160°00'00"
delta lat      =   0°01'00"
delta lon      =   0°01'30"
nrows          =        2101
ncols          =        1601
nodata         =  -9999.0000
creation date  =  01/04/2025
ISG format     =         2.0
end_of_head ==================================================

JPGEO2024（ISG）の読み取りポイント

項目 値の例 説明 備考
解像度（緯度） 0°01'00" 緯度方向 1′ 約 1/60°
解像度（経度） 0°01'30" 経度方向 1′30″ 1.5′ = 1/40°
グリッド数 2101 × 1601 行列サイズ nrows, ncols
緯度範囲 15°N ～ 50°N 日本域をカバー 差 35°
経度範囲 120°E ～ 160°E 日本域をカバー 差 40°
データ順序 N-to-S, W-to-E 北→南、西→東 行優先
単位 meters メートル
欠損値 -9999.0000 無効値 nodata

この形式はテキストベースで読みやすい一方、毎回フルパースするとコストが大きくなります。そのため、ヘッダ解析結果や近傍格子の再利用にはキャッシュ戦略が有効です。

また、初見の読者向けに直感的に言えば、JPGEO2024 の ISG 形式は日本域を対象とした比較的高密度な地域ジオイド格子です。そのため、日本国内では全球モデルよりも整合しやすいことが期待できます。

4.3. バイリニア補間

DTM でもジオイド格子でも、基本的な補間方法は同じです。
4 点（左上・右上・左下・右下）を使い、格子内の位置 (fx, fy) に応じて値を求めます。

v0 = z00 * (1 - fx) + z10 * fx
v1 = z01 * (1 - fx) + z11 * fx
result = v0 * (1 - fy) + v1 * fy

数学的には、バイリニア補間は次式で表せます。

f(x,y)=a_{00}+a_{10}x+a_{01}y+a_{11}xy

ただし、実装上は「x 方向の線形補間を 2 回行い、その結果を y 方向にもう一度補間する」という 2 段階の方が分かりやすく、実装もしやすいです。

また、計算量は 1 点あたり $O(1)$ で一定です。
大量点処理でも扱いやすく、実務での速度と精度のバランスが良い方法です。

4.4. 統合

DTM が MSL 基準（正高） なのか、楕円体高基準 なのかによって、変換方向は逆になります。

  • DTM が MSL 基準のとき、DTM の高さを $H$ とみなし、

    h = H + N

  • DTM が楕円体高基準のとき、DTM の高さを $h$ とみなし、

    H = h - N

この分岐を曖昧にすると、実装は動いていても結果が間違う、という事態が起こります。
そのため、vertical_datum を明示的な型として持たせ、分岐を固定するのが実務的には安全です。

5. 実行結果と精度評価

本章では、精度評価を DTM 側の正高 $H$ と ジオイド高 $N$ に分離して整理します。
最終的な楕円体高は

h = H + N

で与えられるため、誤差要因を分離して見ることで、どこを改善すべきかを判断しやすくなります。

5.1. 評価設計

  • 対象地点：那覇 5 点 + 新宿 5 点（計 10 点）

  • 比較基準：

    • 正高 $H$：国土地理院 API による値
    • ジオイド高 $N$：比較用に採用した基準値

  • DTM：DEM1A から生成した GeoTIFF（日本版 / グローバル版で共通）

  • ジオイド：

    • 日本版：JPGEO2024
    • グローバル版：EGM2008

  • 評価対象：

    • $H$（正高, MSL） = DTM 側の影響
    • $N$（ジオイド高） = ジオイドモデル側の影響
    • $h$（楕円体高） = 統合結果

ここで重要なのは、DTM は両方式で共通であり、差が出るのは主にジオイドモデル側だという点です。

5.2. DTM 高さ（正高 H）の精度評価

ここでは、正高 $H$ の一致度を評価します。
差異は次のように定義します。

差異 = オフライン算出の正高 H - 比較基準の正高
地点 GSI API
正高 $H$		 GSI API
データモデル		 日本版
正高 $H$		 差異 グローバル版
正高 $H$		 差異
那覇_1 2.200 5m（写真測量） 2.820 +0.620 2.820 +0.620
那覇_2 87.300 5m（写真測量） 87.336 +0.036 87.336 +0.036
那覇_3 28.300 5m（写真測量） 28.892 +0.592 28.892 +0.592
那覇_4 28.300 5m（写真測量） 28.257 -0.043 28.257 -0.043
那覇_5 2.300 5m（写真測量） 2.995 +0.695 2.995 +0.695
新宿_1 33.200 5m（レーザ） 33.217 +0.017 33.217 +0.017
新宿_2 30.300 5m（レーザ） 30.101 -0.199 30.101 -0.199
新宿_3 27.000 5m（レーザ） 26.980 -0.020 26.980 -0.020
新宿_4 19.600 5m（レーザ） 20.394 +0.794 20.394 +0.794
新宿_5 6.300 5m（レーザ） 5.890 -0.410 5.890 -0.410

統計（10 地点）

  • MAE [m]：0.3425
  • RMSE [m]：0.4537
  • 最大誤差 [m]：0.7938
  • 最小誤差 [m]：0.0169
  • 標準偏差 [m]：0.2976

地域別

  • 那覇（5 地点）

    • MAE [m]：0.3971
    • RMSE [m]：0.4941

  • 新宿（5 地点）

    • MAE [m]：0.2879
    • RMSE [m]：0.4095

所見

  • 正高誤差は全地点でおおむね ±1 m 以内

つまり、今回の差の主要因は DTM 側よりも、むしろジオイドモデルの選定にあると考えられます。

5.3. ジオイド高 N の精度評価

ここでは、ジオイド高 $N$ のモデル差を評価します。
同一点で比較した場合、日本版（JPGEO2024）は比較基準と整合しやすく、グローバル版（EGM2008）は日本域において系統差を持つことが確認できます。

地点 API
ジオイド高 $N$		 日本版
ジオイド高 $N$		 差異 グローバル版
ジオイド高 $N$		 差異
那覇_1 30.8623 30.8623 0.0000 30.5404 -0.3219
那覇_2 30.7319 30.7319 0.0000 30.4245 -0.3074
那覇_3 30.7434 30.7434 0.0000 30.4185 -0.3249
那覇_4 30.7073 30.7073 0.0000 30.4028 -0.3045
那覇_5 31.0245 31.0245 0.0000 30.6951 -0.3294
新宿_1 36.9830 36.9830 0.0000 36.6048 -0.3782
新宿_2 37.3373 37.3373 0.0000 36.9489 -0.3884
新宿_3 37.1352 37.1352 0.0000 36.7480 -0.3872
新宿_4 37.1185 37.1185 0.0000 36.7392 -0.3793
新宿_5 37.0422 37.0422 0.0000 36.6660 -0.3762

統計（10 地点, 比較基準）

  • EGM2008 の $N$ に対する MAE [m]：0.3497
  • 平均差異（JPGEO2024 - EGM2008） [m]：0.3497
  • 最大差異 [m]：0.3884
  • 最小差異 [m]：0.3044
  • 標準偏差 [m]：0.0331

解釈

  • 日本版の $N$ は比較基準に整合
  • EGM2008 の差は、採用ジオイドモデルの差による系統差

日本域で基準面整合を重視するなら、JPGEO2024 を採用した方がよいと言えます。

5.4. 統合結果（h の評価）

統合結果である楕円体高は、次式で求められます。

h = H + N

このため、$H$ が同じでも、$N$ が異なれば $h$ は変わります。

DEM1A 条件での楕円体高 $h$ の統計は次の通りです。

  • 日本版（JPGEO2024）

    • MAE：0.3425 m
    • RMSE：0.4538 m

  • グローバル版（EGM2008）

    • MAE：0.4106 m
    • RMSE：0.4381 m

補足

  • 今回の比較では $H$ は両方式で同一
  • $h$ の差は主に $N$ の系統差（約 0.35 m）に起因する

つまり、楕円体高を厳密に扱いたい用途ほど、ジオイドモデルの選定が効くということです。

5.5. 感度整理

誤差分解は近似的に次のように整理できます。

\Delta h \approx \Delta H + \Delta N

ここで、

  • $\Delta H$：DTM ソース、補間仕様、再サンプリング、地形条件の影響
  • $\Delta N$：採用ジオイドモデル（JPGEO2024 / EGM2008）の基準面差

実務上の示唆

  1. 日本域で $H$ の整合を重視する場合、まず DTM 品質 を改善する
  2. $h$ を厳密に扱う場合、ジオイドモデル選定 がそのまま成果物の基準面を決める
  3. MAE と RMSE は分けて評価し、大誤差点の個別分析を行う

5.6. 評価の限界

今回の評価には、いくつかの限界もあります。

  • 評価点は 10 地点 であり、全国一律の一般化には十分ではない
  • 都市部、急傾斜地、海岸近傍などでは誤差傾向が変わる可能性がある
  • 比較基準側も、補間仕様や採用モデルに依存する

したがって、本結果は実装方針の妥当性確認と、日本域でのモデル差の傾向把握には有効ですが、全国網羅的な性能評価としては追加検証の余地があります。

6. 実装設計の考え方

6.1. I/O 最適化

本実装の主なボトルネックは、演算そのものよりもデータをどう読むかにあります。
そのため、I/O を減らす設計を優先しています。

DTM（GeoTIFF）

毎回 2×2 ウィンドウのみを読む構成にしています。

  • バイリニア補間に必要な 4 点だけ取得する
  • 1 点評価あたりの読み出し量が一定
  • 全体読み込みや広い近傍探索を避けられる

EGM2008（PGM）

np.memmap により、必要セルだけへアクセスします。

  • ファイル全体を RAM に展開しない
  • OS のページキャッシュを利用できる
  • 巨大格子でも初期ロード時間とメモリ使用量を抑えやすい

再利用キャッシュ

_GEOID_CACHE や @lru_cache を用いて、固定コストを削減します。

  • 同一ファイルの再読込を避ける
  • ヘッダ解析の再実行を防ぐ
  • 複数点連続処理で効果が大きい

I/O 最適化は、単発の 1 点評価を極端に速くするというより、大量点処理時の総処理時間を安定して下げるために効きます。

6.2. 計算最適化

計算側は、単純で定数時間な処理として固定しています。

  • 補間は 4 点に対するバイリニア補間なので $O(1)$
  • 高次補間を避け、速度と実務精度のバランスを優先
  • 座標変換オブジェクトは本番処理では再利用可能

特に高速化では、計算式そのものを微調整するより、

  • 不要な I/O を減らす
  • 使い回せるオブジェクトを再利用する
  • 例外処理で無駄な再試行を防ぐ

といった設計の方が効きやすいです。

7. まとめ

本記事では、DTM とジオイドモデルを組み合わせ、正高 $H$ と楕円体高 $h$ をオフラインで相互変換する方法を解説しました。

ポイントは次の 3 つです。

  1. 高さ変換の本質は $h = H + N$ であること
  2. 誤差は DTM 側とジオイド側に分けて考えるべきこと
  3. 日本域で基準面整合を重視するなら JPGEO2024 が有利であること

精度評価では、DTM（DEM1A）による正高誤差はおおむね ±1 m 以内に収まりました。
一方で、ジオイドモデルの違いはそのまま楕円体高に反映され、日本域では JPGEO2024 が EGM2008 より約 0.35 m 良好に整合することが確認できました。

つまり、日本国内の実務では、単に「ジオイドを使う」だけでなく、どのジオイドモデルを採用するかが成果物の品質を左右します。

また、実装面では 2×2 ウィンドウ読み込みと np.memmap を用いて I/O を最小化しており、オフライン環境でも大量点処理に耐えられる構成となっています。

8. ソースコード

実装の詳細は、GitHub リポジトリを参照してください。

9. 参考文献・関連資料

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?