第3章: WGS84座標系とENU

はじめに

WGS84・ECEF・ENUはGIS技術者に馴染み深い座標系ですが、GISソフトが内部処理してくれるため「中身」を自前で実装する機会は少ないのではないでしょうか。この章では、普段ブラックボックスとして利用している座標変換を一つひとつ組み立てます。

3Dエンジン固有の注意点として、Three.jsのY-up規約とECEFのZ-up（北極方向）の軸マッピングがあります。都市マーカーをクリックすると、各座標系での値を確認でき、変換パイプラインの動作を実感できます。

この章で学ぶこと

• WGS84楕円体の定義とパラメータ —— GISの基礎をThree.jsの球体として可視化
• ECEF座標系: 地球中心固定の直交座標 —— GPSの内部座標系を3D空間で体験
• 測地座標（緯度・経度・高度）からECEFへの変換 —— 卯酉線曲率半径 N の幾何学的意味
• ECEF → Three.jsシーン座標への軸変換 —— Y-up（Three.js）とZ-up（ECEF）の変換ルール
• ENU座標系: 地表上の局所座標系 —— 東西南北の方向感覚を3Dで確認
• ECEF ↔ ENU変換の回転行列 —— 直交変換行列の導出と実装

デモの表示要素

なぜ楕円体が必要か

地球は完全な球ではありません。自転による遠心力で赤道方向にやや膨らんだ回転楕円体に近い形状をしています。 この形状を数学的にモデル化したものが準拠楕円体（reference ellipsoid）です。

GPSが使用するWGS84（World Geodetic System 1984）は、地球の形状を最もよく近似する標準的な準拠楕円体を定めています。 GISで緯度経度を扱う際、実質的にWGS84楕円体上の座標を意味します。

WGS84楕円体パラメータ

WGS84楕円体は、2つの基本パラメータで定義されます。

赤道半径と極半径の差は約21km。地球の直径（約12,756km）に対して0.3%程度の差ですが、 精密な座標変換では無視できません。

export const WGS84 = {
  /** 長半径（赤道半径）[m] */
  a: 6_378_137.0,
  /** 扁平率 */
  f: 1 / 298.257_223_563,
  /** 短半径（極半径）[m] */
  get b(): number {
    return this.a * (1 - this.f);
  },
  /** 第一離心率の二乗 */
  get e2(): number {
    return 2 * this.f - this.f * this.f;
  }
} as const;

e²の計算で a や b を使わず f から直接導出しているのは、浮動小数点精度を維持するためです。

左のデモでは Y軸方向に WGS84.b / WGS84.a のスケールを掛けて、 この扁平形状を再現しています。透明度を下げて内部構造が見えるようにしています。

ECEF座標系とは

ECEF（Earth-Centered, Earth-Fixed）は、地球の中心を原点とし、地球と一緒に回転する直交座標系です。 GPSの衛星から受け取る位置情報は内部的にECEFで処理されています。

測地座標からECEFへの変換

緯度・経度・高度からECEFへの変換では、楕円体の卯酉線曲率半径（N）を使います。 N は楕円体表面のその緯度における「東西方向の曲率半径」で、 楕円体表面のある点から法線方向にZ軸（自転軸）まで延ばした距離です。

Z成分の式に注目してください。X、Y成分が (N + h) であるのに対し、Z成分には (1 - e²) が掛かっています。 これは楕円体が扁平であることに起因します。もし e² = 0（完全な球）なら、(1-e²) = 1 となり、この違いは消えます。

export function geodeticToEcef(
  lat: number, lon: number, alt: number
): EcefCoord {
  const phi = degToRad(lat);
  const lam = degToRad(lon);
  const sinPhi = Math.sin(phi);
  const cosPhi = Math.cos(phi);
  const sinLam = Math.sin(lam);
  const cosLam = Math.cos(lam);

  const N = WGS84.a / Math.sqrt(
    1 - WGS84.e2 * sinPhi * sinPhi
  );

  return {
    x: (N + alt) * cosPhi * cosLam,
    y: (N + alt) * cosPhi * sinLam,
    z: (N * (1 - WGS84.e2) + alt) * sinPhi
  };
}

ECEFから測地座標への逆変換

逆変換（ECEF → 測地座標）は順変換ほど単純ではありません。 実用上はBowringの反復法が広く使われています。 本実装では補助角 θ = atan2(z·a, p·b) による初期推定が十分に高精度であるため、 実質的に1回の計算でmm以下の精度が得られます。

処理の流れ:

• 経度: atan2(Y, X) で直接求まる（ECEFのX-Y平面での角度 = 経度）
• 補助角: atan2(z*a, p*b) で幾何学的に推定（この初期値だけで十分な精度）
• 緯度の計算: 補助角 θ を用いて緯度 φ を算出
• 楕円体高: 赤道付近と極付近で計算式を使い分け

export function ecefToGeodetic(
  x: number, y: number, z: number
): GeodeticCoord {
  const p = Math.sqrt(x*x + y*y);
  const lam = Math.atan2(y, x);

  // Bowring初期値
  let theta = Math.atan2(z*a, p*b);
  let phi = Math.atan2(
    z + ep2 * b * Math.sin(theta)**3,
    p - e2 * a * Math.cos(theta)**3
  );

  // 3回反復で十分な精度
  for (let i = 0; i < 3; i++) {
    theta = Math.atan2(z*a, p*b);
    phi = Math.atan2(/* ... */);
  }
  // ... 高度の計算
}

