地球上のある地点からある距離ある方向へ離れた点の座標を求める（大円距離の逆を解く）

はじめに

より具体的には，ある任意の地点から出発し，ある方位へXXX km移動して辿り着く座標はどこか，という問題である．2点間の大円距離を求める公式を知っている人は，その既知と未知の部分を取り替えているだけなので簡単に解けそうと感じるはず．

しかし，簡単な解法でもあるのだろうとネットを探せども，なかなかヒットしない．

(追記：2023/08/04実は以下の楕円体における移動問題を解くモジュールがあるらしいこちらのページ．ただし，以下のような楕円体における解であり，球体解ではない．)

こちらのページで同様な質問がなされているが，回答では平面の三角形における単純な三角比が説明されている．確かに地球の曲率を無視できるくらい近距離であればそれで良いのだけど，遠くの点の座標は誤差を無視できない．

ちなみに，地球を楕円体として考えるとこれは非線形問題となるらしく，数値解を得る方法はいくつか議論されているらしい．こちらのサイトではある解法を利用しWeb上で計算してくれるようだ．
Direct 距離と方位角から楕円面上で到達点を求める - mapli.net
これはこれで便利だが，地球を球体と見做した場合のシンプルな解を知りたい，

ということで，球面三角法により大円公式の逆を解き，おお手軽なPythonモジュールを作成する．

2点間の大円距離を求める公式

ここでは便宜的に「大円公式」と呼ぶ．この公式は大変有名であり，これは2地点の座標（(緯度，経度)=$ (\phi_1, \lambda_1), (\phi_2, \lambda_2) $）を与えることで大円距離$ D $を返す．

D=Rc=R\cos^{-1}(\sin \phi_1 \sin \phi_2 + \cos \phi_1 \cos \phi_2 \cos(\lambda_2 - \lambda_1))

ここで， $ R $ は球体半径，$ c $ は球体中心から2点を望む中心角である．こちらのページ（地球上の2点間の距離の求め方 - Qiita）で分かりやすく説明されている通り．上記の大円公式は，球面三角法の余弦定理（詳細は後ほど）を応用することで定義される．
（なお，サイトの図A5では中心角は $ x $ とあるが説明にある $ c $ と等価であり，こちらを採用．また緯度，経度も，文章は $\theta, \phi$ だが図は $\phi, \lambda$ である．こちらは図の表記を用いた．）

ひとまず，上の大円公式を変形すると，

c=\frac{D}{R}

を得る．つまり，移動距離$D$を知っていれば，中心角は簡単に求めることができる．$\cos^{-1}$を解くと，

\cos c = \sin \phi_1 \sin \phi_2 + \cos \phi_1 \cos \phi_2 \cos(\lambda_2 - \lambda_1)

となる．これは，球面三角法の余弦定理そのものである．
ここで注意なのは， $ D/R $ は単位がラジアンであるが， $ \theta, \lambda $ は度である．混乱を防ぐため，計算は度に統一し， $ D/R $ を度へ変換した値を

c = \frac{180D}{\pi R}

と再定義し，「cの単位は度」と思って利用する．（ラジアンに統一すると数式が煩雑となるので，許してほしい．．）

球面三角法の余弦定理

この説明はこちらのPDFが詳しい
（https://inertial-navigator.net/wp-content/uploads/2019/10/TB003_00.pdf）．

このPDFの図7にある，3つの大円の作る球面三角形ABCにおいて，弧AB，弧BC，弧CAにそれぞれ対応する中心角c，a，bと， $ A = \angle DAE $ により定義される角Aとの間に，

\cos a = \cos b\ cos c + \sin b \sin c \cos A

が成り立つことが球面三角法の余弦定理であるらしい．つまり，2つの弧に対応する中心角とその2つの弧の間の角が分かれば，残された弧に対応する中心角の余弦を得る．

本題

上記にある2つのサイトの2つの図，これらの図はとても良くできており，これらを組み合わせることで，球面三角法の余弦定理を用いて問題を解くことができた．図に加筆することで説明を試みる（手書きでごめんなさい）．加筆部分は赤線である．

方位角を考える

