Skip to content

Instantly share code, notes, and snippets.

@Fil
Last active Apr 30, 2017
Embed
What would you like to do?
InterruptedSinuMollweide inverse projection
license: gpl-3.0
border: no
height: 484

Computing the inverse of the InterruptedSinuMollweide projection for issue #85.

See Fil's block: ChamberlinAfrica inverse projection for details of the method.

For this specific projection we had to add syntactic sugar, and export d3-geo:solve as d3.geoSolve() to be able to use it directly in d3-geo-projection:interrupt.invert.

// https://d3js.org/d3-geo-projection/ Version 2.1.0. Copyright 2017 Mike Bostock.
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-geo'), require('d3-array')) :
typeof define === 'function' && define.amd ? define(['exports', 'd3-geo', 'd3-array'], factory) :
(factory((global.d3 = global.d3 || {}),global.d3,global.d3));
}(this, (function (exports,d3Geo,d3Array) { 'use strict';
var abs = Math.abs;
var atan = Math.atan;
var atan2 = Math.atan2;
var cos = Math.cos;
var exp = Math.exp;
var floor = Math.floor;
var log = Math.log;
var max = Math.max;
var min = Math.min;
var pow = Math.pow;
var round = Math.round;
var sign = Math.sign || function(x) { return x > 0 ? 1 : x < 0 ? -1 : 0; };
var sin = Math.sin;
var tan = Math.tan;
var epsilon = 1e-6;
var epsilon2 = 1e-12;
var pi = Math.PI;
var halfPi = pi / 2;
var quarterPi = pi / 4;
var sqrt1_2 = Math.SQRT1_2;
var sqrt2 = sqrt(2);
var sqrtPi = sqrt(pi);
var tau = pi * 2;
var degrees = 180 / pi;
var radians = pi / 180;
function sinci(x) {
return x ? x / Math.sin(x) : 1;
}
function asin(x) {
return x > 1 ? halfPi : x < -1 ? -halfPi : Math.asin(x);
}
function acos(x) {
return x > 1 ? 0 : x < -1 ? pi : Math.acos(x);
}
function sqrt(x) {
return x > 0 ? Math.sqrt(x) : 0;
}
function tanh(x) {
x = exp(2 * x);
return (x - 1) / (x + 1);
}
function sinh(x) {
return (exp(x) - exp(-x)) / 2;
}
function cosh(x) {
return (exp(x) + exp(-x)) / 2;
}
function arsinh(x) {
return log(x + sqrt(x * x + 1));
}
function arcosh(x) {
return log(x + sqrt(x * x - 1));
}
function airyRaw(beta) {
var tanBeta_2 = tan(beta / 2),
b = 2 * log(cos(beta / 2)) / (tanBeta_2 * tanBeta_2);
function forward(x, y) {
var cosx = cos(x),
cosy = cos(y),
siny = sin(y),
cosz = cosy * cosx,
k = -((1 - cosz ? log((1 + cosz) / 2) / (1 - cosz) : -0.5) + b / (1 + cosz));
return [k * cosy * sin(x), k * siny];
}
forward.invert = function(x, y) {
var r = sqrt(x * x + y * y),
z = -beta / 2,
i = 50, delta;
if (!r) return [0, 0];
do {
var z_2 = z / 2,
cosz_2 = cos(z_2),
sinz_2 = sin(z_2),
tanz_2 = tan(z_2),
lnsecz_2 = log(1 / cosz_2);
z -= delta = (2 / tanz_2 * lnsecz_2 - b * tanz_2 - r) / (-lnsecz_2 / (sinz_2 * sinz_2) + 1 - b / (2 * cosz_2 * cosz_2));
} while (abs(delta) > epsilon && --i > 0);
var sinz = sin(z);
return [atan2(x * sinz, r * cos(z)), asin(y * sinz / r)];
};
return forward;
}
var airy = function() {
var beta = halfPi,
m = d3Geo.geoProjectionMutator(airyRaw),
p = m(beta);
p.radius = function(_) {
return arguments.length ? m(beta = _ * radians) : beta * degrees;
};
return p
.scale(179.976)
.clipAngle(147);
};
function aitoffRaw(x, y) {
var cosy = cos(y), sincia = sinci(acos(cosy * cos(x /= 2)));
return [2 * cosy * sin(x) * sincia, sin(y) * sincia];
}
// Abort if [x, y] is not within an ellipse centered at [0, 0] with
// semi-major axis pi and semi-minor axis pi/2.
aitoffRaw.invert = function(x, y) {
if (x * x + 4 * y * y > pi * pi + epsilon) return;
var x1 = x, y1 = y, i = 25;
do {
var sinx = sin(x1),
sinx_2 = sin(x1 / 2),
cosx_2 = cos(x1 / 2),
siny = sin(y1),
cosy = cos(y1),
sin_2y = sin(2 * y1),
sin2y = siny * siny,
cos2y = cosy * cosy,
sin2x_2 = sinx_2 * sinx_2,
c = 1 - cos2y * cosx_2 * cosx_2,
e = c ? acos(cosy * cosx_2) * sqrt(f = 1 / c) : f = 0,
f,
fx = 2 * e * cosy * sinx_2 - x,
fy = e * siny - y,
dxdx = f * (cos2y * sin2x_2 + e * cosy * cosx_2 * sin2y),
dxdy = f * (0.5 * sinx * sin_2y - e * 2 * siny * sinx_2),
dydx = f * 0.25 * (sin_2y * sinx_2 - e * siny * cos2y * sinx),
dydy = f * (sin2y * cosx_2 + e * sin2x_2 * cosy),
z = dxdy * dydx - dydy * dxdx;
if (!z) break;
var dx = (fy * dxdy - fx * dydy) / z,
dy = (fx * dydx - fy * dxdx) / z;
x1 -= dx, y1 -= dy;
} while ((abs(dx) > epsilon || abs(dy) > epsilon) && --i > 0);
return [x1, y1];
};
var aitoff = function() {
return d3Geo.geoProjection(aitoffRaw)
.scale(152.63);
};
function armadilloRaw(phi0) {
var sinPhi0 = sin(phi0),
cosPhi0 = cos(phi0),
sPhi0 = phi0 >= 0 ? 1 : -1,
tanPhi0 = tan(sPhi0 * phi0),
k = (1 + sinPhi0 - cosPhi0) / 2;
function forward(lambda, phi) {
var cosPhi = cos(phi),
cosLambda = cos(lambda /= 2);
return [
(1 + cosPhi) * sin(lambda),
(sPhi0 * phi > -atan2(cosLambda, tanPhi0) - 1e-3 ? 0 : -sPhi0 * 10) + k + sin(phi) * cosPhi0 - (1 + cosPhi) * sinPhi0 * cosLambda // TODO D3 core should allow null or [NaN, NaN] to be returned.
];
}
forward.invert = function(x, y) {
var lambda = 0,
phi = 0,
i = 50;
do {
var cosLambda = cos(lambda),
sinLambda = sin(lambda),
cosPhi = cos(phi),
sinPhi = sin(phi),
A = 1 + cosPhi,
fx = A * sinLambda - x,
fy = k + sinPhi * cosPhi0 - A * sinPhi0 * cosLambda - y,
dxdLambda = A * cosLambda / 2,
dxdPhi = -sinLambda * sinPhi,
dydLambda = sinPhi0 * A * sinLambda / 2,
dydPhi = cosPhi0 * cosPhi + sinPhi0 * cosLambda * sinPhi,
denominator = dxdPhi * dydLambda - dydPhi * dxdLambda,
dLambda = (fy * dxdPhi - fx * dydPhi) / denominator / 2,
dPhi = (fx * dydLambda - fy * dxdLambda) / denominator;
lambda -= dLambda, phi -= dPhi;
} while ((abs(dLambda) > epsilon || abs(dPhi) > epsilon) && --i > 0);
return sPhi0 * phi > -atan2(cos(lambda), tanPhi0) - 1e-3 ? [lambda * 2, phi] : null;
};
return forward;
}
var armadillo = function() {
var phi0 = 20 * radians,
sPhi0 = phi0 >= 0 ? 1 : -1,
tanPhi0 = tan(sPhi0 * phi0),
m = d3Geo.geoProjectionMutator(armadilloRaw),
p = m(phi0),
stream_ = p.stream;
p.parallel = function(_) {
if (!arguments.length) return phi0 * degrees;
tanPhi0 = tan((sPhi0 = (phi0 = _ * radians) >= 0 ? 1 : -1) * phi0);
return m(phi0);
};
p.stream = function(stream) {
var rotate = p.rotate(),
rotateStream = stream_(stream),
sphereStream = (p.rotate([0, 0]), stream_(stream));
p.rotate(rotate);
rotateStream.sphere = function() {
sphereStream.polygonStart(), sphereStream.lineStart();
for (var lambda = sPhi0 * -180; sPhi0 * lambda < 180; lambda += sPhi0 * 90) sphereStream.point(lambda, sPhi0 * 90);
while (sPhi0 * (lambda -= phi0) >= -180) { // TODO precision?
sphereStream.point(lambda, sPhi0 * -atan2(cos(lambda * radians / 2), tanPhi0) * degrees);
}
sphereStream.lineEnd(), sphereStream.polygonEnd();
};
return rotateStream;
};
return p
.scale(218.695)
.center([0, 28.0974]);
};
function augustRaw(lambda, phi) {
var tanPhi = tan(phi / 2),
k = sqrt(1 - tanPhi * tanPhi),
c = 1 + k * cos(lambda /= 2),
x = sin(lambda) * k / c,
y = tanPhi / c,
x2 = x * x,
y2 = y * y;
return [
4 / 3 * x * (3 + x2 - 3 * y2),
4 / 3 * y * (3 + 3 * x2 - y2)
];
}
augustRaw.invert = function(x, y) {
x *= 3 / 8, y *= 3 / 8;
if (!x && abs(y) > 1) return null;
var x2 = x * x,
y2 = y * y,
s = 1 + x2 + y2,
sin3Eta = sqrt((s - sqrt(s * s - 4 * y * y)) / 2),
eta = asin(sin3Eta) / 3,
xi = sin3Eta ? arcosh(abs(y / sin3Eta)) / 3 : arsinh(abs(x)) / 3,
cosEta = cos(eta),
coshXi = cosh(xi),
d = coshXi * coshXi - cosEta * cosEta;
return [
sign(x) * 2 * atan2(sinh(xi) * cosEta, 0.25 - d),
sign(y) * 2 * atan2(coshXi * sin(eta), 0.25 + d)
];
};
var august = function() {
return d3Geo.geoProjection(augustRaw)
.scale(66.1603);
};
var sqrt8 = sqrt(8);
var phi0 = log(1 + sqrt2);
function bakerRaw(lambda, phi) {
var phi0 = abs(phi);
return phi0 < quarterPi
? [lambda, log(tan(quarterPi + phi / 2))]
: [lambda * cos(phi0) * (2 * sqrt2 - 1 / sin(phi0)), sign(phi) * (2 * sqrt2 * (phi0 - quarterPi) - log(tan(phi0 / 2)))];
}
bakerRaw.invert = function(x, y) {
if ((y0 = abs(y)) < phi0) return [x, 2 * atan(exp(y)) - halfPi];
var phi = quarterPi, i = 25, delta, y0;
do {
var cosPhi_2 = cos(phi / 2), tanPhi_2 = tan(phi / 2);
phi -= delta = (sqrt8 * (phi - quarterPi) - log(tanPhi_2) - y0) / (sqrt8 - cosPhi_2 * cosPhi_2 / (2 * tanPhi_2));
} while (abs(delta) > epsilon2 && --i > 0);
return [x / (cos(phi) * (sqrt8 - 1 / sin(phi))), sign(y) * phi];
};
var baker = function() {
return d3Geo.geoProjection(bakerRaw)
.scale(112.314);
};
function berghausRaw(lobes) {
var k = 2 * pi / lobes;
function forward(lambda, phi) {
var p = d3Geo.geoAzimuthalEquidistantRaw(lambda, phi);
if (abs(lambda) > halfPi) { // back hemisphere
var theta = atan2(p[1], p[0]),
r = sqrt(p[0] * p[0] + p[1] * p[1]),
theta0 = k * round((theta - halfPi) / k) + halfPi,
α = atan2(sin(theta -= theta0), 2 - cos(theta)); // angle relative to lobe end
theta = theta0 + asin(pi / r * sin(α)) - α;
p[0] = r * cos(theta);
p[1] = r * sin(theta);
}
return p;
}
forward.invert = function(x, y) {
var r = sqrt(x * x + y * y);
if (r > halfPi) {
var theta = atan2(y, x),
theta0 = k * round((theta - halfPi) / k) + halfPi,
s = theta > theta0 ? -1 : 1,
A = r * cos(theta0 - theta),
cotα = 1 / tan(s * acos((A - pi) / sqrt(pi * (pi - 2 * A) + r * r)));
theta = theta0 + 2 * atan((cotα + s * sqrt(cotα * cotα - 3)) / 3);
x = r * cos(theta), y = r * sin(theta);
}
return d3Geo.geoAzimuthalEquidistantRaw.invert(x, y);
};
return forward;
}
var berghaus = function() {
var lobes = 5,
m = d3Geo.geoProjectionMutator(berghausRaw),
p = m(lobes),
projectionStream = p.stream,
epsilon$$1 = 1e-2,
cr = -cos(epsilon$$1 * radians),
sr = sin(epsilon$$1 * radians);
p.lobes = function(_) {
return arguments.length ? m(lobes = +_) : lobes;
};
p.stream = function(stream) {
var rotate = p.rotate(),
rotateStream = projectionStream(stream),
sphereStream = (p.rotate([0, 0]), projectionStream(stream));
p.rotate(rotate);
rotateStream.sphere = function() {
sphereStream.polygonStart(), sphereStream.lineStart();
for (var i = 0, delta = 360 / lobes, delta0 = 2 * pi / lobes, phi = 90 - 180 / lobes, phi0 = halfPi; i < lobes; ++i, phi -= delta, phi0 -= delta0) {
sphereStream.point(atan2(sr * cos(phi0), cr) * degrees, asin(sr * sin(phi0)) * degrees);
if (phi < -90) {
sphereStream.point(-90, -180 - phi - epsilon$$1);
sphereStream.point(-90, -180 - phi + epsilon$$1);
} else {
sphereStream.point(90, phi + epsilon$$1);
sphereStream.point(90, phi - epsilon$$1);
}
}
sphereStream.lineEnd(), sphereStream.polygonEnd();
};
return rotateStream;
};
return p
.scale(87.8076)
.center([0, 17.1875])
.clipAngle(180 - 1e-3);
};
function mollweideBromleyTheta(cp, phi) {
var cpsinPhi = cp * sin(phi), i = 30, delta;
do phi -= delta = (phi + sin(phi) - cpsinPhi) / (1 + cos(phi));
while (abs(delta) > epsilon && --i > 0);
return phi / 2;
}
function mollweideBromleyRaw(cx, cy, cp) {
function forward(lambda, phi) {
return [cx * lambda * cos(phi = mollweideBromleyTheta(cp, phi)), cy * sin(phi)];
}
forward.invert = function(x, y) {
return y = asin(y / cy), [x / (cx * cos(y)), asin((2 * y + sin(2 * y)) / cp)];
};
return forward;
}
var mollweideRaw = mollweideBromleyRaw(sqrt2 / halfPi, sqrt2, pi);
var mollweide = function() {
return d3Geo.geoProjection(mollweideRaw)
.scale(169.529);
};
var k = 2.00276;
var w = 1.11072;
function boggsRaw(lambda, phi) {
var theta = mollweideBromleyTheta(pi, phi);
return [k * lambda / (1 / cos(phi) + w / cos(theta)), (phi + sqrt2 * sin(theta)) / k];
}
boggsRaw.invert = function(x, y) {
var ky = k * y, theta = y < 0 ? -quarterPi : quarterPi, i = 25, delta, phi;
do {
phi = ky - sqrt2 * sin(theta);
theta -= delta = (sin(2 * theta) + 2 * theta - pi * sin(phi)) / (2 * cos(2 * theta) + 2 + pi * cos(phi) * sqrt2 * cos(theta));
} while (abs(delta) > epsilon && --i > 0);
phi = ky - sqrt2 * sin(theta);
return [x * (1 / cos(phi) + w / cos(theta)) / k, phi];
};
var boggs = function() {
return d3Geo.geoProjection(boggsRaw)
.scale(160.857);
};
var parallel1 = function(projectAt) {
var phi0 = 0,
m = d3Geo.geoProjectionMutator(projectAt),
p = m(phi0);
p.parallel = function(_) {
return arguments.length ? m(phi0 = _ * radians) : phi0 * degrees;
};
return p;
};
function sinusoidalRaw(lambda, phi) {
return [lambda * cos(phi), phi];
}
sinusoidalRaw.invert = function(x, y) {
return [x / cos(y), y];
};
var sinusoidal = function() {
return d3Geo.geoProjection(sinusoidalRaw)
.scale(152.63);
};
function bonneRaw(phi0) {
if (!phi0) return sinusoidalRaw;
var cotPhi0 = 1 / tan(phi0);
function forward(lambda, phi) {
var rho = cotPhi0 + phi0 - phi,
e = rho ? lambda * cos(phi) / rho : rho;
return [rho * sin(e), cotPhi0 - rho * cos(e)];
}
forward.invert = function(x, y) {
var rho = sqrt(x * x + (y = cotPhi0 - y) * y),
phi = cotPhi0 + phi0 - rho;
return [rho / cos(phi) * atan2(x, y), phi];
};
return forward;
}
var bonne = function() {
return parallel1(bonneRaw)
.scale(123.082)
.center([0, 26.1441])
.parallel(45);
};
function bottomleyRaw(sinPsi) {
function forward(lambda, phi) {
var rho = halfPi - phi,
eta = rho ? lambda * sinPsi * sin(rho) / rho : rho;
return [rho * sin(eta) / sinPsi, halfPi - rho * cos(eta)];
}
forward.invert = function(x, y) {
var x1 = x * sinPsi,
y1 = halfPi - y,
rho = sqrt(x1 * x1 + y1 * y1),
eta = atan2(x1, y1);
return [(rho ? rho / sin(rho) : 1) * eta / sinPsi, halfPi - rho];
};
return forward;
}
var bottomley = function() {
var sinPsi = 0.5,
m = d3Geo.geoProjectionMutator(bottomleyRaw),
p = m(sinPsi);
p.fraction = function(_) {
return arguments.length ? m(sinPsi = +_) : sinPsi;
};
return p
.scale(158.837);
};
var bromleyRaw = mollweideBromleyRaw(1, 4 / pi, pi);
var bromley = function() {
return d3Geo.geoProjection(bromleyRaw)
.scale(152.63);
};
// Azimuthal distance.
function distance(dPhi, c1, s1, c2, s2, dLambda) {
var cosdLambda = cos(dLambda), r;
if (abs(dPhi) > 1 || abs(dLambda) > 1) {
r = acos(s1 * s2 + c1 * c2 * cosdLambda);
} else {
var sindPhi = sin(dPhi / 2), sindLambda = sin(dLambda / 2);
r = 2 * asin(sqrt(sindPhi * sindPhi + c1 * c2 * sindLambda * sindLambda));
}
return abs(r) > epsilon ? [r, atan2(c2 * sin(dLambda), c1 * s2 - s1 * c2 * cosdLambda)] : [0, 0];
}
// Angle opposite a, and contained between sides of lengths b and c.
function angle(b, c, a) {
return acos((b * b + c * c - a * a) / (2 * b * c));
}
// Normalize longitude.
function longitude(lambda) {
return lambda - 2 * pi * floor((lambda + pi) / (2 * pi));
}
function chamberlinRaw(p0, p1, p2) {
var points = [
[p0[0], p0[1], sin(p0[1]), cos(p0[1])],
[p1[0], p1[1], sin(p1[1]), cos(p1[1])],
[p2[0], p2[1], sin(p2[1]), cos(p2[1])]
];
for (var a = points[2], b, i = 0; i < 3; ++i, a = b) {
b = points[i];
a.v = distance(b[1] - a[1], a[3], a[2], b[3], b[2], b[0] - a[0]);
a.point = [0, 0];
}
var beta0 = angle(points[0].v[0], points[2].v[0], points[1].v[0]),
beta1 = angle(points[0].v[0], points[1].v[0], points[2].v[0]),
beta2 = pi - beta0;
points[2].point[1] = 0;
points[0].point[0] = -(points[1].point[0] = points[0].v[0] / 2);
var mean = [
points[2].point[0] = points[0].point[0] + points[2].v[0] * cos(beta0),
2 * (points[0].point[1] = points[1].point[1] = points[2].v[0] * sin(beta0))
];
function forward(lambda, phi) {
var sinPhi = sin(phi),
cosPhi = cos(phi),
v = new Array(3), i;
// Compute distance and azimuth from control points.
for (i = 0; i < 3; ++i) {
var p = points[i];
v[i] = distance(phi - p[1], p[3], p[2], cosPhi, sinPhi, lambda - p[0]);
if (!v[i][0]) return p.point;
v[i][1] = longitude(v[i][1] - p.v[1]);
}
// Arithmetic mean of interception points.
var point = mean.slice();
for (i = 0; i < 3; ++i) {
var j = i == 2 ? 0 : i + 1;
var a = angle(points[i].v[0], v[i][0], v[j][0]);
if (v[i][1] < 0) a = -a;
if (!i) {
point[0] += v[i][0] * cos(a);
point[1] -= v[i][0] * sin(a);
} else if (i == 1) {
a = beta1 - a;
point[0] -= v[i][0] * cos(a);
point[1] -= v[i][0] * sin(a);
} else {
a = beta2 - a;
point[0] += v[i][0] * cos(a);
point[1] += v[i][0] * sin(a);
}
}
point[0] /= 3, point[1] /= 3;
return point;
}
return forward;
}
function pointRadians(p) {
return p[0] *= radians, p[1] *= radians, p;
}
function chamberlinAfrica() {
return chamberlin([0, 22], [45, 22], [22.5, -22])
.scale(380)
.center([22.5, 2]);
}
function chamberlin(p0, p1, p2) { // TODO order matters!
var c = d3Geo.geoCentroid({type: "MultiPoint", coordinates: [p0, p1, p2]}),
R = [-c[0], -c[1]],
r = d3Geo.geoRotation(R),
p = d3Geo.geoProjection(chamberlinRaw(pointRadians(r(p0)), pointRadians(r(p1)), pointRadians(r(p2)))).rotate(R),
center = p.center;
delete p.rotate;
p.center = function(_) {
return arguments.length ? center(r(_)) : r.invert(center());
};
return p
.clipAngle(90);
}
function collignonRaw(lambda, phi) {
var alpha = sqrt(1 - sin(phi));
return [(2 / sqrtPi) * lambda * alpha, sqrtPi * (1 - alpha)];
}
collignonRaw.invert = function(x, y) {
var lambda = (lambda = y / sqrtPi - 1) * lambda;
return [lambda > 0 ? x * sqrt(pi / lambda) / 2 : 0, asin(1 - lambda)];
};
var collignon = function() {
return d3Geo.geoProjection(collignonRaw)
.scale(95.6464)
.center([0, 30]);
};
function craigRaw(phi0) {
var tanPhi0 = tan(phi0);
function forward(lambda, phi) {
return [lambda, (lambda ? lambda / sin(lambda) : 1) * (sin(phi) * cos(lambda) - tanPhi0 * cos(phi))];
}
forward.invert = tanPhi0 ? function(x, y) {
if (x) y *= sin(x) / x;
var cosλ = cos(x);
return [x, 2 * atan2(sqrt(cosλ * cosλ + tanPhi0 * tanPhi0 - y * y) - cosλ, tanPhi0 - y)];
} : function(x, y) {
return [x, asin(x ? y * tan(x) / x : y)];
};
return forward;
}
var craig = function() {
return parallel1(craigRaw)
.scale(249.828)
.clipAngle(90);
};
var sqrt3 = sqrt(3);
function crasterRaw(lambda, phi) {
return [sqrt3 * lambda * (2 * cos(2 * phi / 3) - 1) / sqrtPi, sqrt3 * sqrtPi * sin(phi / 3)];
}
crasterRaw.invert = function(x, y) {
var phi = 3 * asin(y / (sqrt3 * sqrtPi));
return [sqrtPi * x / (sqrt3 * (2 * cos(2 * phi / 3) - 1)), phi];
};
var craster = function() {
return d3Geo.geoProjection(crasterRaw)
.scale(156.19);
};
function cylindricalEqualAreaRaw(phi0) {
var cosPhi0 = cos(phi0);
function forward(lambda, phi) {
return [lambda * cosPhi0, sin(phi) / cosPhi0];
}
forward.invert = function(x, y) {
return [x / cosPhi0, asin(y * cosPhi0)];
};
return forward;
}
var cylindricalEqualArea = function() {
return parallel1(cylindricalEqualAreaRaw)
.parallel(38.58) // acos(sqrt(width / height / pi)) * radians
.scale(195.044); // width / (sqrt(width / height / pi) * 2 * pi)
};
function cylindricalStereographicRaw(phi0) {
var cosPhi0 = cos(phi0);
function forward(lambda, phi) {
return [lambda * cosPhi0, (1 + cosPhi0) * tan(phi / 2)];
}
forward.invert = function(x, y) {
return [x / cosPhi0, atan(y / (1 + cosPhi0)) * 2];
};
return forward;
}
var cylindricalStereographic = function() {
return parallel1(cylindricalStereographicRaw)
.scale(124.75);
};
function eckert1Raw(lambda, phi) {
var alpha = sqrt(8 / (3 * pi));
return [
alpha * lambda * (1 - abs(phi) / pi),
alpha * phi
];
}
eckert1Raw.invert = function(x, y) {
var alpha = sqrt(8 / (3 * pi)),
phi = y / alpha;
return [
x / (alpha * (1 - abs(phi) / pi)),
phi
];
};
var eckert1 = function() {
return d3Geo.geoProjection(eckert1Raw)
.scale(165.664);
};
function eckert2Raw(lambda, phi) {
var alpha = sqrt(4 - 3 * sin(abs(phi)));
return [
2 / sqrt(6 * pi) * lambda * alpha,
sign(phi) * sqrt(2 * pi / 3) * (2 - alpha)
];
}
eckert2Raw.invert = function(x, y) {
var alpha = 2 - abs(y) / sqrt(2 * pi / 3);
return [
x * sqrt(6 * pi) / (2 * alpha),
sign(y) * asin((4 - alpha * alpha) / 3)
];
};
var eckert2 = function() {
return d3Geo.geoProjection(eckert2Raw)
.scale(165.664);
};
function eckert3Raw(lambda, phi) {
var k = sqrt(pi * (4 + pi));
return [
2 / k * lambda * (1 + sqrt(1 - 4 * phi * phi / (pi * pi))),
4 / k * phi
];
}
eckert3Raw.invert = function(x, y) {
var k = sqrt(pi * (4 + pi)) / 2;
return [
x * k / (1 + sqrt(1 - y * y * (4 + pi) / (4 * pi))),
y * k / 2
];
};
var eckert3 = function() {
return d3Geo.geoProjection(eckert3Raw)
.scale(180.739);
};
function eckert4Raw(lambda, phi) {
var k = (2 + halfPi) * sin(phi);
phi /= 2;
for (var i = 0, delta = Infinity; i < 10 && abs(delta) > epsilon; i++) {
var cosPhi = cos(phi);
phi -= delta = (phi + sin(phi) * (cosPhi + 2) - k) / (2 * cosPhi * (1 + cosPhi));
}
return [
2 / sqrt(pi * (4 + pi)) * lambda * (1 + cos(phi)),
2 * sqrt(pi / (4 + pi)) * sin(phi)
];
}
eckert4Raw.invert = function(x, y) {
var A = y * sqrt((4 + pi) / pi) / 2,
k = asin(A),
c = cos(k);
return [
x / (2 / sqrt(pi * (4 + pi)) * (1 + c)),
asin((k + A * (c + 2)) / (2 + halfPi))
];
};
var eckert4 = function() {
return d3Geo.geoProjection(eckert4Raw)
.scale(180.739);
};
function eckert5Raw(lambda, phi) {
return [
lambda * (1 + cos(phi)) / sqrt(2 + pi),
2 * phi / sqrt(2 + pi)
];
}
eckert5Raw.invert = function(x, y) {
var k = sqrt(2 + pi),
phi = y * k / 2;
return [
k * x / (1 + cos(phi)),
phi
];
};
var eckert5 = function() {
return d3Geo.geoProjection(eckert5Raw)
.scale(173.044);
};
function eckert6Raw(lambda, phi) {
var k = (1 + halfPi) * sin(phi);
for (var i = 0, delta = Infinity; i < 10 && abs(delta) > epsilon; i++) {
phi -= delta = (phi + sin(phi) - k) / (1 + cos(phi));
}
k = sqrt(2 + pi);
return [
lambda * (1 + cos(phi)) / k,
2 * phi / k
];
}
eckert6Raw.invert = function(x, y) {
var j = 1 + halfPi,
k = sqrt(j / 2);
return [
x * 2 * k / (1 + cos(y *= k)),
asin((y + sin(y)) / j)
];
};
var eckert6 = function() {
return d3Geo.geoProjection(eckert6Raw)
.scale(173.044);
};
var eisenlohrK = 3 + 2 * sqrt2;
function eisenlohrRaw(lambda, phi) {
var s0 = sin(lambda /= 2),
c0 = cos(lambda),
k = sqrt(cos(phi)),
c1 = cos(phi /= 2),
t = sin(phi) / (c1 + sqrt2 * c0 * k),
c = sqrt(2 / (1 + t * t)),
v = sqrt((sqrt2 * c1 + (c0 + s0) * k) / (sqrt2 * c1 + (c0 - s0) * k));
return [
eisenlohrK * (c * (v - 1 / v) - 2 * log(v)),
eisenlohrK * (c * t * (v + 1 / v) - 2 * atan(t))
];
}
eisenlohrRaw.invert = function(x, y) {
if (!(p = augustRaw.invert(x / 1.2, y * 1.065))) return null;
var lambda = p[0], phi = p[1], i = 20, p;
x /= eisenlohrK, y /= eisenlohrK;
do {
var _0 = lambda / 2,
_1 = phi / 2,
s0 = sin(_0),
c0 = cos(_0),
s1 = sin(_1),
c1 = cos(_1),
cos1 = cos(phi),
k = sqrt(cos1),
t = s1 / (c1 + sqrt2 * c0 * k),
t2 = t * t,
c = sqrt(2 / (1 + t2)),
v0 = (sqrt2 * c1 + (c0 + s0) * k),
v1 = (sqrt2 * c1 + (c0 - s0) * k),
v2 = v0 / v1,
v = sqrt(v2),
vm1v = v - 1 / v,
vp1v = v + 1 / v,
fx = c * vm1v - 2 * log(v) - x,
fy = c * t * vp1v - 2 * atan(t) - y,
deltatDeltaLambda = s1 && sqrt1_2 * k * s0 * t2 / s1,
deltatDeltaPhi = (sqrt2 * c0 * c1 + k) / (2 * (c1 + sqrt2 * c0 * k) * (c1 + sqrt2 * c0 * k) * k),
deltacDeltat = -0.5 * t * c * c * c,
deltacDeltaLambda = deltacDeltat * deltatDeltaLambda,
deltacDeltaPhi = deltacDeltat * deltatDeltaPhi,
A = (A = 2 * c1 + sqrt2 * k * (c0 - s0)) * A * v,
deltavDeltaLambda = (sqrt2 * c0 * c1 * k + cos1) / A,
deltavDeltaPhi = -(sqrt2 * s0 * s1) / (k * A),
deltaxDeltaLambda = vm1v * deltacDeltaLambda - 2 * deltavDeltaLambda / v + c * (deltavDeltaLambda + deltavDeltaLambda / v2),
deltaxDeltaPhi = vm1v * deltacDeltaPhi - 2 * deltavDeltaPhi / v + c * (deltavDeltaPhi + deltavDeltaPhi / v2),
deltayDeltaLambda = t * vp1v * deltacDeltaLambda - 2 * deltatDeltaLambda / (1 + t2) + c * vp1v * deltatDeltaLambda + c * t * (deltavDeltaLambda - deltavDeltaLambda / v2),
deltayDeltaPhi = t * vp1v * deltacDeltaPhi - 2 * deltatDeltaPhi / (1 + t2) + c * vp1v * deltatDeltaPhi + c * t * (deltavDeltaPhi - deltavDeltaPhi / v2),
denominator = deltaxDeltaPhi * deltayDeltaLambda - deltayDeltaPhi * deltaxDeltaLambda;
if (!denominator) break;
var deltaLambda = (fy * deltaxDeltaPhi - fx * deltayDeltaPhi) / denominator,
deltaPhi = (fx * deltayDeltaLambda - fy * deltaxDeltaLambda) / denominator;
lambda -= deltaLambda;
phi = max(-halfPi, min(halfPi, phi - deltaPhi));
} while ((abs(deltaLambda) > epsilon || abs(deltaPhi) > epsilon) && --i > 0);
return abs(abs(phi) - halfPi) < epsilon ? [0, phi] : i && [lambda, phi];
};
var eisenlohr = function() {
return d3Geo.geoProjection(eisenlohrRaw)
.scale(62.5271);
};
var faheyK = cos(35 * radians);
function faheyRaw(lambda, phi) {
var t = tan(phi / 2);
return [lambda * faheyK * sqrt(1 - t * t), (1 + faheyK) * t];
}
faheyRaw.invert = function(x, y) {
var t = y / (1 + faheyK);
return [x && x / (faheyK * sqrt(1 - t * t)), 2 * atan(t)];
};
var fahey = function() {
return d3Geo.geoProjection(faheyRaw)
.scale(137.152);
};
function foucautRaw(lambda, phi) {
var k = phi / 2, cosk = cos(k);
return [ 2 * lambda / sqrtPi * cos(phi) * cosk * cosk, sqrtPi * tan(k)];
}
foucautRaw.invert = function(x, y) {
var k = atan(y / sqrtPi), cosk = cos(k), phi = 2 * k;
return [x * sqrtPi / 2 / (cos(phi) * cosk * cosk), phi];
};
var foucaut = function() {
return d3Geo.geoProjection(foucautRaw)
.scale(135.264);
};
function gilbertForward(point) {
return [point[0] / 2, asin(tan(point[1] / 2 * radians)) * degrees];
}
function gilbertInvert(point) {
return [point[0] * 2, 2 * atan(sin(point[1] * radians)) * degrees];
}
var gilbert = function(projectionType) {
if (projectionType == null) projectionType = d3Geo.geoOrthographic;
var projection = projectionType(),
equirectangular = d3Geo.geoEquirectangular().scale(degrees).precision(0).clipAngle(null).translate([0, 0]); // antimeridian cutting
function gilbert(point) {
return projection(gilbertForward(point));
}
if (projection.invert) gilbert.invert = function(point) {
return gilbertInvert(projection.invert(point));
};
gilbert.stream = function(stream) {
var s1 = projection.stream(stream), s0 = equirectangular.stream({
point: function(lambda, phi) { s1.point(lambda / 2, asin(tan(-phi / 2 * radians)) * degrees); },
lineStart: function() { s1.lineStart(); },
lineEnd: function() { s1.lineEnd(); },
polygonStart: function() { s1.polygonStart(); },
polygonEnd: function() { s1.polygonEnd(); }
});
s0.sphere = s1.sphere;
return s0;
};
function property(name) {
gilbert[name] = function(_) {
return arguments.length ? (projection[name](_), gilbert) : projection[name]();
};
}
gilbert.rotate = function(_) {
return arguments.length ? (equirectangular.rotate(_), gilbert) : equirectangular.rotate();
};
gilbert.center = function(_) {
return arguments.length ? (projection.center(gilbertForward(_)), gilbert) : gilbertInvert(projection.center());
};
property("clipAngle");
property("clipExtent");
property("scale");
property("translate");
property("precision");
return gilbert
.scale(249.5);
};
function gingeryRaw(rho, n) {
var k = 2 * pi / n,
rho2 = rho * rho;
function forward(lambda, phi) {
var p = d3Geo.geoAzimuthalEquidistantRaw(lambda, phi),
x = p[0],
y = p[1],
r2 = x * x + y * y;
if (r2 > rho2) {
var r = sqrt(r2),
theta = atan2(y, x),
theta0 = k * round(theta / k),
alpha = theta - theta0,
rhoCosAlpha = rho * cos(alpha),
k_ = (rho * sin(alpha) - alpha * sin(rhoCosAlpha)) / (halfPi - rhoCosAlpha),
s_ = gingeryLength(alpha, k_),
e = (pi - rho) / gingeryIntegrate(s_, rhoCosAlpha, pi);
x = r;
var i = 50, delta;
do {
x -= delta = (rho + gingeryIntegrate(s_, rhoCosAlpha, x) * e - r) / (s_(x) * e);
} while (abs(delta) > epsilon && --i > 0);
y = alpha * sin(x);
if (x < halfPi) y -= k_ * (x - halfPi);
var s = sin(theta0),
c = cos(theta0);
p[0] = x * c - y * s;
p[1] = x * s + y * c;
}
return p;
}
forward.invert = function(x, y) {
var r2 = x * x + y * y;
if (r2 > rho2) {
var r = sqrt(r2),
theta = atan2(y, x),
theta0 = k * round(theta / k),
dTheta = theta - theta0;
x = r * cos(dTheta);
y = r * sin(dTheta);
var x_halfPi = x - halfPi,
sinx = sin(x),
alpha = y / sinx,
delta = x < halfPi ? Infinity : 0,
i = 10;
while (true) {
var rhosinAlpha = rho * sin(alpha),
rhoCosAlpha = rho * cos(alpha),
sinRhoCosAlpha = sin(rhoCosAlpha),
halfPi_RhoCosAlpha = halfPi - rhoCosAlpha,
k_ = (rhosinAlpha - alpha * sinRhoCosAlpha) / halfPi_RhoCosAlpha,
s_ = gingeryLength(alpha, k_);
if (abs(delta) < epsilon2 || !--i) break;
alpha -= delta = (alpha * sinx - k_ * x_halfPi - y) / (
sinx - x_halfPi * 2 * (
halfPi_RhoCosAlpha * (rhoCosAlpha + alpha * rhosinAlpha * cos(rhoCosAlpha) - sinRhoCosAlpha) -
rhosinAlpha * (rhosinAlpha - alpha * sinRhoCosAlpha)
) / (halfPi_RhoCosAlpha * halfPi_RhoCosAlpha));
}
r = rho + gingeryIntegrate(s_, rhoCosAlpha, x) * (pi - rho) / gingeryIntegrate(s_, rhoCosAlpha, pi);
theta = theta0 + alpha;
x = r * cos(theta);
y = r * sin(theta);
}
return d3Geo.geoAzimuthalEquidistantRaw.invert(x, y);
};
return forward;
}
function gingeryLength(alpha, k) {
return function(x) {
var y_ = alpha * cos(x);
if (x < halfPi) y_ -= k;
return sqrt(1 + y_ * y_);
};
}
// Numerical integration: trapezoidal rule.
function gingeryIntegrate(f, a, b) {
var n = 50,
h = (b - a) / n,
s = f(a) + f(b);
for (var i = 1, x = a; i < n; ++i) s += 2 * f(x += h);
return s * 0.5 * h;
}
var gingery = function() {
var n = 6,
rho = 30 * radians,
cRho = cos(rho),
sRho = sin(rho),
m = d3Geo.geoProjectionMutator(gingeryRaw),
p = m(rho, n),
stream_ = p.stream,
epsilon$$1 = 1e-2,
cr = -cos(epsilon$$1 * radians),
sr = sin(epsilon$$1 * radians);
p.radius = function(_) {
if (!arguments.length) return rho * degrees;
cRho = cos(rho = _ * radians);
sRho = sin(rho);
return m(rho, n);
};
p.lobes = function(_) {
if (!arguments.length) return n;
return m(rho, n = +_);
};
p.stream = function(stream) {
var rotate = p.rotate(),
rotateStream = stream_(stream),
sphereStream = (p.rotate([0, 0]), stream_(stream));
p.rotate(rotate);
rotateStream.sphere = function() {
sphereStream.polygonStart(), sphereStream.lineStart();
for (var i = 0, delta = 2 * pi / n, phi = 0; i < n; ++i, phi -= delta) {
sphereStream.point(atan2(sr * cos(phi), cr) * degrees, asin(sr * sin(phi)) * degrees);
sphereStream.point(atan2(sRho * cos(phi - delta / 2), cRho) * degrees, asin(sRho * sin(phi - delta / 2)) * degrees);
}
sphereStream.lineEnd(), sphereStream.polygonEnd();
};
return rotateStream;
};
return p
.rotate([90, -40])
.scale(91.7095)
.clipAngle(180 - 1e-3);
};
var ginzburgPolyconicRaw = function(a, b, c, d, e, f, g, h) {
if (arguments.length < 8) h = 0;
function forward(lambda, phi) {
if (!phi) return [a * lambda / pi, 0];
var phi2 = phi * phi,
xB = a + phi2 * (b + phi2 * (c + phi2 * d)),
yB = phi * (e - 1 + phi2 * (f - h + phi2 * g)),
m = (xB * xB + yB * yB) / (2 * yB),
alpha = lambda * asin(xB / m) / pi;
return [m * sin(alpha), phi * (1 + phi2 * h) + m * (1 - cos(alpha))];
}
forward.invert = function(x, y) {
var lambda = pi * x / a,
phi = y,
deltaLambda, deltaPhi, i = 50;
do {
var phi2 = phi * phi,
xB = a + phi2 * (b + phi2 * (c + phi2 * d)),
yB = phi * (e - 1 + phi2 * (f - h + phi2 * g)),
p = xB * xB + yB * yB,
q = 2 * yB,
m = p / q,
m2 = m * m,
dAlphadLambda = asin(xB / m) / pi,
alpha = lambda * dAlphadLambda,
xB2 = xB * xB,
dxBdPhi = (2 * b + phi2 * (4 * c + phi2 * 6 * d)) * phi,
dyBdPhi = e + phi2 * (3 * f + phi2 * 5 * g),
dpdPhi = 2 * (xB * dxBdPhi + yB * (dyBdPhi - 1)),
dqdPhi = 2 * (dyBdPhi - 1),
dmdPhi = (dpdPhi * q - p * dqdPhi) / (q * q),
cosAlpha = cos(alpha),
sinAlpha = sin(alpha),
mcosAlpha = m * cosAlpha,
msinAlpha = m * sinAlpha,
dAlphadPhi = ((lambda / pi) * (1 / sqrt(1 - xB2 / m2)) * (dxBdPhi * m - xB * dmdPhi)) / m2,
fx = msinAlpha - x,
fy = phi * (1 + phi2 * h) + m - mcosAlpha - y,
deltaxDeltaPhi = dmdPhi * sinAlpha + mcosAlpha * dAlphadPhi,
deltaxDeltaLambda = mcosAlpha * dAlphadLambda,
deltayDeltaPhi = 1 + dmdPhi - (dmdPhi * cosAlpha - msinAlpha * dAlphadPhi),
deltayDeltaLambda = msinAlpha * dAlphadLambda,
denominator = deltaxDeltaPhi * deltayDeltaLambda - deltayDeltaPhi * deltaxDeltaLambda;
if (!denominator) break;
lambda -= deltaLambda = (fy * deltaxDeltaPhi - fx * deltayDeltaPhi) / denominator;
phi -= deltaPhi = (fx * deltayDeltaLambda - fy * deltaxDeltaLambda) / denominator;
} while ((abs(deltaLambda) > epsilon || abs(deltaPhi) > epsilon) && --i > 0);
return [lambda, phi];
};
return forward;
};
var ginzburg4Raw = ginzburgPolyconicRaw(2.8284, -1.6988, 0.75432, -0.18071, 1.76003, -0.38914, 0.042555);
var ginzburg4 = function() {
return d3Geo.geoProjection(ginzburg4Raw)
.scale(149.995);
};
var ginzburg5Raw = ginzburgPolyconicRaw(2.583819, -0.835827, 0.170354, -0.038094, 1.543313, -0.411435,0.082742);
var ginzburg5 = function() {
return d3Geo.geoProjection(ginzburg5Raw)
.scale(153.93);
};
var ginzburg6Raw = ginzburgPolyconicRaw(5 / 6 * pi, -0.62636, -0.0344, 0, 1.3493, -0.05524, 0, 0.045);
var ginzburg6 = function() {
return d3Geo.geoProjection(ginzburg6Raw)
.scale(130.945);
};
function ginzburg8Raw(lambda, phi) {
var lambda2 = lambda * lambda,
phi2 = phi * phi;
return [
lambda * (1 - 0.162388 * phi2) * (0.87 - 0.000952426 * lambda2 * lambda2),
phi * (1 + phi2 / 12)
];
}
ginzburg8Raw.invert = function(x, y) {
var lambda = x,
phi = y,
i = 50, delta;
do {
var phi2 = phi * phi;
phi -= delta = (phi * (1 + phi2 / 12) - y) / (1 + phi2 / 4);
} while (abs(delta) > epsilon && --i > 0);
i = 50;
x /= 1 -0.162388 * phi2;
do {
var lambda4 = (lambda4 = lambda * lambda) * lambda4;
lambda -= delta = (lambda * (0.87 - 0.000952426 * lambda4) - x) / (0.87 - 0.00476213 * lambda4);
} while (abs(delta) > epsilon && --i > 0);
return [lambda, phi];
};
var ginzburg8 = function() {
return d3Geo.geoProjection(ginzburg8Raw)
.scale(131.747);
};
var ginzburg9Raw = ginzburgPolyconicRaw(2.6516, -0.76534, 0.19123, -0.047094, 1.36289, -0.13965,0.031762);
var ginzburg9 = function() {
return d3Geo.geoProjection(ginzburg9Raw)
.scale(131.087);
};
var squareRaw = function(project) {
var dx = project(halfPi, 0)[0] - project(-halfPi, 0)[0];
function projectSquare(lambda, phi) {
var s = lambda > 0 ? -0.5 : 0.5,
point = project(lambda + s * pi, phi);
point[0] -= s * dx;
return point;
}
if (project.invert) projectSquare.invert = function(x, y) {
var s = x > 0 ? -0.5 : 0.5,
location = project.invert(x + s * dx, y),
lambda = location[0] - s * pi;
if (lambda < -pi) lambda += 2 * pi;
else if (lambda > pi) lambda -= 2 * pi;
location[0] = lambda;
return location;
};
return projectSquare;
};
function gringortenRaw(lambda, phi) {
var sLambda = sign(lambda),
sPhi = sign(phi),
cosPhi = cos(phi),
x = cos(lambda) * cosPhi,
y = sin(lambda) * cosPhi,
z = sin(sPhi * phi);
lambda = abs(atan2(y, z));
phi = asin(x);
if (abs(lambda - halfPi) > epsilon) lambda %= halfPi;
var point = gringortenHexadecant(lambda > pi / 4 ? halfPi - lambda : lambda, phi);
if (lambda > pi / 4) z = point[0], point[0] = -point[1], point[1] = -z;
return (point[0] *= sLambda, point[1] *= -sPhi, point);
}
gringortenRaw.invert = function(x, y) {
var sx = sign(x),
sy = sign(y),
x0 = -sx * x,
y0 = -sy * y,
t = y0 / x0 < 1,
p = gringortenHexadecantInvert(t ? y0 : x0, t ? x0 : y0),
lambda = p[0],
phi = p[1],
cosPhi = cos(phi);
if (t) lambda = -halfPi - lambda;
return [sx * (atan2(sin(lambda) * cosPhi, -sin(phi)) + pi), sy * asin(cos(lambda) * cosPhi)];
};
function gringortenHexadecant(lambda, phi) {
if (phi === halfPi) return [0, 0];
var sinPhi = sin(phi),
r = sinPhi * sinPhi,
r2 = r * r,
j = 1 + r2,
k = 1 + 3 * r2,
q = 1 - r2,
z = asin(1 / sqrt(j)),
v = q + r * j * z,
p2 = (1 - sinPhi) / v,
p = sqrt(p2),
a2 = p2 * j,
a = sqrt(a2),
h = p * q,
x,
i;
if (lambda === 0) return [0, -(h + r * a)];
var cosPhi = cos(phi),
secPhi = 1 / cosPhi,
drdPhi = 2 * sinPhi * cosPhi,
dvdPhi = (-3 * r + z * k) * drdPhi,
dp2dPhi = (-v * cosPhi - (1 - sinPhi) * dvdPhi) / (v * v),
dpdPhi = (0.5 * dp2dPhi) / p,
dhdPhi = q * dpdPhi - 2 * r * p * drdPhi,
dra2dPhi = r * j * dp2dPhi + p2 * k * drdPhi,
mu = -secPhi * drdPhi,
nu = -secPhi * dra2dPhi,
zeta = -2 * secPhi * dhdPhi,
lambda1 = 4 * lambda / pi,
delta;
// Slower but accurate bisection method.
if (lambda > 0.222 * pi || phi < pi / 4 && lambda > 0.175 * pi) {
x = (h + r * sqrt(a2 * (1 + r2) - h * h)) / (1 + r2);
if (lambda > pi / 4) return [x, x];
var x1 = x, x0 = 0.5 * x;
x = 0.5 * (x0 + x1), i = 50;
do {
var g = sqrt(a2 - x * x),
f = (x * (zeta + mu * g) + nu * asin(x / a)) - lambda1;
if (!f) break;
if (f < 0) x0 = x;
else x1 = x;
x = 0.5 * (x0 + x1);
} while (abs(x1 - x0) > epsilon && --i > 0);
}
// Newton-Raphson.
else {
x = epsilon, i = 25;
do {
var x2 = x * x,
g2 = sqrt(a2 - x2),
zetaMug = zeta + mu * g2,
f2 = x * zetaMug + nu * asin(x / a) - lambda1,
df = zetaMug + (nu - mu * x2) / g2;
x -= delta = g2 ? f2 / df : 0;
} while (abs(delta) > epsilon && --i > 0);
}
return [x, -h - r * sqrt(a2 - x * x)];
}
function gringortenHexadecantInvert(x, y) {
var x0 = 0,
x1 = 1,
r = 0.5,
i = 50;
while (true) {
var r2 = r * r,
sinPhi = sqrt(r),
z = asin(1 / sqrt(1 + r2)),
v = (1 - r2) + r * (1 + r2) * z,
p2 = (1 - sinPhi) / v,
p = sqrt(p2),
a2 = p2 * (1 + r2),
h = p * (1 - r2),
g2 = a2 - x * x,
g = sqrt(g2),
y0 = y + h + r * g;
if (abs(x1 - x0) < epsilon2 || --i === 0 || y0 === 0) break;
if (y0 > 0) x0 = r;
else x1 = r;
r = 0.5 * (x0 + x1);
}
if (!i) return null;
var phi = asin(sinPhi),
cosPhi = cos(phi),
secPhi = 1 / cosPhi,
drdPhi = 2 * sinPhi * cosPhi,
dvdPhi = (-3 * r + z * (1 + 3 * r2)) * drdPhi,
dp2dPhi = (-v * cosPhi - (1 - sinPhi) * dvdPhi) / (v * v),
dpdPhi = 0.5 * dp2dPhi / p,
dhdPhi = (1 - r2) * dpdPhi - 2 * r * p * drdPhi,
zeta = -2 * secPhi * dhdPhi,
mu = -secPhi * drdPhi,
nu = -secPhi * (r * (1 + r2) * dp2dPhi + p2 * (1 + 3 * r2) * drdPhi);
return [pi / 4 * (x * (zeta + mu * g) + nu * asin(x / sqrt(a2))), phi];
}
var gringorten = function() {
return d3Geo.geoProjection(squareRaw(gringortenRaw))
.scale(239.75);
};
// Returns [sn, cn, dn](u + iv|m).
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).
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);
return [sin(phi), t = cos(phi), t / cos(phi - b), phi];
}
// Calculate F(phi+iPsi|m).
// See Abramowitz and Stegun, 17.4.11.
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² = sin²α.
// See Abramowitz and Stegun, 17.6.7.
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);
}
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
];
};
var guyou = function() {
return d3Geo.geoProjection(squareRaw(guyouRaw))
.scale(151.496);
};
function hammerRaw(A, B) {
if (arguments.length < 2) B = A;
if (B === 1) return d3Geo.geoAzimuthalEqualAreaRaw;
if (B === Infinity) return hammerQuarticAuthalicRaw;
function forward(lambda, phi) {
var coordinates = d3Geo.geoAzimuthalEqualAreaRaw(lambda / B, phi);
coordinates[0] *= A;
return coordinates;
}
forward.invert = function(x, y) {
var coordinates = d3Geo.geoAzimuthalEqualAreaRaw.invert(x / A, y);
coordinates[0] *= B;
return coordinates;
};
return forward;
}
function hammerQuarticAuthalicRaw(lambda, phi) {
return [
lambda * cos(phi) / cos(phi /= 2),
2 * sin(phi)
];
}
hammerQuarticAuthalicRaw.invert = function(x, y) {
var phi = 2 * asin(y / 2);
return [
x * cos(phi / 2) / cos(phi),
phi
];
};
var hammer = function() {
var B = 2,
m = d3Geo.geoProjectionMutator(hammerRaw),
p = m(B);
p.coefficient = function(_) {
if (!arguments.length) return B;
return m(B = +_);
};
return p
.scale(169.529);
};
function hammerRetroazimuthalRaw(phi0) {
var sinPhi0 = sin(phi0),
cosPhi0 = cos(phi0),
rotate = hammerRetroazimuthalRotation(phi0);
rotate.invert = hammerRetroazimuthalRotation(-phi0);
function forward(lambda, phi) {
var p = rotate(lambda, phi);
lambda = p[0], phi = p[1];
var sinPhi = sin(phi),
cosPhi = cos(phi),
cosLambda = cos(lambda),
z = acos(sinPhi0 * sinPhi + cosPhi0 * cosPhi * cosLambda),
sinz = sin(z),
K = abs(sinz) > epsilon ? z / sinz : 1;
return [
K * cosPhi0 * sin(lambda),
(abs(lambda) > halfPi ? K : -K) // rotate for back hemisphere
* (sinPhi0 * cosPhi - cosPhi0 * sinPhi * cosLambda)
];
}
forward.invert = function(x, y) {
var rho = sqrt(x * x + y * y),
sinz = -sin(rho),
cosz = cos(rho),
a = rho * cosz,
b = -y * sinz,
c = rho * sinPhi0,
d = sqrt(a * a + b * b - c * c),
phi = atan2(a * c + b * d, b * c - a * d),
lambda = (rho > halfPi ? -1 : 1) * atan2(x * sinz, rho * cos(phi) * cosz + y * sin(phi) * sinz);
return rotate.invert(lambda, phi);
};
return forward;
}
// Latitudinal rotation by phi0.
// Temporary hack until D3 supports arbitrary small-circle clipping origins.
function hammerRetroazimuthalRotation(phi0) {
var sinPhi0 = sin(phi0),
cosPhi0 = cos(phi0);
return function(lambda, phi) {
var cosPhi = cos(phi),
x = cos(lambda) * cosPhi,
y = sin(lambda) * cosPhi,
z = sin(phi);
return [
atan2(y, x * cosPhi0 - z * sinPhi0),
asin(z * cosPhi0 + x * sinPhi0)
];
};
}
var hammerRetroazimuthal = function() {
var phi0 = 0,
m = d3Geo.geoProjectionMutator(hammerRetroazimuthalRaw),
p = m(phi0),
rotate_ = p.rotate,
stream_ = p.stream,
circle = d3Geo.geoCircle();
p.parallel = function(_) {
if (!arguments.length) return phi0 * degrees;
var r = p.rotate();
return m(phi0 = _ * radians).rotate(r);
};
// Temporary hack; see hammerRetroazimuthalRotation.
p.rotate = function(_) {
if (!arguments.length) return (_ = rotate_.call(p), _[1] += phi0 * degrees, _);
rotate_.call(p, [_[0], _[1] - phi0 * degrees]);
circle.center([-_[0], -_[1]]);
return p;
};
p.stream = function(stream) {
stream = stream_(stream);
stream.sphere = function() {
stream.polygonStart();
var epsilon$$1 = 1e-2,
ring = circle.radius(90 - epsilon$$1)().coordinates[0],
n = ring.length - 1,
i = -1,
p;
stream.lineStart();
while (++i < n) stream.point((p = ring[i])[0], p[1]);
stream.lineEnd();
ring = circle.radius(90 + epsilon$$1)().coordinates[0];
n = ring.length - 1;
stream.lineStart();
while (--i >= 0) stream.point((p = ring[i])[0], p[1]);
stream.lineEnd();
stream.polygonEnd();
};
return stream;
};
return p
.scale(79.4187)
.parallel(45)
.clipAngle(180 - 1e-3);
};
var healpixParallel = 41 + 48 / 36 + 37 / 3600;
var healpixLambert = cylindricalEqualAreaRaw(0);
function healpixRaw(H) {
var phi0 = healpixParallel * radians,
dx = collignonRaw(pi, phi0)[0] - collignonRaw(-pi, phi0)[0],
y0 = healpixLambert(0, phi0)[1],
y1 = collignonRaw(0, phi0)[1],
dy1 = sqrtPi - y1,
k = tau / H,
w = 4 / tau,
h = y0 + (dy1 * dy1 * 4) / tau;
function forward(lambda, phi) {
var point,
phi2 = abs(phi);
if (phi2 > phi0) {
var i = min(H - 1, max(0, floor((lambda + pi) / k)));
lambda += pi * (H - 1) / H - i * k;
point = collignonRaw(lambda, phi2);
point[0] = point[0] * tau / dx - tau * (H - 1) / (2 * H) + i * tau / H;
point[1] = y0 + (point[1] - y1) * 4 * dy1 / tau;
if (phi < 0) point[1] = -point[1];
} else {
point = healpixLambert(lambda, phi);
}
point[0] *= w, point[1] /= h;
return point;
}
forward.invert = function(x, y) {
x /= w, y *= h;
var y2 = abs(y);
if (y2 > y0) {
var i = min(H - 1, max(0, floor((x + pi) / k)));
x = (x + pi * (H - 1) / H - i * k) * dx / tau;
var point = collignonRaw.invert(x, 0.25 * (y2 - y0) * tau / dy1 + y1);
point[0] -= pi * (H - 1) / H - i * k;
if (y < 0) point[1] = -point[1];
return point;
}
return healpixLambert.invert(x, y);
};
return forward;
}
function sphere(step) {
return {
type: "Polygon",
coordinates: [
d3Array.range(-180, 180 + step / 2, step).map(function(x, i) { return [x, i & 1 ? 90 - 1e-6 : healpixParallel]; })
.concat(d3Array.range(180, -180 - step / 2, -step).map(function(x, i) { return [x, i & 1 ? -90 + 1e-6 : -healpixParallel]; }))
]
};
}
var healpix = function() {
var H = 4,
m = d3Geo.geoProjectionMutator(healpixRaw),
p = m(H),
stream_ = p.stream;
p.lobes = function(_) {
return arguments.length ? m(H = +_) : H;
};
p.stream = function(stream) {
var rotate = p.rotate(),
rotateStream = stream_(stream),
sphereStream = (p.rotate([0, 0]), stream_(stream));
p.rotate(rotate);
rotateStream.sphere = function() { d3Geo.geoStream(sphere(180 / H), sphereStream); };
return rotateStream;
};
return p
.scale(239.75);
};
function hillRaw(K) {
var L = 1 + K,
sinBt = sin(1 / L),
Bt = asin(sinBt),
A = 2 * sqrt(pi / (B = pi + 4 * Bt * L)),
B,
rho0 = 0.5 * A * (L + sqrt(K * (2 + K))),
K2 = K * K,
L2 = L * L;
function forward(lambda, phi) {
var t = 1 - sin(phi),
rho,
omega;
if (t && t < 2) {
var theta = halfPi - phi, i = 25, delta;
do {
var sinTheta = sin(theta),
cosTheta = cos(theta),
Bt_Bt1 = Bt + atan2(sinTheta, L - cosTheta),
C = 1 + L2 - 2 * L * cosTheta;
theta -= delta = (theta - K2 * Bt - L * sinTheta + C * Bt_Bt1 -0.5 * t * B) / (2 * L * sinTheta * Bt_Bt1);
} while (abs(delta) > epsilon2 && --i > 0);
rho = A * sqrt(C);
omega = lambda * Bt_Bt1 / pi;
} else {
rho = A * (K + t);
omega = lambda * Bt / pi;
}
return [
rho * sin(omega),
rho0 - rho * cos(omega)
];
}
forward.invert = function(x, y) {
var rho2 = x * x + (y -= rho0) * y,
cosTheta = (1 + L2 - rho2 / (A * A)) / (2 * L),
theta = acos(cosTheta),
sinTheta = sin(theta),
Bt_Bt1 = Bt + atan2(sinTheta, L - cosTheta);
return [
asin(x / sqrt(rho2)) * pi / Bt_Bt1,
asin(1 - 2 * (theta - K2 * Bt - L * sinTheta + (1 + L2 - 2 * L * cosTheta) * Bt_Bt1) / B)
];
};
return forward;
}
var hill = function() {
var K = 1,
m = d3Geo.geoProjectionMutator(hillRaw),
p = m(K);
p.ratio = function(_) {
return arguments.length ? m(K = +_) : K;
};
return p
.scale(167.774)
.center([0, 18.67]);
};
var sinuMollweidePhi = 0.7109889596207567;
var sinuMollweideY = 0.0528035274542;
function sinuMollweideRaw(lambda, phi) {
return phi > -sinuMollweidePhi
? (lambda = mollweideRaw(lambda, phi), lambda[1] += sinuMollweideY, lambda)
: sinusoidalRaw(lambda, phi);
}
sinuMollweideRaw.invert = function(x, y) {
return y > -sinuMollweidePhi
? mollweideRaw.invert(x, y - sinuMollweideY)