# Copyright 2013-2016 Mike Bostock
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# * Neither the name of the author nor the names of contributors may be used to
# endorse or promote products derived from this software without specific prior
# written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# https://github.com/d3/d3-geo-projection
# commit fd2886555e46b35163c7b898d43c7d1bcebbba7c 2016-07-02
#
# 2017-11-01 Paul Kienzle
# * converted to python, with degrees rather than radians
"""
Convert between latitude-longitude and Guyou map coordinates.
"""
from __future__ import division, print_function
import numpy as np
from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees
from numpy import sinh, arctan as atan
# scipy version of special functions
from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF
_ = """
# mpmath version of special functions
import mpmath as mp
sn, cn, dn = [mp.ellipfun(v) for v in 'sn', 'cn', 'dn']
def ellipticJi(u, v, m):
z = u+v*1j
return sn(z, m), cn(z, m), dn(z, m)
def ellipticJ(u, m):
s, c, d = sn(u, m), cn(u, m), dn(u, m)
phi = mp.asin(s)
return s, c, d, phi
def ellipticFi(phi, psi, m):
z = phi + psi*1j
return mp.ellipf(z, m)
def ellipticF(phi, m):
return mp.ellipf(phi, m)
"""
[docs]def ellipticJi(u, v, m):
scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m)
u, v, m = np.broadcast_arrays(u, v, m)
result = np.empty_like([u, u, u], 'D')
real = (v == 0)
imag = (u == 0)
mixed = ~(real|imag)
result[:, real] = _ellipticJi_real(u[real], m[real])
result[:, imag] = _ellipticJi_imag(v[imag], m[imag])
result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed])
return result[0, :] if scalar else result
def _ellipticJi_real(u, m):
sn, cn, dn, phi = ellipticJ(u, m)
return sn, cn, dn
def _ellipticJi_imag(v, m):
sn, cn, dn, phi = ellipticJ(v, 1-m)
return 1j*sn/cn, 1/cn, dn/cn
def _ellipticJi(u, v, m):
# Ignoring special cases for now.
# u=0: (1j*b[0]/b[1], 1/b[1], b[2]/b[1])
# v=0: (a[0], a[1], a[2])
a = ellipticJ(u, m)
b = ellipticJ(v, 1 - m)
c = b[1]**2 + m * (a[0] * b[0])**2
return [
(a[0] * b[2] / c) + 1j*(a[1] * a[2] * b[0] * b[1] / c),
(a[1] * b[1] / c) + 1j*(-a[0] * a[2] * b[0] * b[2] / c),
(a[2] * b[1] * b[2] / c) + 1j*(-m * a[0] * a[1] * b[0] / c),
]
# calculate F(phi+ipsi|m).
# see Abramowitz and Stegun, 17.4.11.
[docs]def ellipticFi(phi, psi, m):
if np.any(phi == 0):
scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m)
phi, psi, m = np.broadcast_arrays(phi, psi, m)
result = np.empty_like(phi, 'D')
index = (phi == 0)
result[index] = ellipticF(atan(sinh(abs(phi[index]))),
1-m[index]) * sign(psi[index])
result[~index] = ellipticFi(phi[~index], psi[~index], m[~index])
return result.reshape(1)[0] if scalar else result
r = abs(phi)
i = abs(psi)
sinhpsi2 = sinh(i)**2
cscphi2 = 1 / sin(r)**2
cotphi2 = 1 / tan(r)**2
b = -(cotphi2 + m * (sinhpsi2 * cscphi2) - 1 + m)
c = (m - 1) * cotphi2
cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2
re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi)
im = ellipticF(atan(sqrt(np.maximum(0, (cotlambda2 / cotphi2 - 1) / m))),
1 - m) * sign(psi)
return re + 1j*im
SQRT2 = sqrt(2)
# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion
# cos_u = 0.171572875253809902396622551580603842860656249246103853646...
# sinsq_u = 0.970562748477140585620264690516376942836062504523376878120...
# K = 3.165103454447431823666270142140819753058976299237578486994...
[docs]def guyou(lam, phi):
"""Transform from (latitude, longitude) to point (x, y)"""
# [PAK] wrap into [-pi/2, pi/2] radians
x, y = np.asarray(lam), np.asarray(phi)
xn, x = divmod(x+90, 180)
yn, y = divmod(y+90, 180)
xn, lam = xn*180, radians(x-90)
yn, phi = yn*180, radians(y-90)
# Compute constant K
cos_u = (SQRT2 - 1) / (SQRT2 + 1)
sinsq_u = 1 - cos_u**2
K = ellipticF(pi/2, sinsq_u)
# [PAK] simplify expressions, using the fact that f = -1
# Note: exp(f log(x)) => 1/x, cos(f x) => cos(x), sin(f x) => -sin(x)
r = 1/(tan(pi/4 + abs(phi)/2) * sqrt(cos_u))
at = atan(r * (cos(lam) - 1j*sin(lam)))
t = ellipticFi(at.real, at.imag, sinsq_u)
x, y = (-t.imag, sign(phi + (phi == 0))*(0.5 * K - t.real))
# [PAK] convert to degrees, and return to original tile
return degrees(x)+xn, degrees(y)+yn
[docs]def guyou_invert(x, y):
"""Transform from point (x, y) on plot to (latitude, longitude)"""
# [PAK] wrap into [-pi/2, pi/2] radians
x, y = np.asarray(x), np.asarray(y)
xn, x = divmod(x+90, 180)
yn, y = divmod(y+90, 180)
xn, x = xn*180, radians(x-90)
yn, y = yn*180, radians(y-90)
# compute constant K
cos_u = (SQRT2 - 1) / (SQRT2 + 1)
sinsq_u = 1 - cos_u**2
K = ellipticF(pi/2, sinsq_u)
# [PAK] simplify expressions, using the fact that f = -1
j = ellipticJi(K/2 - y, -x, sinsq_u)
tn = j[0]/j[1] # j[0], j[1] are complex
# Note: -atan2(im(x)/re(x)) => angle(x)
lam = -np.angle(tn)
# Note: exp(0.5/f log(a re(x)^2 + a im(x)^2)) => 1/(sqrt(a) |x|)
phi = 2*atan(1/sqrt(cos_u)/abs(tn)) - pi/2
# [PAK] convert to degrees, and return to original tile
return degrees(lam)+xn, degrees(phi)+yn
[docs]def plot_grid():
"""Plot the latitude-longitude grid for Guyou transform"""
import matplotlib.pyplot as plt
from numpy import linspace
lat_line = linspace(-90, 90, 400)
long_line = linspace(-90, 90, 400)
#scale = 1
limit, step, scale = 90, 10, 2
plt.subplot(211)
for lat in range(-limit, limit+1, step):
x, y = guyou(scale*lat_line, scale*lat)
plt.plot(x, y, 'g')
for longitude in range(-limit, limit+1, step):
x, y = guyou(scale*longitude, scale*long_line)
plt.plot(x, y, 'b')
#plt.xlabel('longitude')
plt.ylabel('latitude')
plt.title('forward transform')
plt.subplot(212)
for lat in range(-limit, limit+1, step):
x, y = guyou_invert(scale*lat_line, scale*lat)
plt.plot(x, y, 'g')
for long in range(-limit, limit+1, step):
x, y = guyou_invert(scale*long, scale*long_line)
plt.plot(x, y, 'b')
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.title('inverse transform')
[docs]def main():
"""Show the Guyou transformation"""
plot_grid()
import matplotlib.pyplot as plt
plt.show()
if __name__ == "__main__":
main()
_ = """
// Javascript source for elliptic functions
//
// Returns [sn, cn, dn](u + iv|m).
export function ellipticJi(u, v, m) {
var a, b, c;
if (!u) {
b = ellipticJ(v, 1 - m);
return [
[0, b[0] / b[1]],
[1 / b[1], 0],
[b[2] / b[1], 0]
];
}
a = ellipticJ(u, m);
if (!v) return [[a[0], 0], [a[1], 0], [a[2], 0]];
b = ellipticJ(v, 1 - m);
c = b[1] * b[1] + m * a[0] * a[0] * b[0] * b[0];
return [
[a[0] * b[2] / c, a[1] * a[2] * b[0] * b[1] / c],
[a[1] * b[1] / c, -a[0] * a[2] * b[0] * b[2] / c],
[a[2] * b[1] * b[2] / c, -m * a[0] * a[1] * b[0] / c]
];
}
// Returns [sn, cn, dn, ph](u|m).
export function ellipticJ(u, m) {
var ai, b, phi, t, twon;
if (m < epsilon) {
t = sin(u);
b = cos(u);
ai = m * (u - t * b) / 4;
return [
t - ai * b,
b + ai * t,
1 - m * t * t / 2,
u - ai
];
}
if (m >= 1 - epsilon) {
ai = (1 - m) / 4;
b = cosh(u);
t = tanh(u);
phi = 1 / b;
twon = b * sinh(u);
return [
t + ai * (twon - u) / (b * b),
phi - ai * t * phi * (twon - u),
phi + ai * t * phi * (twon + u),
2 * atan(exp(u)) - halfPi + ai * (twon - u) / b
];
}
var a = [1, 0, 0, 0, 0, 0, 0, 0, 0],
c = [sqrt(m), 0, 0, 0, 0, 0, 0, 0, 0],
i = 0;
b = sqrt(1 - m);
twon = 1;
while (abs(c[i] / a[i]) > epsilon && i < 8) {
ai = a[i++];
c[i] = (ai - b) / 2;
a[i] = (ai + b) / 2;
b = sqrt(ai * b);
twon *= 2;
}
phi = twon * a[i] * u;
do {
t = c[i] * sin(b = phi) / a[i];
phi = (asin(t) + phi) / 2;
} while (--i);
// [PAK] Cephes uses dn = sqrt(1 - m*sin^2 phi) rather than cos(phi)/cos(phi-b)
// DLMF says the second version is unstable near x = K.
return [sin(phi), t = cos(phi), t / cos(phi - b), phi];
}
// calculate F(phi+ipsi|m).
// see Abramowitz and Stegun, 17.4.11.
export function ellipticFi(phi, psi, m) {
var r = abs(phi),
i = abs(psi),
sinhpsi = sinh(i);
if (r) {
var cscphi = 1 / sin(r),
cotphi2 = 1 / (tan(r) * tan(r)),
b = -(cotphi2 + m * (sinhpsi * sinhpsi * cscphi * cscphi) - 1 + m),
c = (m - 1) * cotphi2,
cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2;
return [
ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi),
ellipticF(atan(sqrt((cotlambda2 / cotphi2 - 1) / m)), 1 - m) * sign(psi)
];
}
return [
0,
ellipticF(atan(sinhpsi), 1 - m) * sign(psi)
];
}
// Calculate F(phi|m) where m = k^2 = sin(alpha)^2.
// See Abramowitz and Stegun, 17.6.7.
export function ellipticF(phi, m) {
if (!m) return phi;
if (m === 1) return log(tan(phi / 2 + quarterPi));
var a = 1,
b = sqrt(1 - m),
c = sqrt(m);
for (var i = 0; abs(c) > epsilon; i++) {
if (phi % pi) {
var dPhi = atan(b * tan(phi) / a);
if (dPhi < 0) dPhi += pi;
phi += dPhi + ~~(phi / pi) * pi;
} else phi += phi;
c = (a + b) / 2;
b = sqrt(a * b);
c = ((a = c) - b) / 2;
}
return phi / (pow(2, i) * a);
export function guyouRaw(lambda, phi) {
var k_ = (sqrt2 - 1) / (sqrt2 + 1),
k = sqrt(1 - k_ * k_),
K = ellipticF(halfPi, k * k),
f = -1,
psi = log(tan(pi / 4 + abs(phi) / 2)),
r = exp(f * psi) / sqrt(k_),
at = guyouComplexAtan(r * cos(f * lambda), r * sin(f * lambda)),
t = ellipticFi(at[0], at[1], k * k);
return [-t[1], (phi >= 0 ? 1 : -1) * (0.5 * K - t[0])];
}
function guyouComplexAtan(x, y) {
var x2 = x * x,
y_1 = y + 1,
t = 1 - x2 - y * y;
return [
0.5 * ((x >= 0 ? halfPi : -halfPi) - atan2(t, 2 * x)),
-0.25 * log(t * t + 4 * x2) +0.5 * log(y_1 * y_1 + x2)
];
}
function guyouComplexDivide(a, b) {
var denominator = b[0] * b[0] + b[1] * b[1];
return [
(a[0] * b[0] + a[1] * b[1]) / denominator,
(a[1] * b[0] - a[0] * b[1]) / denominator
];
}
guyouRaw.invert = function(x, y) {
var k_ = (sqrt2 - 1) / (sqrt2 + 1),
k = sqrt(1 - k_ * k_),
K = ellipticF(halfPi, k * k),
f = -1,
j = ellipticJi(0.5 * K - y, -x, k * k),
tn = guyouComplexDivide(j[0], j[1]),
lambda = atan2(tn[1], tn[0]) / f;
return [
lambda,
2 * atan(exp(0.5 / f * log(k_ * tn[0] * tn[0] + k_ * tn[1] * tn[1]))) - halfPi
];
};
"""