図A5内の地点Aを出発点（(緯度,経度)=$(\phi_1,\lambda_1)$），地点Bを到達点（(緯度,経度)=$(\phi_2,\lambda_2)$），その間の移動距離を$d$，地球半径を$R$，出発点における方位角を$ \theta $（東向きを0，反時計回りを正とする）とする．
既知の値は$ (\phi_1,\lambda_1), d, R, \theta $であり，これらを使い，$(\phi_2,\lambda_2)$を求める．

点Aを通る緯度円をぐるっと書いてみたのだが，この緯度円に対し，進向する方向との間の角が$\theta$であり，$\angle BAN$は$90º-\theta$となる（弧ANは真北を向いているため）．また地球中心付近にはAとBを望む中心角cを書き加えてある．

中心角と緯度軽度の対応を考える

こちらの図7のAとBを先ほどの図A5の地点Aと地点Bに対応させることができる．中心角cも対応している．また，北極点NはCに対応する．ここで，図7の中心角cと中心角bを赤ぺんでなぞっているが，これは「既知」であることを示している．というのも，cは最初の方で求めた通り$c=180d/\pi R$（Dをdとした）であり，bは天頂角ゆえ図A5にある「$90-\phi_1$」に対応しており，そもそも$\phi_1$は出発点の緯度であるからである．

今，弧ABと弧ACに対応するそれぞれの中心角を知ったわけだが，弧ABと弧ACの間の角$A = \angle DAE$（$=\angle BAC$，多分）をもし知っていれば球面三角法の余弦定理により弧BCに対応する中心角aを得る．Cは北極点であるので，中心角aはやはり地点Bの天頂角，つまり図A5の「$90 - \phi_2$」に対応するので，つまり地点Bの緯度を得る．

よく見ると，図7の$\angle BAC$は図A5の地点Aの $ \angle BAN $ に対応するので，実は角Aも既知であり，方位角$\theta$を用いて$90º-\theta$と表される．よって上段の議論により地点Bの緯度 $\phi_2$ を求めることができる．球面三角法の余弦定理より，

\cos a = \cos b \cos c + \sin b \sin c \cos A \\
\Leftrightarrow \cos (90º-\phi_2) = \cos (90º-\phi_1) \cos c + \sin (90º - \phi_1) \sin c \cos (90º-\theta) \\

\Leftrightarrow

 90º-\phi_2 = \cos^{-1} \bigl( \cos (90º-\phi_1) \cos c + \sin (90º - \phi_1) \sin c \cos (90º-\theta) \bigr)

よって

\phi_2 = 90º - \cos^{-1} \bigl( \cos (90º-\phi_1) \cos c + \sin (90º - \phi_1) \sin c \cos (90º-\theta) \bigr)

ただし，$c=180d/\pi R$

経度については，図A5では$\angle ANB = \lambda_2 - \lambda_1$であり，それは図7では$\angle ACB$に対応している．今度は角C（$= \angle ACB$）に注目して球面三角形の余弦定理を考えると，$\phi_2$も使って，

\cos c = \cos a \cos b + \sin a \sin b \cos C

\Leftrightarrow

\cos c = \cos (90º - \phi_2) \cos (90º - \phi_1) + \sin (90º - \phi_2) \sin (90º - \phi_1) \cos (\lambda_2 - \lambda_1)

\Leftrightarrow

\cos (\lambda_2 - \lambda_1) = \frac{\cos c - \cos (90º - \phi_2) \cos (90º - \phi_1)}{\sin (90º - \phi_2) \sin (90º - \phi_1)}

よって

\lambda_2 = \lambda_1 \pm \cos ^{-1} \bigl( \frac{\cos c - \cos (90º - \phi_2) \cos (90º - \phi_1)}{\sin (90º - \phi_2) \sin (90º - \phi_1)} \bigr)

ただし，$\phi_1 \ne \pm 90º$ ， $\phi_2 \ne \pm 90º$ ， $-180º < \lambda_2 - \lambda_1 < 180º$，$c=180d/\pi R$

$\cos$は偶関数なので$\cos^{-1}$は正も負もとりうる．進行方向の方位角$\theta$を知っているので，東に進むなら正，西に進むなら負とすれば良いだろう．（実際にプログラム上はそれで動作する）

Pythonによるサンプルコード

参考までにpythonモジュールを晒しておきます．保証はしませんが（ちょっと精度悪い気がする）．注！入力は度，内部の計算はラジアン，出力を度としています．

