メルカトル投影計算（回転楕円体モデル）

こんにちは
地図の[メルカトル投影](https://en.wikipedia.org/wiki/Mercator_projection#Generalization_to_the_ellipsoid)を計算しました。下記間の変換です（回転楕円体モデル[^1]）。Python 利用です。

* EPSG:4326  WGS 84 geographical:  $(\lambda, \varphi)$
* EPSG:3395  WGS 84 World Mercator  (ellipsoidal earth):  $(x, y)$

[^1]: 真球モデルの場合は、「[webメルカトル投影計算（web地図タイル）](https://qiita.com/kkdd/items/e3fa99aec45f87aca0f4#_reference-7a34090d55e21e64997a)」。

### 緯度

緯度 $\varphi$ ⇔ $y$ の変換式は（離心率 $e$、地球長半径 $R$）：

```math
\begin{align}
\frac{2 y}{R} &= \ln \left( \frac {1 + \sin\varphi}{1 - \sin\varphi} \right) - e \ln \left( \frac {1 + e\sin\varphi}{1 - e\sin\varphi} \right) \\
&= 2 \left( \tanh^{-1} \sin \varphi - e \tanh^{-1} \left(e \sin \varphi \right) \right)
\end{align}
```

逆変換は下記の変形式を使い、[Newton 法](https://en.wikipedia.org/wiki/Newton%27s_method)で $\sin \varphi$ を求めています（$\varphi$ を求めるには、それから $\varphi = \arcsin \(\sin \varphi\)$ を計算します[^2]）。

```math
\sin \varphi - 1 + \frac{2}{\left( \frac{1+e\sin\varphi}{1-e\sin\varphi}\right)^e \tau^2 +1} = 0,\ \\
\tau \triangleq \exp\left(\frac{y}{R}\right),\ \\
\tau^2=\exp\left(\frac{2 y}{R}\right)
```

その初期値としては、解の展開形（下記）を $e^2$ の項で打ち切って使えば、反復は１回で十分です。

```math
\sin \varphi = \frac{\tau^2-1}{\tau^2+1} + \frac{4 e^2 \tau^2 \left(\tau^2-1\right)}{\left(\tau^2+1\right)^3} + \mathcal{O} (e^4)
```

[^2]: ただし極付近では $\arcsin$ の計算が精度劣化を増大させます。

#### Braun投影

なお上記を[変数変換](http://qiita.com/kkdd/items/864aa75041b77e2b4c28)（Braun投影）すると、

```math
\chi \triangleq \tan \frac{\varphi}{2}\ ,\\
\sin \varphi = \frac{2 \chi}{1+\chi^2}\ ,\\
\frac{y}{R} = \ln \tau = \ln \left( \frac {1 + \chi}{1 - \chi} \right) - \frac{e}{2} \ln \left( \frac{1+\chi^2+2e\chi}{1+\chi^2-2e\chi} \right)
```

逆変換に用いる変形式は、

```math
\chi - 1 + \frac{2}{\left( \frac{1+\chi^2+2e\chi}{1+\chi^2-2e\chi} \right)^{\frac{e}{2}} \tau + 1} = 0,\ \\
\chi = \frac{\tau-1}{\tau+1} + \frac{2 e^2 \tau \left(\tau-1\right) }{ \left(\tau^2+1\right) \left(\tau+1\right) } + \mathcal{O} (e^4)
```

### 経度

経度は、

```math
\lambda = \frac{x}{R}
```

### スケール係数

スケール係数は

```math
k = \sqrt{\frac{1-e^2 \sin^2 \varphi}{1-\sin^2 \varphi}}
```

### 実行例・ソース

```console
$ ./mercator_lat.py
lat=  0.0001: error=  1.00e-16
lat= 20.0000: error= -1.11e-16
lat= 40.0000: error=  1.11e-16
lat= 60.0000: error=  0.00e+00
lat= 80.0000: error= -2.22e-16
lat= 89.0000: error=  1.78e-15
lat= 89.9000: error= -3.55e-15
lat= 89.9900: error=  3.03e-13
lat= 89.9990: error= -6.47e-13
lat= 89.9999: error= -1.60e-11
```

```mercator_lat.py
#!/usr/bin/env python

# ellipsoidal Mercator projection and its reverse (inverse)
#  latitude: EPSG:4326 WGS 84 geographical
#  y: EPSG:3395 WGS 84 World Mercator (ellipsoidal earth)

from __future__ import print_function
import numpy as np

DEG = np.pi / 180
RE = 6378137.0  # GRS80
FE = (1/298.257223563) # IS-GPS
E2 = FE * (2 - FE)
E1 = np.sqrt(E2)

def mercator(lat):
    slat = np.sin(lat)
    y2 = np.log((1+slat)/(1-slat)) - E1*np.log((1+E1*slat)/(1-E1*slat))
    return y2/2

def mercator_inverse(y):
    niter = 1
    ey2 = np.exp(y*2)    
    slat = 1 - 2/(1+ey2) # initial value (conformal latitude)
    slat += 4*E2*ey2*(ey2-1)/(ey2+1)**3 # E2-term
    for i in range(niter): # Newton's method
        w = ey2*((1+E1*slat)/(1-E1*slat))**E1
        dslat = slat - 1 + 2/(1+w)
        slat -= dslat / (1 - 4*E2*w/((1-E2*slat**2)*(1+w)**2))
    return np.arcsin(slat)

def mercator_scale(slat):
    return np.sqrt((1-E2*slat**2)/(1-slat**2))

for lat in [0.0001, 20.0, 40.0, 60.0, 80.0, 89.0, 89.9, 89.99, 89.999, 89.9999]:
    print("lat= %7.4f: error= %9.2e" % (lat, mercator_inverse(mercator(lat*DEG)) - lat*DEG))
exit(0)
```