緯度・経度の差分から移動距離を計算する方法（２）

前回は、計算や考え方を簡単にするために、完全球体で考えましたが、今回は回転楕円体と向き合って考えてみることにします。

２．きちんと回転楕円体として考える例

　地球は、楕円が回転している、回転楕円体として考えていきます。赤道半径を$a$、極半径を$b$とします。

　赤道面が$xy$平面にあり、極が$z$軸方向にあるとして軸を取ると、回転楕円体の表面を表す方程式は

$$\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1$$

として表すことができます。

　一般的には、赤道半径と極半径の比率、扁平率を

$$f=\frac{a-b}{a}$$

を用いて、

$$\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{a^2(1-f)^2}=1$$

$${x^2}+{y^2}+\frac{z^2}{(1-f)^2}=a^2\tag{4}$$

として表しておきます。

２－１．生真面目に方程式を解く方法

　(4)の回転楕円体の表面上にある点$A$から点$B$への移動として考えます。点$A$の緯度を$\theta$、点$A$と点$B$との間の緯度の差を$\gamma$、経度の差を$\delta$とします。点$A$の経度は結果には反映されないので、ゼロとします。すなわち、点$A$は\x\軸上にあることになります。

　点$A$の座標は、

$$A(x_a, y_a, z_a)=(x_a, 0, x_a\tan\theta)\tag{5}$$

となり、点$B$の座標は、

$$B(x_b, y_b, z_b)=(x_b, x_b\tan\delta, x_b\tan(\theta+\gamma))\tag{6}$$

と表すことができます。これらの座標を(4)の方程式に当てはめ、それぞれ解を求めると（結構しんどい連立方程式ですが、解けます）、

$$x_a=\frac{a(1-f)}{\sqrt{\tan^2\theta+(1-f)^2}}\tag{7}$$

$$y_a=0\tag{8}$$

$$z_a=\frac{a(1-f)\tan\theta}{\sqrt{\tan^2\theta+(1-f)^2}}\tag{9}$$

$$x_b=\frac{a(1-f)\cos\delta}{\sqrt{\cos^2\delta\tan^2(\theta+\gamma)+(1-f)^2}}\tag{10}$$

$$y_b=\frac{a(1-f)\sin\delta}{\sqrt{\cos^2\delta\tan^2(\theta+\gamma)+(1-f)^2}}\tag{11}$$

$$z_b=\frac{a(1-f)\cos\delta\tan(\theta+\gamma)}{\sqrt{\cos^2\delta\tan^2(\theta+\gamma)+(1-f)^2}}\tag{12}$$

となります。

これより、線分$AB$の座標を三平方の定理を使って求めると、

$$L=\sqrt{(x_a-x_b)^2+(y_a-y_b)^2+(z_a+z_b)^2}\tag{13}$$

として２点間の距離が求められます。

pythonの関数として作ると

def calculateLength(lat_a, lng_a, lat_b, lng_b):
    theta = lat_a / 180 * np.pi
    gamma = (lat_b - lat_a) / 180 * np.pi
    delta = (lng_b - lng_a) / 180 * np.pi
    xa = a * (1 - f) / np.sqrt(np.tan(theta)**2 + (1 - f)**2)
    ya = 0
    za = a * (1 - f) * np.tan(theta) / np.sqrt(np.tan(theta)**2 + (1 - f)**2)
    xb = a * (1 - f) * np.cos(delta) / np.sqrt(np.cos(delta)**2 * np.tan(theta + gamma)**2 + (1 - f)**2)
    yb = a * (1 - f) * np.sin(delta) / np.sqrt(np.cos(delta)**2 * np.tan(theta + gamma)**2 + (1 - f)**2)
    zb = a * (1 - f) * np.cos(delta) * np.tan(theta + gamma) / np.sqrt(np.cos(delta)**2 * np.tan(theta + gamma)**2 + (1 - f)**2)
    return np.sqrt((xb - xa)**2 + (yb - ya)**2 + (zb - za)**2)

このようになります。（引数は、点$A$、点$B$の緯度、経度です）

遅くはないですが、計算量が大きいですね。もう少し、しなくていい計算を減らしてあげたくなります。

次は、もう少し合理的な近似をして、簡単に求める方法を説明したいと思います。

CHAR-IN-CO.SE