GPSなどの衛星測位計算を行っていると、地球上の位置を表すのに

ECEF (earth-centered earth-fixed)で定められるＸＹＺの直交座標系
geodetic datum と呼ばれる緯度、経度、高さで表す方法

の２つが用いられます。高さは地球を楕円体でモデル化した時の球面からの鉛直方向で計算したものを「楕円体高」と言い、ジオイド（平均海水面に相当）からの高さを標高と言います。

参考資料はたくさんあり、

WGS84と座標変換のはなし
WikiPediaの　経緯度

などたくさんあります。

おじさんは今まで自分で実装したライブラリを用いていたのですが、今日見たらpython のライブラリがあったのでメモします。

pymap3d を用いた座標変換の動作確認

pymap3d

https://github.com/geospace-code/pymap3d でソースあり
https://pypi.org/project/pymap3d/ に登録あり。pip で問題なくインストールできた。

目的の座標変換の関数は pymap3d.ecef.ecef2gedetic にあります。名前の通り。

動作確認

電子基準点のつくば１の座標値を使用します。

from math import isclose
from pytest import approx
from pymap3d.ecef import ecef2geodetic

def test_ecef2geodetic():
    x_ecef = [-3957162.4119, 3310203.4927, 3737752.2980]
    x_geod = [36.106112803, 140.08720184, 70.336]
    x_geod_from_ecef = ecef2geodetic(x_ecef[0], x_ecef[1], x_ecef[2])
    print( "x_geod_from_ecef=", x_geod_from_ecef)
    print( "x_geod=", x_geod)
    assert x_geod[0] == approx( x_geod_from_ecef[0] )
    assert x_geod[1] == approx( x_geod_from_ecef[1] )
    assert isclose(x_geod[2], x_geod_from_ecef[2], rel_tol=0.01)

いざ pytest --capture=no -v をすると警告がなぜか出ますが、ＯＫでした。

x_geod_from_ecef= (36.10611280179599, 140.08720184317684, 70.33612007893962)
x_geod= [36.106112803, 140.08720184, 70.336]
PASSED

これからはこれを薦めよう。