import numpy as np

def invert_gcd(lon=135., lat=35., theta=0., d=1000.):
    '''
    invert great circle distance and get new coordinate
    using Spherical Trigonometry
    kasuga 2022.04.10

    Parameters
    ----------
    lon, lat: coordinate of the start point [degree]
    theta: angle from east toward the end point [degree]
    d: travel distance [km]

    Return
    ------
    coordinate of the end point (lon_end, lat_end) [degree]
    '''

    rr = 6371.  # earth radius [km]

    # convert angles within +-360
    theta = np.divmod(theta, 360.)[1]

    A = np.deg2rad(90.-theta)
    c = d/rr  # center angle between pos. vectors of start and end [radian]
    b = np.deg2rad(90.-lat)  # zenith of start point

    # cosine formula of Spherical Trigonometry
    cos_a = np.cos(b)*np.cos(c) + np.sin(b)*np.sin(c)*np.cos(A)
    a = np.arccos(cos_a)  # zenith of end point [radian]
    lat_end = 90. - np.rad2deg(a)

    _N = np.cos(c) - np.cos(a)*np.cos(b)
    _D = np.sin(a)*np.sin(b)
    cos_C = np.divide(_N, _D, where=_D!=0.)  # ignore RuntimeWarning

    cos_C = np.where(cos_C>1., 1., cos_C)
    cos_C = np.where(cos_C<-1., -1., cos_C)
    C = np.arccos(cos_C)  # memo: np.arccos() retruns positive value

    lon_end = np.where((-90.<theta)&(theta<90.)|(270.<theta),
                       lon+np.rad2deg(C),
                       lon-np.rad2deg(C))

    # convert angles within +-360
    lon_end = np.divmod(lon_end, 360.)[1]

    lon_end = np.where(np.sin(a)==0., 0., lon_end)  # end point is North/South pole
    lon_end = np.where(np.sin(b)==0., theta, lon_end)  # start point is North/South pole

    return (lon_end, lat_end)

作図して確認してみます．緯度35ºN，経度135ºEから1000 km移動した地点を，方位角を15度ずつ回転させて一周させてみます．

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

from invert_gcd import invert_gcd  # これが自作モジュール 

lat, lon = 35, 135
radius = 1000

data_proj = ccrs.PlateCarree()
map_proj = ccrs.Orthographic(central_latitude=lat, central_longitude=lon)

fig = plt.figure()
ax = plt.axes(projection=map_proj)

ax.coastlines()
ax.gridlines()

# center point plot
ax.plot(lon, lat, marker='.', c='r', transform=data_proj)

for theta in np.arange(0, 360, 15):

    lon_end, lat_end = invert_gcd(lon=lon, lat=lat, theta=theta, d=radius)
    print(theta, lon_end, lat_end)

    # around point plots
    ax.plot(lon_end, lat_end, marker='.', c='r', transform=data_proj)

plt.savefig('test.png')

こうして見る分にはいい感じそうです．一応，極点，赤道，子午線を含む場合も試しましたが問題ありませんでした．

まとめ

既知の値$ (\phi_1,\lambda_1), d, R, \theta $を使い，未知の座標$(\phi_2,\lambda_2)$は，

\phi_2 = 90º - \cos^{-1} \bigl( \cos (90º-\phi_1) \cos c + \sin (90º - \phi_1) \sin c \cos (90º-\theta) \bigr) \\

\lambda_2 = \lambda_1 \pm \cos ^{-1} \bigl( \frac{\cos c - \cos (90º - \phi_2) \cos (90º - \phi_1)}{\sin (90º - \phi_2) \sin (90º - \phi_1)} \bigr)

のように与えられる．ただし，$c=180d/\pi R$

間違い，勘違いありましたらご指摘していただきますと幸いです．m(_ _)m

（2022/04/11）$\phi_2$に関し，絶対値の表記を$\pm$へ修正しました．
（2022/11/24）モジュールの修正，経度が360度を越すバグの修正

参考文献

Direct 距離と方位角から楕円面上で到達点を求める - mapli.net
地球上の2点間の距離の求め方 - Qiita
球面三角法（PDF）https://inertial-navigator.net/wp-content/uploads/2019/10/TB003_00.pdf