楕円体は a と f だけで決まる — 離心率・曲率半径の導出

以前、姉妹アプリGeoPrism JPのブログで「GRS80とWGS84は同じ楕円体？ — 似ているようで違う2つの基準楕円体を比べる」という記事を書きました。あちらは2つの楕円体を概念として比べる内容でしたが、今回は視点を変えて、GeoCoreJPの座標変換を根っこで支えている数式そのものを扱います。今回から、座標変換の計算式を1つずつ分解していく連載を始めます。第0回は、すべての出発点になる楕円体の定数の話です。

楕円体パラメータは、実はaとfの2つだけ

測地系の話をすると「長半径」「扁平率」「離心率」「短半径」「曲率半径」といろいろな用語が出てきますが、これらはすべて独立した値ではありません。楕円体の形を決める根本のパラメータは、実は次の2つだけです。

長半径 a（赤道半径。GRS80・WGS84ともに 6378137.0 m）
扁平率 f（赤道半径に対して、極方向にどれだけ潰れているか）

離心率も短半径も曲率半径も、すべてこの a と f から計算式1本で導けます。「たくさんの定数を覚える」のではなく「2つの根本値から導く」と捉えると、楕円体まわりの数式がぐっと見通しやすくなります。

パラメータの関係図

       a（長半径）
       f（扁平率）        ← 根本の2値
        │
        ├─ e² = 2f − f²         （離心率の2乗）
        ├─ b  = a(1 − f)         （短半径）
        └─ e'²= e² / (1 − e²)    （第2離心率の2乗）
             │
             ├─ N（卯酉線曲率半径）= a / √(1 − e²sin²φ)
             └─ M（子午線曲率半径）= a(1 − e²) / (1 − e²sin²φ)^(3/2)

a と f を起点に、離心率→短半径→曲率半径という順に、すべてが式でつながっています。

aとfから導く：e²・b・第2離心率

離心率の2乗 e²

e² = 2f − f²

扁平率 f がそのまま離心率になるわけではなく、この式で変換します。GRS80では f = 1/298.257222101 なので、e² ≈ 0.00669438002290 です。

短半径 b

b = a(1 − f)

長半径から扁平率の分だけ引いた値が短半径（極半径）です。GRS80では b ≈ 6356752.31414 m となり、赤道半径との差はおよそ21.4kmです。地球はよく「ほぼ球」と言われますが、これは全周4万kmに対する比率で見れば非常に小さい潰れ方だから、という意味です。

第2離心率の2乗 e’²

e'² = e² / (1 − e²)

第2離心率は、曲率半径の式にも登場しますが、それ以上によく使われるのがガウス・クリューゲル投影（平面直角座標系の投影式）です。次回以降の連載で扱う投影式の下準備として、ここで押さえておきます。

卯酉線曲率半径Nと子午線曲率半径M — 「東西」と「南北」で曲がり方が違う

楕円体は完全な球ではないため、同じ緯度でも東西方向の曲がり方と南北方向の曲がり方が異なります。この2つを表すのが、卯酉線曲率半径 N と子午線曲率半径 M です。

卯酉線曲率半径 N：東西方向（緯度円に直交する方向）の曲率半径。N = a / √(1 − e²sin²φ)
子午線曲率半径 M：南北方向（子午線に沿った方向）の曲率半径。M = a(1 − e²) / (1 − e²sin²φ)^(3/2)

どちらも緯度φの関数になっているのがポイントです。赤道（φ=0）ではN、Mともにaに近い値を取りますが、極に近づくほど値が変化していきます。緯度・経度から距離を計算したり、平面直角座標系に投影したりする式には、必ずこのNとMのどちらか（あるいは両方）が登場します。次回予告している「子午線弧長」の記事は、このMを緯度0から積分することで求まる値です。

GRS80とWGS84の定数差 — fのわずかな違い

GeoCoreJPが内部で扱う楕円体は主にGRS80ですが、GPSの世界標準であるWGS84と比較されることがよくあります。両者の長半径 a はまったく同じ 6378137.0 m です。違いは扁平率 f にあります。

楕円体	長半径 a (m)	扁平率 f
GRS80	6378137.0	1 / 298.257222101
WGS84	6378137.0	1 / 298.257223563

fの差は小数点以下9桁目でようやく現れる程度で、これを短半径bの差に換算すると0.1mmのオーダーです。実務上の座標計算では両者を区別する必要がほぼないほど近い値ですが、「同じ」ではなく「定義上は極めて近い別の値」という点は、厳密な精度検証をする際には意識しておく必要があります。

コード実装：定数テーブルの持ち方

GeoCoreJPでは、楕円体ごとの根本パラメータ（aとf）だけを定数として持ち、e²・b・e’²といった派生値は必要な計算のたびに導出する、あるいは初期化時に一度だけ計算してキャッシュする形にしています。派生値をハードコードしてしまうと、桁数の丸め方によって計算結果がわずかにブレる可能性があるため、「根本の2値だけを信頼できる定数として持ち、それ以外は式で導く」という設計方針を徹底しています。

struct EllipsoidConstants {
    let a: Double   // 長半径
    let f: Double   // 扁平率

    var e2: Double { 2 * f - f * f }               // 離心率の2乗
    var b: Double { a * (1 - f) }                    // 短半径
    var e2prime: Double { e2 / (1 - e2) }            // 第2離心率の2乗

    func primeVerticalRadius(latRad: Double) -> Double {
        let sinPhi = sin(latRad)
        return a / sqrt(1 - e2 * sinPhi * sinPhi)    // N
    }

    func meridianRadius(latRad: Double) -> Double {
        let sinPhi = sin(latRad)
        let denom = pow(1 - e2 * sinPhi * sinPhi, 1.5)
        return a * (1 - e2) / denom                  // M
    }
}

static let grs80 = EllipsoidConstants(a: 6378137.0, f: 1.0 / 298.257222101)
static let wgs84 = EllipsoidConstants(a: 6378137.0, f: 1.0 / 298.257223563)

こうしておくと、楕円体を追加・変更したいとき（たとえば旧測地系で使われるベッセル楕円体を扱う場合など）も、aとfの2値を差し替えるだけで済み、e²・b・N・Mの計算式には一切手を入れずに済みます。