sw1227’s diary

Visualization, GIS, Machine Learning, Generative Art, Simulation, Math

緯度経度と平面直角座標の相互変換を実装するための数式

GIS Math Visualization

1. はじめに

1.1. 背景

緯度経度はご存知のように地球上の位置を角度で表現するもので、Google Mapsなど様々なサービスで目にすることも多いでしょう。 しかしながら、角度ではなくメートル単位の x, y座標で地球上の位置を表現したいケースも存在します。そんな時に使えるのが平面直角座標系です。

経緯度と平面直角座標の相互変換式は国土地理院のWebサイトに公開されているのですが、角度の「ラジアン」と「度」に混同があったり変数間の依存関係が見辛かったりする*1ため、この記事で解説したいと思います。 なお、Pythonによる実装は別途Qiitaに投稿しています。

緯度経度と平面直角座標の相互変換をPythonで実装する - Qiita

1.2. 平面直角座標系の概要

平面直角座標系の大まかなイメージとしては、ある地点を原点とし、その近辺で地球が平坦であるとみなしたものであると考えれば良いでしょう。もちろん地球の楕円体を歪みなく平面に変換することは不可能ですが、なるべく「ずれ」が無いように工夫してあります。

ただ、日本全土のスケールになってくると地球の丸みは無視できないものになるため、「ずれ」を抑えるためには原点を複数点用意する必要が生じます。実際、日本国内では平成十四年国土交通省告示第九号によって19の平面直角座標系が定められており、それぞれの座標系原点の経緯度と適用範囲が定められています。平面直角座標はこの適用範囲それぞれの内部においてのみ比較可能なものであることに留意しましょう。*2

なお、適用範囲に関してはわかりやすい平面直角座標系に詳しい記載があります。

2. 平面直角座標から経緯度への変換式

基本的には国土地理院の公開している計算式に基づきます。 ただ、冒頭で述べたようにこの式では角度の表現に混同がみられたため、以下の式では全てラジアンに統一するものとします。すなわち、経緯度はラジアンで表現されており、三角関数はラジアンを引数として受け取り、逆三角関数はラジアンを返すものであるとします。

2.1. 記号

準拠楕円体に関しては以下を参照します。

2001年以前の測地基準点成果は、緯度・経度においては日本測地系に基づいた数値で準拠楕円体はベッセル楕円体でしたが、現在の測地測量成果は世界測地系（測地成果2011）と呼び、準拠楕円体はITRF座標系GRS80楕円体です。
日本測地系と世界測地系

2.2. 計算詳細

目標は、パラメータ \phi_0, \lambda_0, a, F, m_0 を既知として、平面直角座標 x, yを緯度 \phiと経度 \lambdaに変換することです。

しかし、国土地理院の計算式を見ると計算過程に中間的な変数が多数存在します。それぞれの依存関係がよく分からず、「この変数を計算するためにはこっちの変数を計算しておく必要があって、その前に…」と途方に暮れてしまうかもしれません。

そこで、まずは変数間の依存関係を可視化してみます*3。

入力 x, y（平面直角座標）が出力 \phi, \lambda（緯度経度）に変換される過程を、国土地理院の計算式の記号に合わせてグラフ状に表現したものが下の画像です。入力 x, y を青背景・出力 \phi, \lambdaを赤背景・定数（地球の形状などに依存）を黒背景のノードとして塗り分けています。

変数間の依存関係
変数間の依存関係

この図を見ながら、順番に計算式を追っていきます。

2.2.1. n, A_j, \beta_j, \delta_j の計算

先ほどの図を見ると、以下の赤線のように n, A_j, \beta_j, \delta_j は F （逆扁平率）だけに依存していることが分かります。

変数間の依存関係(2)
変数間の依存関係(2)

具体的な数式に移りましょう。まずは単純に nを計算します。

\displaystyle { n = \frac{1}{2F-1} }

この nから以下のように A_i(i=0, 1, ..., 5), \; \beta _i(i=1, 2, ..., 5), \; \delta_i(i=1, 2, ..., 6) が定まります。

  • A_i(i=0, 1, ..., 5)
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    A_0 =  1 + \frac{n^2}{4} + \frac{n^4}{64} \\
    A_1 =  -\frac{3}{2} \left( n-\frac{n^3}{8}-\frac{n^5}{64}  \right) \\
    A_2 =  \frac{15}{16} \left( n^2-\frac{n^4}{4}  \right) \\
    A_3 =  -\frac{35}{48} \left( n^3-\frac{5}{16}n^5  \right) \\
    A_4 =  \frac{315}{512} n^4 \\
    A_5 =  -\frac{693}{1280}n^5
    \end{array}
  \right.
