経路断面の内挿

なぜ内挿が要るのか

経緯度が 2 次元の格子では、格子の行・列は緯線・経線と一致しません。たとえば ClimCORE のランベルト格子 (817 × 661、5 km) では、1 本の行の中で緯度が 30.3〜35.2°N と約 5° 変わり、1 本の列の中で経度が 10° 近く変わります。「格子の行に沿う」断面は内挿なしで値が正確ですが、「35°N の断面」の代わりにはなりません。

緯線や任意の 2 点間に沿った断面を描くには、経路の上の各点について「その点が格子のどこにあるか」を求め、周りの格子点の値から内挿する必要があります。経緯度が 1 次元の格子なら緯度・経度の軸ごとに 1 次元で位置が求まりますが、2 次元座標の格子では格子線が曲がっているので、経緯度から格子番号への逆算が要ります。これがこの章の中心です。

全体の流れ — 3 段階

内挿は「経路の点 → 小数の格子番号 → 値」の 3 段階で、どの段階も地図投影の知識を使いません。格子の経度・緯度の値 (と変数の値) だけから計算するので、投影パラメータを持たないファイルでも、どの投影の格子でも同じように動きます。

小数の格子番号とは、格子点 (行 j, 列 i) の番号を連続量に広げたものです。図の P の (fj, fi) = (26.27, 27.53) なら、行 26〜27・列 27〜28 のマスの中で、行方向に 0.27、列方向に 0.53 の位置にあることを表します (整数部分がマス、小数部分がマスの中の位置)。

1. 経路の上に点を取る

• 等緯度線 (緯度を固定): 指定した経度の始点・終点の間を等分し、緯度は固定値です。横軸は経度になります。
• 等経度線 (経度を固定): 緯度の始点・終点の間を等分し、経度は固定値です。横軸は緯度です。
• 2 点間の大円: 始点・終点を単位球上のベクトル p1, p2 にして、球面線形補間 (slerp) で等間隔の点を作ります。横軸は始点からの距離 (km) です。

• t を 0〜1 で等分するので、点は大円の上で等間隔になります。ω を atan2 で求めるのは、ごく近い 2 点でも大きく離れた 2 点でも精度が落ちないためです (arccos は ω が 0 や π の近くで精度が落ちます)。
• 経度は経路に沿って連続にします。日付変更線をまたいでも 179 → −179 と跳ばず 170 → 190 のように続くので、地図に経路を引くときや経度を横軸にするときに跳びが出ません。
• 始点と終点が地球の反対側 (対蹠点) だと大円が一つに決まらないので、エラーになります。
• 点の数は、既定では「経路の長さ ÷ 格子間隔」を切り上げて 1 を足した数 (2〜5000) です。格子間隔は、2 次元座標の格子では隣り合う格子点の球面距離の中央値、1 次元の格子では緯度間隔と「経度間隔 × cos(経路の緯度)」の小さい方です。ClimCORE (5 km) で 1000 km の経路なら約 200 点になります。「点の数を自動にする」を外せば数値で指定できます。

2. 点の位置を格子番号で表す

経緯度が 2 次元の格子

格子の経度・緯度の 2 次元配列だけを使って、各点の小数の格子番号を求めます。手順は 3 段で、図は点 P (真の格子番号は投影から作った (23.37, 31.62)) の例です。

手順 1 — 間引いた格子で最寄りを粗く探す (a)。経度・緯度を単位球上の 3 次元ベクトルにすると、2 点の球面上の距離は内積が大きいほど短くなります。そこで、格子を行・列とも間引いた点 (間引き幅 = 格子の大きい方の辺の点数 ÷ 100 を切り上げ) と P の内積をとり、最大のものを粗い最寄りとします。ClimCORE (817 × 661) では間引き幅 9、91 × 74 = 6,734 点との内積で済みます。3 次元ベクトルで比べるので、経度の規約 (0〜360 / −180〜180)・日付変更線・極の影響を受けません。

手順 2 — 窓を動かして本当の最寄りの格子点へ (b)。粗い最寄りを中心に、間引き幅の範囲の窓 (中心 ± 間引き幅) の全格子点と内積をとります。窓の中の最寄りが中心と違えばそこへ窓を動かしてやり直し、中心が窓の中の最寄りになったら確定します。図の例では粗い最寄り (25, 30) から最寄りの格子点 (23, 32) へ 1 回動いて止まります。

