グリッド補間の実装 — バイリニア内挿の中身

⚠ 本記事は、GeoPrism JPの「測地系を変換する2つの流儀 — パラメータ変換とグリッド変換をやさしく学ぶ」で紹介した「グリッド変換」という考え方の実装特化版です。前回は「格子点に補正値を置いて、間を補間する」という概念をやさしく紹介しましたが、今回は「間を補間する」の中身、つまりバイリニア補間の数式と実装上の注意点を掘り下げます。

TKY2JGD・PatchJGD・SemiDynaEXE・POS2JGDといった国土地理院方式の変換や、ジオイドモデルの内挿も、根っこの計算方法は同じバイリニア補間です。

おさらい：グリッド補正はなぜ「格子＋補間」なのか

TKY2JGDやPatchJGDのような地殻変動補正は、実際の地殻変動という非線形で複雑な現象を扱います。これを1本の数式で表すことはできないため、国土地理院は緯度経度の格子（メッシュ）の各交点に「その地点での補正量」をあらかじめ計算して格納したパラメータファイル（.par）を配布しています。

利用者が変換したい地点は、格子の交点ちょうどに乗ることはほとんどありません。そこで必要になるのが、周囲4点の補正量から、格子内の任意の1点の補正量を求める補間計算——これがバイリニア（双線形）補間です。

バイリニア補間の数式

対象地点が、緯度経度の格子の4つの角 $(x_0, y_0)$・$(x_1, y_0)$・$(x_0, y_1)$・$(x_1, y_1)$ に囲まれているとします。各角の補正値を $f_{00}, f_{10}, f_{01}, f_{11}$ とし、対象地点が格子内のどの位置にあるかを、辺の長さで正規化した2つの比率 $u, v$（それぞれ0〜1の範囲）で表すと、補間値 $f(u, v)$ は次の式で求まります。

$$
f(u, v) = f_{00}(1-u)(1-v) + f_{10}\,u(1-v) + f_{01}(1-u)v + f_{11}\,uv
$$

これは「経度方向に1回線形補間（2回）→ その結果を緯度方向にもう1回線形補間」と分解しても同じ結果になります（先に緯度方向→経度方向の順でも同様）。2方向の線形補間を組み合わせているため「バイ（2重）リニア」と呼ばれます。

TKY2JGD補正量の全国分布ヒートマップ

上の図はTKY2JGDの補正量を全国メッシュで可視化したものですが、実際にアプリが1点を変換するときは、この格子データから該当する1マスを特定し、その4隅の値をバイリニア補間しています。

実装で気をつける3つのポイント

1. メッシュコードから4隅を求める

.parファイルはメッシュコード（例：3次メッシュなら南北30秒×東西45秒）で区切られた格子点の値を格納しています。対象の緯度経度から、まず「どのメッシュに属するか」を計算し、そのメッシュの南西・南東・北西・北東の4隅に対応するレコードをパラメータファイルから引く必要があります。緯度経度をメッシュ番号に変換する式自体は単純な除算・切り捨てですが、南北半球や経度180度をまたぐような特殊ケースは考慮不要な国内専用データである分、実装は比較的シンプルです。

2. 正規化座標 $(u, v)$ への変換

対象地点の緯度経度から、そのメッシュ内での相対位置を0〜1の比率に変換します。単純に「メッシュ内での位置 ÷ メッシュの辺の長さ」ですが、緯度と経度でメッシュの辺の長さの単位（秒）が異なる点に注意が必要です（3次メッシュは南北30秒・東西45秒で、正方形ではありません）。

3. 欠損格子点のフォールバック

補正パラメータが定義されていない海域や、補正エリアの境界付近では、4隅のうち一部が欠損しているケースがあります。この場合にどう振る舞うか（4隅のうち有効な点だけで代用するか、変換不可としてエラーを返すか）は、精度の信頼性に直結する設計判断です。CCQ（座標変換クオリティー）という考え方で、変換結果がどれだけ信頼できるかを利用者に示す仕組みも、この欠損への対応と表裏一体です。

精度への影響

バイリニア補間の誤差は、格子の間隔（メッシュサイズ）が細かいほど、また補正量の変化が滑らかであるほど小さくなります。国土地理院のグリッドデータは、地殻変動という比較的滑らかに変化する現象を対象としているため、標準的なメッシュサイズであれば補間誤差は実務上問題にならない水準に収まるよう設計されています。実際の変換精度がどの程度国土地理院の公式値と一致するかは、全テストポイントでの比較検証で確認済みです。