\end{eqnarray}
}
  • \beta _i(i=1, 2, ..., 5)
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    \beta_1 =  \frac{1}{2}n - \frac{2}{3}n^2 + \frac{37}{96}n^3 -\frac{1}{360}n^4 - \frac{81}{512}n^5 & \\
    \beta_2 =  \frac{1}{48}n^2 + \frac{1}{15}n^3 - \frac{437}{1440}n^4 + \frac{46}{105}n^5 & \\
    \beta_3 =  \frac{17}{480}n^3 - \frac{37}{840}n^4 - \frac{209}{4480}n^5 & \\
    \beta_4 =  \frac{4397}{161280}n^4 - \frac{11}{504}n^5 & \\
    \beta_5 =  \frac{4583}{161280}n^5 &
    \end{array}
  \right.
\end{eqnarray}
}
  • \delta_i(i=1, 2, ..., 6)
 \displaystyle {
\begin{eqnarray}
  \begin{cases}
    \delta_1 =  2n - \frac{2}{3}n^2 - 2n^3 + \frac{116}{45}n^4 + \frac{26}{45}n^5 - \frac{2854}{675}n^6 & \\
    \delta_2 =  \frac{7}{3}n^2 - \frac{8}{5}n^3 - \frac{227}{45}n^4 + \frac{2704}{315}n^5 + \frac{2323}{945}n^6 & \\
    \delta_3 =  \frac{56}{15}n^3 - \frac{136}{35}n^4 - \frac{1262}{105}n^5 + \frac{73814}{2835}n^6 & \\
    \delta_4 =  \frac{4279}{630}n^4 - \frac{332}{35}n^5 - \frac{399572}{14175}n^6 & \\
    \delta_5 =  \frac{4174}{315}n^5 - \frac{144838}{6237}n^6 & \\
    \delta_6 =  \frac{601676}{22275}n^6   &
  \end{cases}
\end{eqnarray}
}

2.2.2. \overline{S}_{\phi_0}, \overline{A} の計算

すでに計算できた n, A_j, \beta_j, \delta_j は定数と同じ黒背景に塗りつぶしてしまいます。すると、今度は下図のように \overline{S}_{\phi_0}, \overline{A} が計算できますね。

変数間の依存関係(3)
変数間の依存関係(3)