手順 3 — P で接する平面の上でニュートン法 (c・d)。最寄りの格子点の周りの格子点を、P で球に接する平面へ心射図法で写します。単位ベクトル v の写る先は次のとおりで、p は P の単位ベクトル、e・n は P での東向き・北向きの単位ベクトルです。この平面では P が原点に来ます。

格子番号の空間で正方形のマス (c) は、平面の上では四辺形になります (d)。マスの 4 つの角の平面座標 x00, x01, x10, x11 (y も同様) から、マスの中の位置 (a, b) (行方向 a、列方向 b、どちらも 0〜1) を平面の上の点へ写す双一次の写像を考え、X = Y = 0 (= P の位置) となる (a, b) をニュートン法で解きます。

• 2 × 2 のヤコビ行列 ∂(X, Y)/∂(a, b) は角の座標から解析的に書けます。最寄りの格子点 (a = b = 0) から始め、1 回の更新は各方向 1 マスまでに制限し、更新のたびに今いるマスを選び直すので、P が隣のマスにあっても移っていけます。最大 8 回、更新量が 1e-10 未満になったら終えます。解は格子番号 (fj, fi) = (マスの行 + a, マスの列 + b) です。
• 接平面で解く理由: 経度・緯度の平面で解くと、日付変更線での経度の跳び・極での経度の特異性・高緯度での歪みに個別の対処が要ります。接平面なら P の近くは歪みが小さく、東向き・北向きの単位ベクトルは極でも定義できます。心射図法では大円が直線になるので、マスの辺も平面の上でほぼ直線になります。
• 格子の外・収束しない点: 解いた後の残差 (平面の上での P とのずれ) がマスの大きさの 1e-6 を超えたら収束しなかったとして欠損にします。格子番号が格子の範囲 (0〜行数 − 1、0〜列数 − 1) を外れた点は格子の外なので欠損です (境界ちょうどの点は中に入れます)。経路の全点が格子の外ならエラー、一部だけ外なら外の点は欠損のまま描き、図の端が空きます。
• 精度: 投影面で等間隔の真の格子は接平面の上で厳密には双一次ではないので、わずかな誤差が出ます。誤差はマスの大きさに比例して小さくなり、実測で 100 km 格子で 8e-4 マス、5 km 格子 (ClimCORE と同じ) で 3e-5 マス (約 15 cm) です。
• 速さ: ClimCORE の格子で経路 260 点の格子番号が 0.02 秒です。

経緯度が 1 次元の格子

経緯度が 1 次元の通常の格子では、緯度・経度の軸ごとに独立に求まります。座標の値を昇順に並べ替えてから、点の値に対応する番号を線形に内挿し (numpy の interp)、範囲外は欠損にします。

• 昇順・降順: 番号は元の並びのものを使うので、北から南へ並ぶ緯度 (90 → −90) などもそのまま扱えます。
• 経度の規約: 点の経度を、格子の経度の最小値から 360° の範囲へ剰余で移します (格子が 0〜357.5° なら −1° は 359° として扱います)。
• 全球格子の継ぎ目: 「最後の経度 − 最初の経度 + 格子間隔 = 360°」なら経度が一周している格子とみなし、最後の列の次を最初の列 (+360°) として内挿に加えます。継ぎ目をまたぐ点 (0〜357.5° の格子での 359° など) は、最後の列と最初の列の間で内挿されます。

2 次元の格子に展開して上の方法を使っても、ほぼ同じ番号になります (2.5° 格子で差 2.7e-3 マス)。完全には一致しないので、通常の格子では 1 次元の方法を使います。

3. 周りの 4 格子点から双一次内挿する

小数の格子番号 (fj, fi) = (j0 + a, i0 + b) の位置の値を、マスの 4 つの角の値から重み付き平均で求めます。f00 は角 (j0, i0)、f01 は (j0, i0 + 1)、f10 は (j0 + 1, i0)、f11 は (j0 + 1, i0 + 1) の値です。

各角の重みは、点をはさんで反対側にある長方形の面積に等しくなります (下の図 a)。点が角の 1 つに重なればその角の値そのもの、マスの辺の上なら辺の両端の値の線形内挿になります。