座標の確認

左のデモで都市マーカーをクリックすると、 各座標系での値を確認できます。変換パイプラインの動作を実感してください。

都市マーカーをクリックして座標を表示

なぜ軸変換が必要か

ECEF座標系とThree.jsのシーン座標系では、軸の向きが異なります。

これを合わせるため、ECEF座標をThree.jsのシーン座標に変換する軸変換が必要です。

シーン座標系へのマッピング

符号の反転（-Y）は、右手座標系を維持するために必要です。 ECEF（右手系）からThree.js（右手系）への変換で、2軸を入れ替える場合、 1つの軸を反転しないと左手系になってしまいます。

export function ecefToScenePosition(
  x: number, y: number, z: number
): EcefCoord {
  return {
    x: x,     // ECEF X → シーン X
    y: z,     // ECEF Z → シーン Y（上）
    z: -y     // ECEF Y → シーン Z（反転）
  };
}

本書では 1 Three.js単位 = 1メートル とし、スケーリングを行いません。 CesiumJSも同様に1:1のスケールを採用しています。 スケーリングを導入すると、タイルテクスチャの解像度計算がスケール依存になったり、 3D Tilesなど外部データとの座標整合が複雑になるためです。

行列形式での軸変換

軸変換は4x4行列として表現することもできます。 3D Tilesアダプターではこの行列形式が使われています。

// ECEF → シーン座標の軸変換行列
const ecefToScene = new THREE.Matrix4().set(
  1,  0,  0,  0,   // X → X
  0,  0,  1,  0,   // Z → Y
  0, -1,  0,  0,   // -Y → Z
  0,  0,  0,  1
);
this.group.applyMatrix4(ecefToScene);

第9章で実装する独自の3D TilesアダプターはECEF座標でモデルを配置するため、Groupに軸変換行列を適用することで、すべての子オブジェクトが自動的にシーン座標に変換されます。

ENU座標系

ENU（East-North-Up）は、地球上のある基準点を中心としたローカル直交座標系です。

基準点が変わるとENU座標系も変わるため、「ローカル」座標系と呼ばれます。 ENUは以下の場面で使われます:

• 3D Tilesのtilesetが持つ局所座標のECEFへの変換
• 地表付近でのオブジェクト配置（建物、道路など）
• 測量データの処理、GPSの相対位置計算

ECEF ↔ ENU変換の回転行列

ECEF→ENU変換は、差分ベクトルに回転行列 R を掛けることで実現します。

各行の意味:

• 第1行（East）: 東方向の単位ベクトル。経度方向に直交する水平ベクトル
• 第2行（North）: 北方向の単位ベクトル。緯度方向かつ水平な方向
• 第3行（Up）: 楕円体の法線方向。ECEFでの法線ベクトルに一致

逆変換（ENU → ECEF）は回転行列の転置（直交行列なので転置 = 逆行列）を使います。

export function enuToEcef(
  e: number, n: number, u: number,
  latRef: number, lonRef: number,
  altRef: number
): EcefCoord {
  const ref = geodeticToEcef(
    latRef, lonRef, altRef
  );
  const phi = degToRad(latRef);
  const lam = degToRad(lonRef);
  // 回転行列の転置（逆変換）
  return {
    x: ref.x
      + (-sinLam*e - sinPhi*cosLam*n
         + cosPhi*cosLam*u),
    y: ref.y
      + (cosLam*e - sinPhi*sinLam*n
         + cosPhi*sinLam*u),
    z: ref.z
      + (cosPhi*n + sinPhi*u)
  };
}

ENU軸の描画

デモでは enuToEcef で各軸の単位ベクトルをECEFに変換し、 原点との差分から方向を求めてラインを描画しています。

// 東方向の単位ベクトルをECEFに変換
const tipEcef = enuToEcef(
  1, 0, 0,       // E=1, N=0, U=0
  tokyoLat, tokyoLon, 0
);
// 原点との差分で方向を取得
const dir = tip.sub(origin).normalize();
// 方向にenuAxisLenを掛けてラインの端点に
const end = tokyoOrigin.clone()
  .add(dir.multiplyScalar(enuAxisLen));

往復変換の精度

実装の正しさは、往復変換（round-trip）テストで検証できます。 テストでは、赤道上、極、東京、南半球の各点で往復変換を行い、 緯度・経度の誤差が 10^-10 度以下であることを確認しています。

it('東京のラウンドトリップ', () => {
  const ecef =
    geodeticToEcef(35.6762, 139.6503, 100);
  const result =
    ecefToGeodetic(ecef.x, ecef.y, ecef.z);
  expect(result.lat)
    .toBeCloseTo(35.6762, 9);
  expect(result.lon)
    .toBeCloseTo(139.6503, 9);
  expect(result.alt)
    .toBeCloseTo(100, 3); // 1mm以内
});

まとめ

この章では、GISの基盤となる座標系と座標変換を実装しました。

• WGS84パラメータ: 長半径 a, 扁平率 f から短半径 b, 離心率 e² を導出
• 順変換: 卯酉線曲率半径 N を用いた測地座標 → ECEF変換
• 逆変換: Bowringの反復法によるECEF → 測地座標変換（3回反復）
• 軸変換: X→X, Z→Y, -Y→Z でECEF → シーン座標
• ENU座標系: 基準点を中心とした東北上のローカル座標と回転行列

次章では、地図タイルの基盤となるWebメルカトル投影法とタイル座標系を解説します。