【proj】地心直交座標系(X,Y,Z)と地理座標+楕円体高(lon,lat,height)の相互変換

projというライブラリはprojctionという名前の通り、地図投影法の変換(projコマンド)と空間参照系同士の変換(cs2csコマンド)を実行するものである。
空間参照系同士の変換の一部に地心直交座標系と地理座標+楕円体高に変換するコードがあるのでそれを利用する。

使い方

詳しい使い方はこちらにある

(lon,lat,height)->(X,Y,Z)への変換

cct +proj=cartを使う

echo 140.0878550027 36.1037747783 65.8397 |
    cct -t 0 +proj=cart |
    awk '{print $1,$2,$3}'
# -3957314.6200 3310254.1370 3737540.0430

解説
cctコマンドは3次元座標に加えて時刻がある４次元座標を扱うコマンドになる、-t 0オプションにより時間を0にして3列で入力が可能になる。それでも出力が4列(4列目は時刻0)となるのでawkで4列目を落とす。

(X,Y,Z)->(lon,lat,height)への変換

cct -I +proj=cartを使う

echo -3957314.6200 3310254.1370 3737540.0430 |
    cct -t 0 -I +proj=cart |
    awk '{print $1,$2,$3}'
# 140.0878550027 36.1037747783 65.8397

解説
-Iで逆変換を行う。それ以外は同じ。

おまけ

WebAPIを使う

ここの緯度・経度と地心直交座標の相互換算のAPIを使う

※ただし大量リクエストはご法度。同一IPアドレスからのリクエストは、10秒間で10回。10行ぐらいのデータをいれてください。

プログラム(xyz->lbh)の変換。フィルタコマンドとして作成した。

xyz2lbh.py

#!/usr/bin/env python3
import sys

import requests

baseurl = 'http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/trans.pl'

for line in sys.stdin:
    x, y, z = line.split()
    params = {
        "outputType": "json",
        "cnv_type": 1,
        "geocentricX": float(x),
        "geocentricY": float(y),
        "geocentricZ": float(z)
    }
    r = requests.get(baseurl, params=params)
    outputJson = r.json()['OutputData']
    # これらはすべて文字列
    longitude = outputJson['longitude']
    latitude = outputJson['latitude']
    ellipsoidHeight = outputJson['ellipsoidHeight']
    print(f'{longitude} {latitude} {ellipsoidHeight}')

echo -3957314.620 3310254.137 3737540.043 | xyz2lbh.py
# 140.08785500  36.10377478  65.840