• 鉛直・時刻などの次元はそのまま: 水平の 2 次元だけが経路方向の 1 次元に置き換わります (例: (lev, y, x) → (lev, path))。鉛直方向・時間方向の内挿はしません。
• 欠損の規則: 重みが正の角に 1 点でも欠損があれば、その点は欠損です。重みが 0 の角の欠損は影響させません。これで、地面の下が欠損のデータや海洋の陸の近くで存在しない値を作らず、しかも格子点ちょうどの点が隣の欠損に巻き込まれて消えることもありません。
• 格子の外: 格子番号が欠損の点 (格子の外) は結果も欠損です。
• 全球格子の継ぎ目: 1 次元の経緯度の全球格子では、最後の列の次を最初の列として継ぎ目をまたいで内挿します。
• 地形マスク: 経路断面で地形マスクを使うと、地上気圧や地形高度も同じ双一次内挿で経路の上へ写してから地面を決めます。

制約と注意

• 周期的な 2 次元座標格子の継ぎ目: 経度方向に一周している全球の 2 次元座標格子 (海洋の三極格子など) では、最後の列と最初の列の間のマスに入る点は欠損になります (格子の外として扱います)。1 次元の経緯度の全球格子は継ぎ目をまたいで内挿できます。
• 経緯度が欠損の格子点の近く: 経緯度が欠損の格子点を角に持つマスの点は欠損です。最寄りの格子点から始めたニュートン法がそういうマスを通ると、隣の正常なマスにある点も欠損になることがあります (値を作らない側に倒しています)。
• ベクトル・流線の成分: 大円の断面ではベクトル・流線の成分は選んだ変数のまま内挿して描き、断面の向きへの射影はしません (枠に注記が出ます)。
• 横軸の範囲と平均: 経路の断面では横軸の範囲は経路で決まり、水平方向の範囲平均は使えません。鉛直の範囲と、時刻などほかの次元の固定・平均は従来どおりです。

再現スクリプトでの確認

内挿に使う関数は、再現スクリプトの冒頭に本文がそのまま埋め込まれます (使うものだけ: 大円なら great_circle_points、2 次元座標の格子なら grid_fractional_indices、1 次元の格子なら grid_fractional_indices_1d、経路の断面では常に sample_bilinear、目盛の経緯度の併記では lonlat_tick_label)。アプリの描画と再現スクリプトは文字どおり同じコードで内挿するので図は一致し、スクリプトは numpy と xarray だけで動きます。2 次元座標の格子で大円の断面を描いたときの該当部分は次のようになります。

始点・終点や点の数はリテラルで書かれるので、スクリプトの側で書き換えて断面を動かせます。関数の中身を読めば、この章の手順をそのまま追えます。

検証していること

ClimCanvas のテストでは、この内挿について次を自動で確かめています (tests/test_section_path.py ほか)。

• 大円の点: 端点が一致し、距離が haversine の公式と一致し、点の間隔が等しく、全点が大円の面の上にある。日付変更線をまたいでも経度が連続。極を通る経路、始点 = 終点、対蹠点のエラー。
• 2 次元の格子番号: cartopy の投影で作った真の格子番号を復元する (100 km 格子で許容 2e-3 マス、5 km 格子で 1e-4 マス)。関数自体は投影を使わないので、独立な検算になる。180° の継ぎ目をまたぐ格子、極ステレオ格子の極そのもの、格子の外、経緯度の欠損。
• 1 次元の格子番号: 全球 (降順の緯度・経度の規約・継ぎ目)、降順の経度、領域格子の外。1 次元と 2 次元の方法の一致。
• 双一次内挿: 双一次の場 (2 + 3j − i + 0.5ij) を厳密に再現する。重み 0 の角の欠損は無視し、正の重みの角の欠損は欠損。全球の継ぎ目の内挿。
• 断面の描画: 等緯度線断面の描画配列が、経緯度が格子番号の 1 次式になる格子で解析的に逆算した値と一致する。アプリの図と再現スクリプトの図が画素一致する。
• 実データ: ClimCORE の Z (817 × 661、17 気圧面) で、等緯度線・等経度線・大円・格子の行の断面がアプリと再現スクリプトで画素一致した (1 断面あたり約 2 秒)。