緯度,経度,GPS高度を三次元直交座標に変換する

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

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

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

#緯度経度の表現について
まず前提というか，前置きとして緯度経度の表現について触れておきます．
緯度経度の記述には主に以下の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](https://qiita.com/sw1227/items/e7a590994ad7dcd0e8ab)

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

#求めかた
まず以下の図のように変数を置きます．こういうの思いつく人ってすごいよね．
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/284083/9a05abbf-a907-5fa7-bf49-9ac3f32b1c6d.png" width=70%>
`緯度Latitude φ`, `経度Longitude λ`, `楕円体高 h`, あと点Pの座標`x`, `y`, `z`です．中央の楕円体が地球ですね．この楕円体にGPSはWGS-84を使っているというわけです．

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

###補足:楕円体高とは
その名の通り，地球に見立てた**楕円体の表面から測った**高度です．GPSから返される高度はこの楕円体高．
ここで注意が必要なのは，楕円体高はよく言う**標高とは別物**だということ．標高は**平均海面から測った**高度です．平均海面を陸まで伸ばしたラインをジオイドといい，**標高**は目標物とジオイドの間の距離だといえます．そんで楕円体とジオイドの間の距離を**ジオイド高**っていうんですね．つまり**楕円体高は標高とジオイド高の和で表せる**ってこと．ややこしいね．
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/284083/5cc1845f-ecdf-47c0-1d37-64cce2c73d66.png" width=100%>
図解したらものすごい色遣いになりました．
ジオイド高は地球上の地点によって異なり，気象庁の計算サイトなどから得ることが出来ます．
https://vldb.gsi.go.jp/sokuchi/surveycalc/geoid/calcgh/calc_f.html

###式は
導出過程を全部書くのは大変なので**全部端折ります**と，点Pの座標は次のように表すことが出来ます．
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/284083/015511ef-0f19-1e51-0475-ddc43384d87d.png" width=100%>
ここで`aは長半径`, `eは離心率`で，WGS-84準拠楕円体の場合は以下の定数を用います(ほんとはeは扁平率から計算します)．

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

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

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

```python
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]
```
以上です．駄文失礼しました．