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

こんにちは
地図のメルカトル投影を Python で計算しました。下記間の変換です（回転楕円体モデル¹）。

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

緯度

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

\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) 

逆変換は下記の変形式を使い、Newton 法で $\sin \varphi$ を求めています（$\varphi$ を求めるには、それから $\varphi = \arcsin (\sin \varphi)$ を計算します²）。

\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$ の項で打ち切って使えば、反復は１回で十分です。

\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)

Braun投影

なお上記を変数変換（Braun投影）すると、

\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)

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

\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)

経度

経度は、

\lambda = \frac{x}{R}

スケール係数

スケール係数は

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

実行例・ソース

$ ./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)

1. 真球モデルの場合は、「webメルカトル投影計算（web地図タイル）」。 ↩

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