緯度経度を直交座標系における位置に変換したり，緯度経度から直線距離を出したりしたいシチュエーションって結構ありますよね(ない)．

実際この手の記事はほかの方も多く書かれていると思います．ただこう，用途によって使いたい座標系が違ったりすると思うんですね．ということで，僕はこのたび地球の重心を原点にした移動しない三次元直交座標系に変換する方法を学んだので備忘録の意味も込めてここにまとめておきたいと思います．

ややアバウトですのでその道の方でご指摘があればお願いいたします．

緯度経度の表現について

まず前提というか，前置きとして緯度経度の表現について触れておきます．
緯度経度の記述には主に以下の2種類の表現が使われています．混同しないように注意が必要です．

度分秒による表記

60進法による135° 12' 34.56"というような形式．
よく耳にするやつですね．「135度12分34秒56」などと読みます．

135° 12' 34" 56と表記したり，「135度12分34.56秒」と読んだりすることもあります．それから1351234.56とか135.12.34.56と書くこともあるようです．
個人的には本当に紛らわしいのでやめてくれの一言ですが，ピリオドの付き方で次の十進法表記と区別することが出来ます．

十進法表記

こちらはその名の通り10進法による表記で，135.2096というようなもの．「分」「秒」の部分を「度」の少数で表しています．
上の例の135° 12' 34.56"は，

135 + (12 / 60) + (34.56 / 3600)

という計算をして135.2096という十進数表記に変換されます．読み方は「135.2096度」です．

一部の度分秒表記はとても十進数表記に似た姿をしていますが，小数点の付き方を見ればどちらで表現しているか見分けられることがお分かりいただけると思います．

ここで求める座標について

この記事ではGPSで測位された緯度経度の情報を，絶対的な三次元座標系における位置に変換する方法について解説しています．重要なポイントはただの三次元座標を得る手順を説明しているということです．

つまり何が言いたいかというと実際に地球がそこにあるかは関係ないってこと．
たとえこの方法で2つの任意の点の座標を得たとして，それらの2点間の距離を求めても地球を貫通した直線距離になります．だから飛行機とか宇宙機のことを考えるときにこの記事の方法は有用です(地表に沿って移動しないからね)．

地上の建造物や自動車とかのことを考えるときは，たとえば任意の点を原点としてある地点の二次元座標を得る(しかも地球の局面を考慮する)方法などが有用でしょう．こちらの記事がすごく参考になります．

緯度経度と平面直角座標の相互変換をPythonで実装する - Qiita

WGS-84座標系とは

今回計算に用いるのは，GPSの衛星たちが使っているWGS-84という三次元直交座標系です．
地球の重心を原点とする地球固定直交座標系の一種です．ようするに地球と一緒にぐるぐる回るってこと．
WGS-84準拠楕円体のパラメータを以下に示...そうと思ったけど必要なら各自おググりください．

求めかた

まず以下の図のように変数を置きます．こういうの思いつく人ってすごいよね．

緯度Latitude φ, 経度Longitude λ, 楕円体高 h, あと点Pの座標x, y, zです．中央の楕円体が地球ですね．この楕円体にGPSはWGS-84を使っているというわけです．

2020.04.01編集
楕円体高と標高の英訳はゆらぐ場合があるようなので削除しました．知見のある方がいらっしゃいましたらコメント等いただけますと幸いです．

補足:楕円体高とは

その名の通り，地球に見立てた楕円体の表面から測った高度です．GPSから返される高度はこの楕円体高．
ここで注意が必要なのは，楕円体高はよく言う標高とは別物だということ．標高は平均海面から測った高度です．平均海面を陸まで伸ばしたラインをジオイドといい，標高は目標物とジオイドの間の距離だといえます．そんで楕円体とジオイドの間の距離をジオイド高っていうんですね．つまり楕円体高は標高とジオイド高の和で表せるってこと．ややこしいね．

図解したらものすごい色遣いになりました．
ジオイド高は地球上の地点によって異なり，気象庁の計算サイトなどから得ることが出来ます．
https://vldb.gsi.go.jp/sokuchi/surveycalc/geoid/calcgh/calc_f.html

式は

導出過程を全部書くのは大変なので全部端折りますと，点Pの座標は次のように表すことが出来ます．

ここでaは長半径, eは離心率で，WGS-84準拠楕円体の場合は以下の定数を用います(ほんとはeは扁平率から計算します)．

• a = 6378137 [m]
• e = 0.0818191908426 [-]

以上の式はその道の方の間では有名(?)なものだそうで，探せば導出はいろいろ出てくると思います．必要に応じておググりください．この記事を参考文献にしたらだめだよ．

Pythonで実装すると

僕はこんな感じで関数にしておくことにしたよ．
　引数 : 緯度φ[deg], 経度λ[deg], 楕円体高h[m]
返り値 : x座標[m], y座標[m], z座標[m]

import numpy as np

def xyz(phi, lamb, h):  #引数の緯度経度は10進数表記

    # 緯度経度をラジアンに直す
    phi = np.deg2rad(phi)
    lamb = np.deg2rad(lamb)

    #定数(WGS-84準拠楕円体の場合)
    a = 6378137          #長半径 [m]
    e = 0.0818191908426  #離心率 [-]

    N = a / np.sqrt(1 - (e ** 2 * np.sin(phi) ** 2))
    x = (N + h) * np.cos(phi) * np.cos(lamb)
    y = (N + h) * np.cos(phi) * np.sin(lamb)
    z = (N * (1 - e ** 2) + h) * np.sin(phi)

    return x, y, z  #[m]

以上です．駄文失礼しました．