数式は以下のようになります。 この記事では \phi_0 をラジアンとしているため、国土地理院の計算式の \rho^{"} を除いていることに注意しましょう。

\displaystyle { \overline{S}_{\phi_0} = \frac{m_0 a}{1+n} \left( A_0 \phi_0 + \sum_{j=1}^5 A_j sin(2j\phi_0) \right) }

\displaystyle { \overline{A} = \frac{m_0 a}{1+n} A_0 }

2.2.3. \xi, \eta の計算

今まで通り計算済みの変数を黒背景に塗りつぶすと、計算済みの変数と入力 x, yから \xi, \eta が求まります。

 \displaystyle {
\begin{eqnarray}
  \begin{cases}
    \xi = \frac{x + \overline{S}_{\phi_0}}{\overline{A}} & \\
    \eta = \frac{y}{\overline{A}} & 
  \end{cases}
\end{eqnarray}
}

変数間の依存関係(4)
変数間の依存関係(4)

2.2.4. \xi^{'}, \eta^{'} の計算

同様に、 \xi^{'}, \eta^{'} が計算できます。

 \displaystyle {
\begin{eqnarray}
  \begin{cases}
    \xi^{'} = \xi - \sum_{j=1}^5 \beta_j sin(2j\xi) cosh(2j\eta) & \\
    \eta^{'} = \eta - \sum_{j=1}^5 \beta_j cos(2j\xi) sinh(2j\eta) &
  \end{cases}
\end{eqnarray}
}

変数間の依存関係(5)
変数間の依存関係(5)

2.2.5. \chi の計算

\chi を計算します（単位はrad）。あともう少し。

\displaystyle { \chi = sin^{-1} \left( \frac{sin \; \xi^{'}}{cosh \; \eta{'}} \right) \;\; [rad] }

変数間の依存関係(6)
変数間の依存関係(6)

2.2.6. 緯度 \phi ・経度 \lambda の計算

目的だった緯度 \phi ・経度 \lambda を計算する準備が整いました。 本記事では \phi, \lambda をラジアンに統一している関係上、ここでも国土地理院の計算式の \rho^{"}を除いていることに注意してください。

  • 緯度

\displaystyle { \phi = \chi + \sum_{j=1}^6 \delta_j sin(2j\chi) \;\; [rad] }

  • 経度

\displaystyle { \lambda = \lambda_0 + tan^{-1} \left( \frac{sinh \; \eta^{'}}{cos \; \xi{'}} \right) \;\; [rad] }

変数間の依存関係(7)
変数間の依存関係(7)

依存関係を整理することによって順序立てて計算を進めることができました。

3. 経緯度から平面直角座標への変換式

先ほどまでとは逆に、経緯度から平面直角座標への変換式を考えます。

こちらも基本的に国土地理院の公開している計算式に基づきます。 以下の注意点はこの章にもあてはまります。

この式では角度の表現に混同がみられたため、以下の式では全てラジアンに統一するものとします。すなわち、経緯度はラジアンで表現されており、三角関数はラジアンを引数として受け取り、逆三角関数はラジアンを返すものであるとします。

　　

この記事では \phi_0 をラジアンとしているため、国土地理院の計算式の \rho^{"} を除いていることに注意しましょう。

3.1. 記号

「2.1. 記号」と全く同じです。

3.2. 計算詳細

目標は、 \phi_0, \lambda_0, a, F, m_0 を既知として、緯度 \phiと経度 \lambdaを平面直角座標 x, yに変換することです。

逆変換に関しても、国土地理院の式における変数間の依存関係を可視化してみましょう。

変数間の依存関係
変数間の依存関係

あとは、平面直角座標から経緯度への変換式の時と同様に順番に求めていきます。

  • n
 \displaystyle {
n = \frac{1}{2F-1}
}
  • A_i(i=0, 1, ..., 5)
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    A_0 =  1 + \frac{n^2}{4} + \frac{n^4}{64} \\
    A_1 =  -\frac{3}{2} \left( n-\frac{n^3}{8}-\frac{n^5}{64}  \right) \\
    A_2 =  \frac{15}{16} \left( n^2-\frac{n^4}{4}  \right) \\
    A_3 =  -\frac{35}{48} \left( n^3-\frac{5}{16}n^5  \right) \\
    A_4 =  \frac{315}{512} n^4 \\
    A_5 =  -\frac{693}{1280}n^5
    \end{array}
  \right.
\end{eqnarray}
}
  • \alpha_i(i=1, 2, ..., 5)
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
     \alpha_1 = \frac{1}{2} n - \frac{2}{3} n^2 +  \frac{5}{16} n^3  + \frac{41}{180} n^4 -  \frac{127}{288} n^5   \\
     \alpha_2 = \frac{13}{48} n^2 - \frac{3}{5} n^3  + \frac{557}{1440} n^4 +  \frac{281}{630} n^5   \\
     \alpha_3 = \frac{61}{240} n^3  - \frac{103}{140} n^4 +  \frac{15061}{26880} n^5   \\
     \alpha_4 = \frac{49561}{161280} n^4 - \frac{179}{168} n^5   \\
     \alpha_5 = \frac{34729}{80640} n^5
    \end{array}
  \right.
\end{eqnarray}
}
  • \overline{S}_{\phi_0}, \overline{A}
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
\overline{S}_{\phi_0} = \frac{m_0 a}{1+n} \left( A_0 \phi_0 + \sum_{j=1}^5 A_j sin(2j\phi_0) \right) \\
\overline{A} = \frac{m_0 a}{1+n} A_0
    \end{array}
  \right.
\end{eqnarray}
}
  • \lambda_c, \lambda_s
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    \lambda_c = cos( \lambda - \lambda_0 ) \\
    \lambda_s = sin( \lambda - \lambda_0 )
    \end{array}
  \right.
\end{eqnarray}
}
  • t, \overline{t}
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    t = sinh \left( tanh^{-1} sin \phi - \frac{2 \sqrt{n}}{1+n} tanh^{-1} \left[  \frac{2 \sqrt{n}}{1+n}  sin \phi \right] \right) \\
    \overline{t} = \sqrt{1+t^{2}}
    \end{array}
  \right.
\end{eqnarray}
}
  • \xi^{'}, \eta^{'}
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    \xi^{'} = tan^{-1} \left( \frac{t}{\lambda_c} \right) \\
    \eta^{'} = tan^{-1} \left( \frac{\lambda_s}{\overline{t}} \right)
    \end{array}
  \right.
\end{eqnarray}
}
  • x, y
 \displaystyle {
\begin{eqnarray}
  \left\{
    \begin{array}{l}
    x = \overline{A} \left( \xi^{'} + \sum_{j=1}^{5} \alpha_j sin (2j \xi^{'}) cosh(2j \eta^{'}) \right) -  \overline{S}_{\phi_0} \\
    y =  \overline{A}  \left( \eta^{'} +  \sum_{j=1}^{5}  \alpha_j cos (2j \xi^{'}) sinh(2j \eta^{'}) \right) 
    \end{array}
  \right.
\end{eqnarray}
}

4. 実装

この記事に記した順序でプログラムに書き起こすだけです。 Qiitaで公開しています。 qiita.com

5. まとめ

  • 緯度経度と平面直角座標系の相互変換式を紹介した
  • 計算式と記号は国土地理院のWebサイトに基づいているが、「ラジアン」と「度」の混同を無くす修正を行った
  • 変数間の依存関係を可視化して整理することにより、そのままプログラムに書き起こせば実装できるような順序で数式を展開した

6. 参考

以上

地図投影法

地図投影法

*1:変換式実装時点

*2:誤差が問題にならないなら適用範囲を無視してもいいかもしれませんが…

*3:Deep LearningにおけるComputational Graphみたいなイメージ