Skip to content

Instantly share code, notes, and snippets.

@Fil

Fil/.block Secret

Last active August 31, 2017 06:50
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Fil/f511f9e493994fd29c70cf4e24178eaf to your computer and use it in GitHub Desktop.
Save Fil/f511f9e493994fd29c70cf4e24178eaf to your computer and use it in GitHub Desktop.
Download ZIP
flashy voronoi projection badge for the unconf 2017
Raw
.block
license: mit
Display the source blob
Display the rendered blob
Raw
countries.geojson
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Raw
d3-geo-projection.js
// https://d3js.org/d3-geo-projection/ Version 2.2.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,
alpha = atan2(sin(theta -= theta0), 2 - cos(theta)); // angle relative to lobe end
theta = theta0 + asin(pi / r * sin(alpha)) - alpha;
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),
cotAlpha = 1 / tan(s * acos((A - pi) / sqrt(pi * (pi - 2 * A) + r * r)));
theta = theta0 + 2 * atan((cotAlpha + s * sqrt(cotAlpha * cotAlpha - 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 cosLambda = cos(x);
return [x, 2 * atan2(sqrt(cosLambda * cosLambda + tanPhi0 * tanPhi0 - y * y) - cosLambda, 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) {
if (abs(x) > 1) x = sign(x) * 2 - x;
if (abs(y) > 1) y = sign(y) * 2 - 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)
: sinusoidalRaw.invert(x, y);
};
var sinuMollweide = function() {
return d3Geo.geoProjection(sinuMollweideRaw)
.rotate([-20, -55])
.scale(164.263)
.center([0, -5.4036]);
};
function homolosineRaw(lambda, phi) {
return abs(phi) > sinuMollweidePhi
? (lambda = mollweideRaw(lambda, phi), lambda[1] -= phi > 0 ? sinuMollweideY : -sinuMollweideY, lambda)
: sinusoidalRaw(lambda, phi);
}
homolosineRaw.invert = function(x, y) {
return abs(y) > sinuMollweidePhi
? mollweideRaw.invert(x, y + (y > 0 ? sinuMollweideY : -sinuMollweideY))
: sinusoidalRaw.invert(x, y);
};
var homolosine = function() {
return d3Geo.geoProjection(homolosineRaw)
.scale(152.63);
};
function pointEqual(a, b) {
return abs(a[0] - b[0]) < epsilon && abs(a[1] - b[1]) < epsilon;
}
function interpolateLine(coordinates, m) {
var i = -1,
n = coordinates.length,
p0 = coordinates[0],
p1,
dx,
dy,
resampled = [];
while (++i < n) {
p1 = coordinates[i];
dx = (p1[0] - p0[0]) / m;
dy = (p1[1] - p0[1]) / m;
for (var j = 0; j < m; ++j) resampled.push([p0[0] + j * dx, p0[1] + j * dy]);
p0 = p1;
}
resampled.push(p1);
return resampled;
}
function interpolateSphere(lobes) {
var coordinates = [],
lobe,
lambda0, phi0, phi1,
lambda2, phi2,
i, n = lobes[0].length;
// Northern Hemisphere
for (i = 0; i < n; ++i) {
lobe = lobes[0][i];
lambda0 = lobe[0][0], phi0 = lobe[0][1], phi1 = lobe[1][1];
lambda2 = lobe[2][0], phi2 = lobe[2][1];
coordinates.push(interpolateLine([
[lambda0 + epsilon, phi0 + epsilon],
[lambda0 + epsilon, phi1 - epsilon],
[lambda2 - epsilon, phi1 - epsilon],
[lambda2 - epsilon, phi2 + epsilon]
], 30));
}
// Southern Hemisphere
for (i = lobes[1].length - 1; i >= 0; --i) {
lobe = lobes[1][i];
lambda0 = lobe[0][0], phi0 = lobe[0][1], phi1 = lobe[1][1];
lambda2 = lobe[2][0], phi2 = lobe[2][1];
coordinates.push(interpolateLine([
[lambda2 - epsilon, phi2 - epsilon],
[lambda2 - epsilon, phi1 + epsilon],
[lambda0 + epsilon, phi1 + epsilon],
[lambda0 + epsilon, phi0 - epsilon]
], 30));
}
return {
type: "Polygon",
coordinates: [d3Array.merge(coordinates)]
};
}
var interrupt = function(project, lobes) {
var sphere, bounds;
function forward(lambda, phi) {
var sign$$1 = phi < 0 ? -1 : +1, lobe = lobes[+(phi < 0)];
for (var i = 0, n = lobe.length - 1; i < n && lambda > lobe[i][2][0]; ++i);
var p = project(lambda - lobe[i][1][0], phi);
p[0] += project(lobe[i][1][0], sign$$1 * phi > sign$$1 * lobe[i][0][1] ? lobe[i][0][1] : phi)[0];
return p;
}
// Assumes mutually exclusive bounding boxes for lobes.
if (project.invert) forward.invert = function(x, y) {
var bound = bounds[+(y < 0)], lobe = lobes[+(y < 0)];
for (var i = 0, n = bound.length; i < n; ++i) {
var b = bound[i];
if (b[0][0] <= x && x < b[1][0] && b[0][1] <= y && y < b[1][1]) {
var p = project.invert(x - project(lobe[i][1][0], 0)[0], y);
p[0] += lobe[i][1][0];
return pointEqual(forward(p[0], p[1]), [x, y]) ? p : null;
}
}
};
var p = d3Geo.geoProjection(forward),
stream_ = p.stream;
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, sphereStream); };
return rotateStream;
};
p.lobes = function(_) {
if (!arguments.length) return lobes.map(function(lobe) {
return lobe.map(function(l) {
return [
[l[0][0] * degrees, l[0][1] * degrees],
[l[1][0] * degrees, l[1][1] * degrees],
[l[2][0] * degrees, l[2][1] * degrees]
];
});
});
sphere = interpolateSphere(_);
lobes = _.map(function(lobe) {
return lobe.map(function(l) {
return [
[l[0][0] * radians, l[0][1] * radians],
[l[1][0] * radians, l[1][1] * radians],
[l[2][0] * radians, l[2][1] * radians]
];
});
});
bounds = lobes.map(function(lobe) {
return lobe.map(function(l) {
var x0 = project(l[0][0], l[0][1])[0],
x1 = project(l[2][0], l[2][1])[0],
y0 = project(l[1][0], l[0][1])[1],
y1 = project(l[1][0], l[1][1])[1],
t;
if (y0 > y1) t = y0, y0 = y1, y1 = t;
return [[x0, y0], [x1, y1]];
});
});
return p;
};
if (lobes != null) p.lobes(lobes);
return p;
};
var lobes = [[ // northern hemisphere
[[-180, 0], [-100, 90], [ -40, 0]],
[[ -40, 0], [ 30, 90], [ 180, 0]]
], [ // southern hemisphere
[[-180, 0], [-160, -90], [-100, 0]],
[[-100, 0], [ -60, -90], [ -20, 0]],
[[ -20, 0], [ 20, -90], [ 80, 0]],
[[ 80, 0], [ 140, -90], [ 180, 0]]
]];
var boggs$1 = function() {
return interrupt(boggsRaw, lobes)
.scale(160.857);
};
var lobes$1 = [[ // northern hemisphere
[[-180, 0], [-100, 90], [ -40, 0]],
[[ -40, 0], [ 30, 90], [ 180, 0]]
], [ // southern hemisphere
[[-180, 0], [-160, -90], [-100, 0]],
[[-100, 0], [ -60, -90], [ -20, 0]],
[[ -20, 0], [ 20, -90], [ 80, 0]],
[[ 80, 0], [ 140, -90], [ 180, 0]]
]];
var homolosine$1 = function() {
return interrupt(homolosineRaw, lobes$1)
.scale(152.63);
};
var lobes$2 = [[ // northern hemisphere
[[-180, 0], [-100, 90], [ -40, 0]],
[[ -40, 0], [ 30, 90], [ 180, 0]]
], [ // southern hemisphere
[[-180, 0], [-160, -90], [-100, 0]],
[[-100, 0], [ -60, -90], [ -20, 0]],
[[ -20, 0], [ 20, -90], [ 80, 0]],
[[ 80, 0], [ 140, -90], [ 180, 0]]
]];
var mollweide$1 = function() {
return interrupt(mollweideRaw, lobes$2)
.scale(169.529);
};
var lobes$3 = [[ // northern hemisphere
[[-180, 0], [ -90, 90], [ 0, 0]],
[[ 0, 0], [ 90, 90], [ 180, 0]]
], [ // southern hemisphere
[[-180, 0], [ -90, -90], [ 0, 0]],
[[ 0, 0], [ 90, -90], [ 180, 0]]
]];
var mollweideHemispheres = function() {
return interrupt(mollweideRaw, lobes$3)
.scale(169.529)
.rotate([20, 0]);
};
var lobes$4 = [[ // northern hemisphere
[[-180, 35], [ -30, 90], [ 0, 35]],
[[ 0, 35], [ 30, 90], [ 180, 35]]
], [ // southern hemisphere
[[-180, -10], [-102, -90], [ -65, -10]],
[[ -65, -10], [ 5, -90], [ 77, -10]],
[[ 77, -10], [ 103, -90], [ 180, -10]]
]];
var sinuMollweide$1 = function() {
return interrupt(sinuMollweideRaw, lobes$4)
.rotate([-20, -55])
.scale(164.263)
.center([0, -5.4036]);
};
var lobes$5 = [[ // northern hemisphere
[[-180, 0], [-110, 90], [ -40, 0]],
[[ -40, 0], [ 0, 90], [ 40, 0]],
[[ 40, 0], [ 110, 90], [ 180, 0]]
], [ // southern hemisphere
[[-180, 0], [-110, -90], [ -40, 0]],
[[ -40, 0], [ 0, -90], [ 40, 0]],
[[ 40, 0], [ 110, -90], [ 180, 0]]
]];
var sinusoidal$1 = function() {
return interrupt(sinusoidalRaw, lobes$5)
.scale(152.63)
.rotate([-20, 0]);
};
function kavrayskiy7Raw(lambda, phi) {
return [3 / tau * lambda * sqrt(pi * pi / 3 - phi * phi), phi];
}
kavrayskiy7Raw.invert = function(x, y) {
return [tau / 3 * x / sqrt(pi * pi / 3 - y * y), y];
};
var kavrayskiy7 = function() {
return d3Geo.geoProjection(kavrayskiy7Raw)
.scale(158.837);
};
function lagrangeRaw(n) {
function forward(lambda, phi) {
if (abs(abs(phi) - halfPi) < epsilon) return [0, phi < 0 ? -2 : 2];
var sinPhi = sin(phi),
v = pow((1 + sinPhi) / (1 - sinPhi), n / 2),
c = 0.5 * (v + 1 / v) + cos(lambda *= n);
return [
2 * sin(lambda) / c,
(v - 1 / v) / c
];
}
forward.invert = function(x, y) {
var y0 = abs(y);
if (abs(y0 - 2) < epsilon) return x ? null : [0, sign(y) * halfPi];
if (y0 > 2) return null;
x /= 2, y /= 2;
var x2 = x * x,
y2 = y * y,
t = 2 * y / (1 + x2 + y2); // tanh(nPhi)
t = pow((1 + t) / (1 - t), 1 / n);
return [
atan2(2 * x, 1 - x2 - y2) / n,
asin((t - 1) / (t + 1))
];
};
return forward;
}
var lagrange = function() {
var n = 0.5,
m = d3Geo.geoProjectionMutator(lagrangeRaw),
p = m(n);
p.spacing = function(_) {
return arguments.length ? m(n = +_) : n;
};
return p
.scale(124.75);
};
var pi_sqrt2 = pi / sqrt2;
function larriveeRaw(lambda, phi) {
return [
lambda * (1 + sqrt(cos(phi))) / 2,
phi / (cos(phi / 2) * cos(lambda / 6))
];
}
larriveeRaw.invert = function(x, y) {
var x0 = abs(x),
y0 = abs(y),
lambda = epsilon,
phi = halfPi;
if (y0 < pi_sqrt2) phi *= y0 / pi_sqrt2;
else lambda += 6 * acos(pi_sqrt2 / y0);
for (var i = 0; i < 25; i++) {
var sinPhi = sin(phi),
sqrtcosPhi = sqrt(cos(phi)),
sinPhi_2 = sin(phi / 2),
cosPhi_2 = cos(phi / 2),
sinLambda_6 = sin(lambda / 6),
cosLambda_6 = cos(lambda / 6),
f0 = 0.5 * lambda * (1 + sqrtcosPhi) - x0,
f1 = phi / (cosPhi_2 * cosLambda_6) - y0,
df0dPhi = sqrtcosPhi ? -0.25 * lambda * sinPhi / sqrtcosPhi : 0,
df0dLambda = 0.5 * (1 + sqrtcosPhi),
df1dPhi = (1 +0.5 * phi * sinPhi_2 / cosPhi_2) / (cosPhi_2 * cosLambda_6),
df1dLambda = (phi / cosPhi_2) * (sinLambda_6 / 6) / (cosLambda_6 * cosLambda_6),
denom = df0dPhi * df1dLambda - df1dPhi * df0dLambda,
dPhi = (f0 * df1dLambda - f1 * df0dLambda) / denom,
dLambda = (f1 * df0dPhi - f0 * df1dPhi) / denom;
phi -= dPhi;
lambda -= dLambda;
if (abs(dPhi) < epsilon && abs(dLambda) < epsilon) break;
}
return [x < 0 ? -lambda : lambda, y < 0 ? -phi : phi];
};
var larrivee = function() {
return d3Geo.geoProjection(larriveeRaw)
.scale(97.2672);
};
function laskowskiRaw(lambda, phi) {
var lambda2 = lambda * lambda, phi2 = phi * phi;
return [
lambda * (0.975534 + phi2 * (-0.119161 + lambda2 * -0.0143059 + phi2 * -0.0547009)),
phi * (1.00384 + lambda2 * (0.0802894 + phi2 * -0.02855 + lambda2 * 0.000199025) + phi2 * (0.0998909 + phi2 * -0.0491032))
];
}
laskowskiRaw.invert = function(x, y) {
var lambda = sign(x) * pi,
phi = y / 2,
i = 50;
do {
var lambda2 = lambda * lambda,
phi2 = phi * phi,
lambdaPhi = lambda * phi,
fx = lambda * (0.975534 + phi2 * (-0.119161 + lambda2 * -0.0143059 + phi2 * -0.0547009)) - x,
fy = phi * (1.00384 + lambda2 * (0.0802894 + phi2 * -0.02855 + lambda2 * 0.000199025) + phi2 * (0.0998909 + phi2 * -0.0491032)) - y,
deltaxDeltaLambda = 0.975534 - phi2 * (0.119161 + 3 * lambda2 * 0.0143059 + phi2 * 0.0547009),
deltaxDeltaPhi = -lambdaPhi * (2 * 0.119161 + 4 * 0.0547009 * phi2 + 2 * 0.0143059 * lambda2),
deltayDeltaLambda = lambdaPhi * (2 * 0.0802894 + 4 * 0.000199025 * lambda2 + 2 * -0.02855 * phi2),
deltayDeltaPhi = 1.00384 + lambda2 * (0.0802894 + 0.000199025 * lambda2) + phi2 * (3 * (0.0998909 - 0.02855 * lambda2) - 5 * 0.0491032 * phi2),
denominator = deltaxDeltaPhi * deltayDeltaLambda - deltayDeltaPhi * deltaxDeltaLambda,
deltaLambda = (fy * deltaxDeltaPhi - fx * deltayDeltaPhi) / denominator,
deltaPhi = (fx * deltayDeltaLambda - fy * deltaxDeltaLambda) / denominator;
lambda -= deltaLambda, phi -= deltaPhi;
} while ((abs(deltaLambda) > epsilon || abs(deltaPhi) > epsilon) && --i > 0);
return i && [lambda, phi];
};
var laskowski = function() {
return d3Geo.geoProjection(laskowskiRaw)
.scale(139.98);
};
function littrowRaw(lambda, phi) {
return [
sin(lambda) / cos(phi),
tan(phi) * cos(lambda)
];
}
littrowRaw.invert = function(x, y) {
var x2 = x * x,
y2 = y * y,
y2_1 = y2 + 1,
cosPhi = x
? sqrt1_2 * sqrt((y2_1 - sqrt(x2 * x2 + 2 * x2 * (y2 - 1) + y2_1 * y2_1)) / x2 + 1)
: 1 / sqrt(y2_1);
return [
asin(x * cosPhi),
sign(y) * acos(cosPhi)
];
};
var littrow = function() {
return d3Geo.geoProjection(littrowRaw)
.scale(144.049)
.clipAngle(90 - 1e-3);
};
function loximuthalRaw(phi0) {
var cosPhi0 = cos(phi0),
tanPhi0 = tan(quarterPi + phi0 / 2);
function forward(lambda, phi) {
var y = phi - phi0,
x = abs(y) < epsilon ? lambda * cosPhi0
: abs(x = quarterPi + phi / 2) < epsilon || abs(abs(x) - halfPi) < epsilon
? 0 : lambda * y / log(tan(x) / tanPhi0);
return [x, y];
}
forward.invert = function(x, y) {
var lambda,
phi = y + phi0;
return [
abs(y) < epsilon ? x / cosPhi0
: (abs(lambda = quarterPi + phi / 2) < epsilon || abs(abs(lambda) - halfPi) < epsilon) ? 0
: x * log(tan(lambda) / tanPhi0) / y,
phi
];
};
return forward;
}
var loximuthal = function() {
return parallel1(loximuthalRaw)
.parallel(40)
.scale(158.837);
};
function millerRaw(lambda, phi) {
return [lambda, 1.25 * log(tan(quarterPi + 0.4 * phi))];
}
millerRaw.invert = function(x, y) {
return [x, 2.5 * atan(exp(0.8 * y)) - 0.625 * pi];
};
var miller = function() {
return d3Geo.geoProjection(millerRaw)
.scale(108.318);
};
function modifiedStereographicRaw(C) {
var m = C.length - 1;
function forward(lambda, phi) {
var cosPhi = cos(phi),
k = 2 / (1 + cosPhi * cos(lambda)),
zr = k * cosPhi * sin(lambda),
zi = k * sin(phi),
i = m,
w = C[i],
ar = w[0],
ai = w[1],
t;
while (--i >= 0) {
w = C[i];
ar = w[0] + zr * (t = ar) - zi * ai;
ai = w[1] + zr * ai + zi * t;
}
ar = zr * (t = ar) - zi * ai;
ai = zr * ai + zi * t;
return [ar, ai];
}
forward.invert = function(x, y) {
var i = 20,
zr = x,
zi = y;
do {
var j = m,
w = C[j],
ar = w[0],
ai = w[1],
br = 0,
bi = 0,
t;
while (--j >= 0) {
w = C[j];
br = ar + zr * (t = br) - zi * bi;
bi = ai + zr * bi + zi * t;
ar = w[0] + zr * (t = ar) - zi * ai;
ai = w[1] + zr * ai + zi * t;
}
br = ar + zr * (t = br) - zi * bi;
bi = ai + zr * bi + zi * t;
ar = zr * (t = ar) - zi * ai - x;
ai = zr * ai + zi * t - y;
var denominator = br * br + bi * bi, deltar, deltai;
zr -= deltar = (ar * br + ai * bi) / denominator;
zi -= deltai = (ai * br - ar * bi) / denominator;
} while (abs(deltar) + abs(deltai) > epsilon * epsilon && --i > 0);
if (i) {
var rho = sqrt(zr * zr + zi * zi),
c = 2 * atan(rho * 0.5),
sinc = sin(c);
return [atan2(zr * sinc, rho * cos(c)), rho ? asin(zi * sinc / rho) : 0];
}
};
return forward;
}
var alaska = [[0.9972523, 0], [0.0052513, -0.0041175], [0.0074606, 0.0048125], [-0.0153783, -0.1968253], [0.0636871, -0.1408027], [0.3660976, -0.2937382]];
var gs48 = [[0.98879, 0], [0, 0], [-0.050909, 0], [0, 0], [0.075528, 0]];
var gs50 = [[0.9842990, 0], [0.0211642, 0.0037608], [-0.1036018, -0.0575102], [-0.0329095, -0.0320119], [0.0499471, 0.1223335], [0.0260460, 0.0899805], [0.0007388, -0.1435792], [0.0075848, -0.1334108], [-0.0216473, 0.0776645], [-0.0225161, 0.0853673]];
var miller$1 = [[0.9245, 0], [0, 0], [0.01943, 0]];
var lee = [[0.721316, 0], [0, 0], [-0.00881625, -0.00617325]];
function modifiedStereographicAlaska() {
return modifiedStereographic(alaska, [152, -64])
.scale(1500)
.center([-160.908, 62.4864])
.clipAngle(25);
}
function modifiedStereographicGs48() {
return modifiedStereographic(gs48, [95, -38])
.scale(1000)
.clipAngle(55)
.center([-96.5563, 38.8675]);
}
function modifiedStereographicGs50() {
return modifiedStereographic(gs50, [120, -45])
.scale(359.513)
.clipAngle(55)
.center([-117.474, 53.0628]);
}
function modifiedStereographicMiller() {
return modifiedStereographic(miller$1, [-20, -18])
.scale(209.091)
.center([20, 16.7214])
.clipAngle(82);
}
function modifiedStereographicLee() {
return modifiedStereographic(lee, [165, 10])
.scale(250)
.clipAngle(130)
.center([-165, -10]);
}
function modifiedStereographic(coefficients, rotate) {
var p = d3Geo.geoProjection(modifiedStereographicRaw(coefficients)).rotate(rotate).clipAngle(90),
r = d3Geo.geoRotation(rotate),
center = p.center;
delete p.rotate;
p.center = function(_) {
return arguments.length ? center(r(_)) : r.invert(center());
};
return p;
}
var sqrt6 = sqrt(6);
var sqrt7 = sqrt(7);
function mtFlatPolarParabolicRaw(lambda, phi) {
var theta = asin(7 * sin(phi) / (3 * sqrt6));
return [
sqrt6 * lambda * (2 * cos(2 * theta / 3) - 1) / sqrt7,
9 * sin(theta / 3) / sqrt7
];
}
mtFlatPolarParabolicRaw.invert = function(x, y) {
var theta = 3 * asin(y * sqrt7 / 9);
return [
x * sqrt7 / (sqrt6 * (2 * cos(2 * theta / 3) - 1)),
asin(sin(theta) * 3 * sqrt6 / 7)
];
};
var mtFlatPolarParabolic = function() {
return d3Geo.geoProjection(mtFlatPolarParabolicRaw)
.scale(164.859);
};
function mtFlatPolarQuarticRaw(lambda, phi) {
var k = (1 + sqrt1_2) * sin(phi),
theta = phi;
for (var i = 0, delta; i < 25; i++) {
theta -= delta = (sin(theta / 2) + sin(theta) - k) / (0.5 * cos(theta / 2) + cos(theta));
if (abs(delta) < epsilon) break;
}
return [
lambda * (1 + 2 * cos(theta) / cos(theta / 2)) / (3 * sqrt2),
2 * sqrt(3) * sin(theta / 2) / sqrt(2 + sqrt2)
];
}
mtFlatPolarQuarticRaw.invert = function(x, y) {
var sinTheta_2 = y * sqrt(2 + sqrt2) / (2 * sqrt(3)),
theta = 2 * asin(sinTheta_2);
return [
3 * sqrt2 * x / (1 + 2 * cos(theta) / cos(theta / 2)),
asin((sinTheta_2 + sin(theta)) / (1 + sqrt1_2))
];
};
var mtFlatPolarQuartic = function() {
return d3Geo.geoProjection(mtFlatPolarQuarticRaw)
.scale(188.209);
};
function mtFlatPolarSinusoidalRaw(lambda, phi) {
var A = sqrt(6 / (4 + pi)),
k = (1 + pi / 4) * sin(phi),
theta = phi / 2;
for (var i = 0, delta; i < 25; i++) {
theta -= delta = (theta / 2 + sin(theta) - k) / (0.5 + cos(theta));
if (abs(delta) < epsilon) break;
}
return [
A * (0.5 + cos(theta)) * lambda / 1.5,
A * theta
];
}
mtFlatPolarSinusoidalRaw.invert = function(x, y) {
var A = sqrt(6 / (4 + pi)),
theta = y / A;
if (abs(abs(theta) - halfPi) < epsilon) theta = theta < 0 ? -halfPi : halfPi;
return [
1.5 * x / (A * (0.5 + cos(theta))),
asin((theta / 2 + sin(theta)) / (1 + pi / 4))
];
};
var mtFlatPolarSinusoidal = function() {
return d3Geo.geoProjection(mtFlatPolarSinusoidalRaw)
.scale(166.518);
};
function naturalEarthRaw(lambda, phi) {
var phi2 = phi * phi, phi4 = phi2 * phi2;
return [
lambda * (0.8707 - 0.131979 * phi2 + phi4 * (-0.013791 + phi4 * (0.003971 * phi2 - 0.001529 * phi4))),
phi * (1.007226 + phi2 * (0.015085 + phi4 * (-0.044475 + 0.028874 * phi2 - 0.005916 * phi4)))
];
}
naturalEarthRaw.invert = function(x, y) {
var phi = y, i = 25, delta;
do {
var phi2 = phi * phi, phi4 = phi2 * phi2;
phi -= delta = (phi * (1.007226 + phi2 * (0.015085 + phi4 * (-0.044475 + 0.028874 * phi2 - 0.005916 * phi4))) - y) /
(1.007226 + phi2 * (0.015085 * 3 + phi4 * (-0.044475 * 7 + 0.028874 * 9 * phi2 - 0.005916 * 11 * phi4)));
} while (abs(delta) > epsilon && --i > 0);
return [
x / (0.8707 + (phi2 = phi * phi) * (-0.131979 + phi2 * (-0.013791 + phi2 * phi2 * phi2 * (0.003971 - 0.001529 * phi2)))),
phi
];
};
var naturalEarth = function() {
return d3Geo.geoProjection(naturalEarthRaw)
.scale(175.295);
};
function naturalEarth2Raw(lambda, phi) {
var phi2 = phi * phi, phi4 = phi2 * phi2, phi6 = phi2 * phi4;
return [
lambda * (0.84719 - 0.13063 * phi2 + phi6 * phi6 * (-0.04515 + 0.05494 * phi2 - 0.02326 * phi4 + 0.00331 * phi6)),
phi * (1.01183 + phi4 * phi4 * (-0.02625 + 0.01926 * phi2 - 0.00396 * phi4))
];
}
naturalEarth2Raw.invert = function(x, y) {
var phi = y, i = 25, delta, phi2, phi4, phi6;
do {
phi2 = phi * phi; phi4 = phi2 * phi2;
phi -= delta = ((phi * (1.01183 + phi4 * phi4 * (-0.02625 + 0.01926 * phi2 - 0.00396 * phi4))) - y) /
(1.01183 + phi4 * phi4 * ((9 * -0.02625) + (11 * 0.01926) * phi2 + (13 * -0.00396) * phi4));
} while (abs(delta) > epsilon2 && --i > 0);
phi2 = phi * phi; phi4 = phi2 * phi2; phi6 = phi2 * phi4;
return [
x / (0.84719 - 0.13063 * phi2 + phi6 * phi6 * (-0.04515 + 0.05494 * phi2 - 0.02326 * phi4 + 0.00331 * phi6)),
phi
];
};
var naturalEarth2 = function() {
return d3Geo.geoProjection(naturalEarth2Raw)
.scale(175.295);
};
function nellHammerRaw(lambda, phi) {
return [
lambda * (1 + cos(phi)) / 2,
2 * (phi - tan(phi / 2))
];
}
nellHammerRaw.invert = function(x, y) {
var p = y / 2;
for (var i = 0, delta = Infinity; i < 10 && abs(delta) > epsilon; ++i) {
var c = cos(y / 2);
y -= delta = (y - tan(y / 2) - p) / (1 - 0.5 / (c * c));
}
return [
2 * x / (1 + cos(y)),
y
];
};
var nellHammer = function() {
return d3Geo.geoProjection(nellHammerRaw)
.scale(152.63);
};
// Based on Java implementation by Bojan Savric.
// https://github.com/OSUCartography/JMapProjLib/blob/master/src/com/jhlabs/map/proj/PattersonProjection.java
var pattersonK1 = 1.0148;
var pattersonK2 = 0.23185;
var pattersonK3 = -0.14499;
var pattersonK4 = 0.02406;
var pattersonC1 = pattersonK1;
var pattersonC2 = 5 * pattersonK2;
var pattersonC3 = 7 * pattersonK3;
var pattersonC4 = 9 * pattersonK4;
var pattersonYmax = 1.790857183;
function pattersonRaw(lambda, phi) {
var phi2 = phi * phi;
return [
lambda,
phi * (pattersonK1 + phi2 * phi2 * (pattersonK2 + phi2 * (pattersonK3 + pattersonK4 * phi2)))
];
}
pattersonRaw.invert = function(x, y) {
if (y > pattersonYmax) y = pattersonYmax;
else if (y < -pattersonYmax) y = -pattersonYmax;
var yc = y, delta;
do { // Newton-Raphson
var y2 = yc * yc;
yc -= delta = ((yc * (pattersonK1 + y2 * y2 * (pattersonK2 + y2 * (pattersonK3 + pattersonK4 * y2)))) - y) / (pattersonC1 + y2 * y2 * (pattersonC2 + y2 * (pattersonC3 + pattersonC4 * y2)));
} while (abs(delta) > epsilon);
return [x, yc];
};
var patterson = function() {
return d3Geo.geoProjection(pattersonRaw)
.scale(139.319);
};
function polyconicRaw(lambda, phi) {
if (abs(phi) < epsilon) return [lambda, 0];
var tanPhi = tan(phi),
k = lambda * sin(phi);
return [
sin(k) / tanPhi,
phi + (1 - cos(k)) / tanPhi
];
}
polyconicRaw.invert = function(x, y) {
if (abs(y) < epsilon) return [x, 0];
var k = x * x + y * y,
phi = y * 0.5,
i = 10, delta;
do {
var tanPhi = tan(phi),
secPhi = 1 / cos(phi),
j = k - 2 * y * phi + phi * phi;
phi -= delta = (tanPhi * j + 2 * (phi - y)) / (2 + j * secPhi * secPhi + 2 * (phi - y) * tanPhi);
} while (abs(delta) > epsilon && --i > 0);
tanPhi = tan(phi);
return [
(abs(y) < abs(phi + 1 / tanPhi) ? asin(x * tanPhi) : sign(x) * (acos(abs(x * tanPhi)) + halfPi)) / sin(phi),
phi
];
};
var polyconic = function() {
return d3Geo.geoProjection(polyconicRaw)
.scale(103.74);
};
// Note: 6-element arrays are used to denote the 3x3 affine transform matrix:
// [a, b, c,
// d, e, f,
// 0, 0, 1] - this redundant row is left out.
// Transform matrix for [a0, a1] -> [b0, b1].
var matrix = function(a, b) {
var u = subtract(a[1], a[0]),
v = subtract(b[1], b[0]),
phi = angle$1(u, v),
s = length(u) / length(v);
return multiply([
1, 0, a[0][0],
0, 1, a[0][1]
], multiply([
s, 0, 0,
0, s, 0
], multiply([
cos(phi), sin(phi), 0,
-sin(phi), cos(phi), 0
], [
1, 0, -b[0][0],
0, 1, -b[0][1]
])));
};
// Inverts a transform matrix.
function inverse(m) {
var k = 1 / (m[0] * m[4] - m[1] * m[3]);
return [
k * m[4], -k * m[1], k * (m[1] * m[5] - m[2] * m[4]),
-k * m[3], k * m[0], k * (m[2] * m[3] - m[0] * m[5])
];
}
// Multiplies two 3x2 matrices.
function multiply(a, b) {
return [
a[0] * b[0] + a[1] * b[3],
a[0] * b[1] + a[1] * b[4],
a[0] * b[2] + a[1] * b[5] + a[2],
a[3] * b[0] + a[4] * b[3],
a[3] * b[1] + a[4] * b[4],
a[3] * b[2] + a[4] * b[5] + a[5]
];
}
// Subtracts 2D vectors.
function subtract(a, b) {
return [a[0] - b[0], a[1] - b[1]];
}
// Magnitude of a 2D vector.
function length(v) {
return sqrt(v[0] * v[0] + v[1] * v[1]);
}
// Angle between two 2D vectors.
function angle$1(a, b) {
return atan2(a[0] * b[1] - a[1] * b[0], a[0] * b[0] + a[1] * b[1]);
}
// Creates a polyhedral projection.
// * root: a spanning tree of polygon faces. Nodes are automatically
// augmented with a transform matrix.
// * face: a function that returns the appropriate node for a given {lambda, phi}
// point (radians).
// * r: rotation angle for final polyhedral net. Defaults to -pi / 6 (for
// butterflies).
var polyhedral = function(root, face, r) {
r = r == null ? -pi / 6 : r; // TODO automate
recurse(root, {transform: [
cos(r), sin(r), 0,
-sin(r), cos(r), 0
]});
function recurse(node, parent) {
node.edges = faceEdges(node.face);
// Find shared edge.
if (parent.face) {
var shared = node.shared = sharedEdge(node.face, parent.face),
m = matrix(shared.map(parent.project), shared.map(node.project));
node.transform = parent.transform ? multiply(parent.transform, m) : m;
// Replace shared edge in parent edges array.
var edges = parent.edges;
for (var i = 0, n = edges.length; i < n; ++i) {
if (pointEqual$1(shared[0], edges[i][1]) && pointEqual$1(shared[1], edges[i][0])) edges[i] = node;
if (pointEqual$1(shared[0], edges[i][0]) && pointEqual$1(shared[1], edges[i][1])) edges[i] = node;
}
edges = node.edges;
for (i = 0, n = edges.length; i < n; ++i) {
if (pointEqual$1(shared[0], edges[i][0]) && pointEqual$1(shared[1], edges[i][1])) edges[i] = parent;
if (pointEqual$1(shared[0], edges[i][1]) && pointEqual$1(shared[1], edges[i][0])) edges[i] = parent;
}
} else {
node.transform = parent.transform;
}
if (node.children) {
node.children.forEach(function(child) {
recurse(child, node);
});
}
return node;
}
function forward(lambda, phi) {
var node = face(lambda, phi),
point = node.project([lambda * degrees, phi * degrees]),
t;
if (t = node.transform) {
return [
t[0] * point[0] + t[1] * point[1] + t[2],
-(t[3] * point[0] + t[4] * point[1] + t[5])
];
}
point[1] = -point[1];
return point;
}
// Naive inverse! A faster solution would use bounding boxes, or even a
// polygonal quadtree.
if (hasInverse(root)) forward.invert = function(x, y) {
var coordinates = faceInvert(root, [x, -y]);
return coordinates && (coordinates[0] *= radians, coordinates[1] *= radians, coordinates);
};
function faceInvert(node, coordinates) {
var invert = node.project.invert,
t = node.transform,
point = coordinates;
if (t) {
t = inverse(t);
point = [
t[0] * point[0] + t[1] * point[1] + t[2],
(t[3] * point[0] + t[4] * point[1] + t[5])
];
}
if (invert && node === faceDegrees(p = invert(point))) return p;
var p,
children = node.children;
for (var i = 0, n = children && children.length; i < n; ++i) {
if (p = faceInvert(children[i], coordinates)) return p;
}
}
function faceDegrees(coordinates) {
return face(coordinates[0] * radians, coordinates[1] * radians);
}
var proj = d3Geo.geoProjection(forward),
stream_ = proj.stream;
// run around the mesh of faces and stream all vertices in order to create the clip polygon
var polygon = [];
outline(1e-3, {point: function(lambda, phi) { polygon.push([lambda, phi]); }}, root);
polygon.push(polygon[0]);
if (proj.clipPolygon) proj.clipPolygon([ polygon ]);
proj.stream = function(stream) {
var rotate = proj.rotate(),
clipPolygon = proj.clipPolygon ? proj.clipPolygon() : null,
rotateStream = stream_(stream),
sphereStream = ((clipPolygon ? proj.clipPolygon([]) : proj).rotate([0, 0]), stream_(stream));
proj.rotate(rotate);
if (clipPolygon) proj.clipPolygon(clipPolygon);
rotateStream.sphere = function() {
sphereStream.polygonStart();
sphereStream.lineStart();
outline(1e-5, sphereStream, root);
sphereStream.lineEnd();
sphereStream.polygonEnd();
};
return rotateStream;
};
return proj;
};
function outline(eps, stream, node, parent) {
var point,
edges = node.edges,
n = edges.length,
edge,
multiPoint = {type: "MultiPoint", coordinates: node.face},
notPoles = node.face.filter(function(d) { return abs(d[1]) !== 90; }),
b = d3Geo.geoBounds({type: "MultiPoint", coordinates: notPoles}),
inside = false,
j = -1,
dx = b[1][0] - b[0][0];
// TODO
var c = dx === 180 || dx === 360
? [(b[0][0] + b[1][0]) / 2, (b[0][1] + b[1][1]) / 2]
: d3Geo.geoCentroid(multiPoint);
// First find the shared edge…
if (parent) while (++j < n) {
if (edges[j] === parent) break;
}
++j;
for (var i = 0; i < n; ++i) {
edge = edges[(i + j) % n];
if (Array.isArray(edge)) {
if (!inside) {
stream.point((point = d3Geo.geoInterpolate(edge[0], c)(eps))[0], point[1]);
inside = true;
}
stream.point((point = d3Geo.geoInterpolate(edge[1], c)(eps))[0], point[1]);
} else {
inside = false;
if (edge !== parent) outline(eps, stream, edge, node);
}
}
}
// Tests equality of two spherical points.
function pointEqual$1(a, b) {
return a && b && a[0] === b[0] && a[1] === b[1];
}
// Finds a shared edge given two clockwise polygons.
function sharedEdge(a, b) {
var x, y, n = a.length, found = null;
for (var i = 0; i < n; ++i) {
x = a[i];
for (var j = b.length; --j >= 0;) {
y = b[j];
if (x[0] === y[0] && x[1] === y[1]) {
if (found) return [found, x];
found = x;
}
}
}
}
// Converts an array of n face vertices to an array of n + 1 edges.
function faceEdges(face) {
var n = face.length,
edges = [];
for (var a = face[n - 1], i = 0; i < n; ++i) edges.push([a, a = face[i]]);
return edges;
}
function hasInverse(node) {
return node.project.invert || node.children && node.children.some(hasInverse);
}
// TODO generate on-the-fly to avoid external modification.
var octahedron = [
[0, 90],
[-90, 0], [0, 0], [90, 0], [180, 0],
[0, -90]
];
var octahedron$1 = [
[0, 2, 1],
[0, 3, 2],
[5, 1, 2],
[5, 2, 3],
[0, 1, 4],
[0, 4, 3],
[5, 4, 1],
[5, 3, 4]
].map(function(face) {
return face.map(function(i) {
return octahedron[i];
});
});
var butterfly = function(faceProjection) {
faceProjection = faceProjection || function(face) {
var c = d3Geo.geoCentroid({type: "MultiPoint", coordinates: face});
return d3Geo.geoGnomonic().scale(1).translate([0, 0]).rotate([-c[0], -c[1]]);
};
var faces = octahedron$1.map(function(face) {
return {face: face, project: faceProjection(face)};
});
[-1, 0, 0, 1, 0, 1, 4, 5].forEach(function(d, i) {
var node = faces[d];
node && (node.children || (node.children = [])).push(faces[i]);
});
return polyhedral(faces[0], function(lambda, phi) {
return faces[lambda < -pi / 2 ? phi < 0 ? 6 : 4
: lambda < 0 ? phi < 0 ? 2 : 0
: lambda < pi / 2 ? phi < 0 ? 3 : 1
: phi < 0 ? 7 : 5];
})
.scale(101.858)
.center([0, 45]);
};
var kx = 2 / sqrt(3);
function collignonK(a, b) {
var p = collignonRaw(a, b);
return [p[0] * kx, p[1]];
}
collignonK.invert = function(x,y) {
return collignonRaw.invert(x / kx, y);
};
var collignon$1 = function(faceProjection) {
faceProjection = faceProjection || function(face) {
var c = d3Geo.geoCentroid({type: "MultiPoint", coordinates: face});
return d3Geo.geoProjection(collignonK).translate([0, 0]).scale(1).rotate(c[1] > 0 ? [-c[0], 0] : [180 - c[0], 180]);
};
var faces = octahedron$1.map(function(face) {
return {face: face, project: faceProjection(face)};
});
[-1, 0, 0, 1, 0, 1, 4, 5].forEach(function(d, i) {
var node = faces[d];
node && (node.children || (node.children = [])).push(faces[i]);
});
return polyhedral(faces[0], function(lambda, phi) {
return faces[lambda < -pi / 2 ? phi < 0 ? 6 : 4
: lambda < 0 ? phi < 0 ? 2 : 0
: lambda < pi / 2 ? phi < 0 ? 3 : 1
: phi < 0 ? 7 : 5];
})
.scale(121.906)
.center([0, 48.5904]);
};
var waterman = function(faceProjection) {
faceProjection = faceProjection || function(face) {
var c = face.length === 6 ? d3Geo.geoCentroid({type: "MultiPoint", coordinates: face}) : face[0];
return d3Geo.geoGnomonic().scale(1).translate([0, 0]).rotate([-c[0], -c[1]]);
};
var w5 = octahedron$1.map(function(face) {
var xyz = face.map(cartesian),
n = xyz.length,
a = xyz[n - 1],
b,
hexagon = [];
for (var i = 0; i < n; ++i) {
b = xyz[i];
hexagon.push(spherical([
a[0] * 0.9486832980505138 + b[0] * 0.31622776601683794,
a[1] * 0.9486832980505138 + b[1] * 0.31622776601683794,
a[2] * 0.9486832980505138 + b[2] * 0.31622776601683794
]), spherical([
b[0] * 0.9486832980505138 + a[0] * 0.31622776601683794,
b[1] * 0.9486832980505138 + a[1] * 0.31622776601683794,
b[2] * 0.9486832980505138 + a[2] * 0.31622776601683794
]));
a = b;
}
return hexagon;
});
var cornerNormals = [];
var parents = [-1, 0, 0, 1, 0, 1, 4, 5];
w5.forEach(function(hexagon, j) {
var face = octahedron$1[j],
n = face.length,
normals = cornerNormals[j] = [];
for (var i = 0; i < n; ++i) {
w5.push([
face[i],
hexagon[(i * 2 + 2) % (2 * n)],
hexagon[(i * 2 + 1) % (2 * n)]
]);
parents.push(j);
normals.push(cross(
cartesian(hexagon[(i * 2 + 2) % (2 * n)]),
cartesian(hexagon[(i * 2 + 1) % (2 * n)])
));
}
});
var faces = w5.map(function(face) {
return {
project: faceProjection(face),
face: face
};
});
parents.forEach(function(d, i) {
var parent = faces[d];
parent && (parent.children || (parent.children = [])).push(faces[i]);
});
function face(lambda, phi) {
var cosphi = cos(phi),
p = [cosphi * cos(lambda), cosphi * sin(lambda), sin(phi)];
var hexagon = lambda < -pi / 2 ? phi < 0 ? 6 : 4
: lambda < 0 ? phi < 0 ? 2 : 0
: lambda < pi / 2 ? phi < 0 ? 3 : 1
: phi < 0 ? 7 : 5;
var n = cornerNormals[hexagon];
return faces[dot(n[0], p) < 0 ? 8 + 3 * hexagon
: dot(n[1], p) < 0 ? 8 + 3 * hexagon + 1
: dot(n[2], p) < 0 ? 8 + 3 * hexagon + 2
: hexagon];
}
return polyhedral(faces[0], face)
.scale(110.625)
.center([0,45]);
};
function dot(a, b) {
for (var i = 0, n = a.length, s = 0; i < n; ++i) s += a[i] * b[i];
return s;
}
function cross(a, b) {
return [
a[1] * b[2] - a[2] * b[1],
a[2] * b[0] - a[0] * b[2],
a[0] * b[1] - a[1] * b[0]
];
}
// Converts 3D Cartesian to spherical coordinates (degrees).
function spherical(cartesian) {
return [
atan2(cartesian[1], cartesian[0]) * degrees,
asin(max(-1, min(1, cartesian[2]))) * degrees
];
}
// Converts spherical coordinates (degrees) to 3D Cartesian.
function cartesian(coordinates) {
var lambda = coordinates[0] * radians,
phi = coordinates[1] * radians,
cosphi = cos(phi);
return [
cosphi * cos(lambda),
cosphi * sin(lambda),
sin(phi)
];
}
var noop = function() {};
var clockwise = function(ring) {
if ((n = ring.length) < 4) return false;
var i = 0,
n,
area = ring[n - 1][1] * ring[0][0] - ring[n - 1][0] * ring[0][1];
while (++i < n) area += ring[i - 1][1] * ring[i][0] - ring[i - 1][0] * ring[i][1];
return area <= 0;
};
var contains = function(ring, point) {
var x = point[0],
y = point[1],
contains = false;
for (var i = 0, n = ring.length, j = n - 1; i < n; j = i++) {
var pi = ring[i], xi = pi[0], yi = pi[1],
pj = ring[j], xj = pj[0], yj = pj[1];
if (((yi > y) ^ (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi)) contains = !contains;
}
return contains;
};
var index = function(object, projection) {
var stream = projection.stream, project;
if (!stream) throw new Error("invalid projection");
switch (object && object.type) {
case "Feature": project = projectFeature; break;
case "FeatureCollection": project = projectFeatureCollection; break;
default: project = projectGeometry; break;
}
return project(object, stream);
};
function projectFeatureCollection(o, stream) {
return {
type: "FeatureCollection",
features: o.features.map(function(f) {
return projectFeature(f, stream);
})
};
}
function projectFeature(o, stream) {
return {
type: "Feature",
id: o.id,
properties: o.properties,
geometry: projectGeometry(o.geometry, stream)
};
}
function projectGeometryCollection(o, stream) {
return {
type: "GeometryCollection",
geometries: o.geometries.map(function(o) {
return projectGeometry(o, stream);
})
};
}
function projectGeometry(o, stream) {
if (!o) return null;
if (o.type === "GeometryCollection") return projectGeometryCollection(o, stream);
var sink;
switch (o.type) {
case "Point": sink = sinkPoint; break;
case "MultiPoint": sink = sinkPoint; break;
case "LineString": sink = sinkLine; break;
case "MultiLineString": sink = sinkLine; break;
case "Polygon": sink = sinkPolygon; break;
case "MultiPolygon": sink = sinkPolygon; break;
case "Sphere": sink = sinkPolygon; break;
default: return null;
}
d3Geo.geoStream(o, stream(sink));
return sink.result();
}
var points = [];
var lines = [];
var sinkPoint = {
point: function(x, y) {
points.push([x, y]);
},
result: function() {
var result = !points.length ? null
: points.length < 2 ? {type: "Point", coordinates: points[0]}
: {type: "MultiPoint", coordinates: points};
points = [];
return result;
}
};
var sinkLine = {
lineStart: noop,
point: function(x, y) {
points.push([x, y]);
},
lineEnd: function() {
if (points.length) lines.push(points), points = [];
},
result: function() {
var result = !lines.length ? null
: lines.length < 2 ? {type: "LineString", coordinates: lines[0]}
: {type: "MultiLineString", coordinates: lines};
lines = [];
return result;
}
};
var sinkPolygon = {
polygonStart: noop,
lineStart: noop,
point: function(x, y) {
points.push([x, y]);
},
lineEnd: function() {
var n = points.length;
if (n) {
do points.push(points[0].slice()); while (++n < 4);
lines.push(points), points = [];
}
},
polygonEnd: noop,
result: function() {
if (!lines.length) return null;
var polygons = [],
holes = [];
// https://github.com/d3/d3/issues/1558
lines.forEach(function(ring) {
if (clockwise(ring)) polygons.push([ring]);
else holes.push(ring);
});
holes.forEach(function(hole) {
var point = hole[0];
polygons.some(function(polygon) {
if (contains(polygon[0], point)) {
polygon.push(hole);
return true;
}
}) || polygons.push([hole]);
});
lines = [];
return !polygons.length ? null
: polygons.length > 1 ? {type: "MultiPolygon", coordinates: polygons}
: {type: "Polygon", coordinates: polygons[0]};
}
};
var quincuncial = function(project) {
var dx = project(halfPi, 0)[0] - project(-halfPi, 0)[0];
function projectQuincuncial(lambda, phi) {
var t = abs(lambda) < halfPi,
p = project(t ? lambda : lambda > 0 ? lambda - pi : lambda + pi, phi),
x = (p[0] - p[1]) * sqrt1_2,
y = (p[0] + p[1]) * sqrt1_2;
if (t) return [x, y];
var d = dx * sqrt1_2,
s = x > 0 ^ y > 0 ? -1 : 1;
return [s * x - sign(y) * d, s * y - sign(x) * d];
}
if (project.invert) projectQuincuncial.invert = function(x0, y0) {
var x = (x0 + y0) * sqrt1_2,
y = (y0 - x0) * sqrt1_2,
t = abs(x) < 0.5 * dx && abs(y) < 0.5 * dx;
if (!t) {
var d = dx * sqrt1_2,
s = x > 0 ^ y > 0 ? -1 : 1,
x1 = -s * x0 + (y > 0 ? 1 : -1) * d,
y1 = -s * y0 + (x > 0 ? 1 : -1) * d;
x = (-x1 - y1) * sqrt1_2;
y = (x1 - y1) * sqrt1_2;
}
var p = project.invert(x, y);
if (!t) p[0] += x > 0 ? pi : -pi;
return p;
};
return d3Geo.geoProjection(projectQuincuncial)
.rotate([-90, -90, 45])
.clipAngle(180 - 1e-3);
};
var gringorten$1 = function() {
return quincuncial(gringortenRaw)
.scale(176.423);
};
var peirce = function() {
return quincuncial(guyouRaw)
.scale(111.48);
};
var quantize = function(input, digits) {
if (!(0 <= (digits = +digits) && digits <= 20)) throw new Error("invalid digits");
function quantizePoint(input) {
var n = input.length, i = 2, output = new Array(n);
output[0] = +input[0].toFixed(digits);
output[1] = +input[1].toFixed(digits);
while (i < n) output[i] = input[i], ++i;
return output;
}
function quantizePoints(input) {
return input.map(quantizePoint);
}
function quantizePolygon(input) {
return input.map(quantizePoints);
}
function quantizeGeometry(input) {
if (input == null) return input;
var output;
switch (input.type) {
case "GeometryCollection": output = {type: "GeometryCollection", geometries: input.geometries.map(quantizeGeometry)}; break;
case "Point": output = {type: "Point", coordinates: quantizePoint(input.coordinates)}; break;
case "MultiPoint": case "LineString": output = {type: input.type, coordinates: quantizePoints(input.coordinates)}; break;
case "MultiLineString": case "Polygon": output = {type: input.type, coordinates: quantizePolygon(input.coordinates)}; break;
case "MultiPolygon": output = {type: "MultiPolygon", coordinates: input.coordinates.map(quantizePolygon)}; break;
default: return input;
}
if (input.bbox != null) output.bbox = input.bbox;
return output;
}
function quantizeFeature(input) {
var output = {type: "Feature", properties: input.properties, geometry: quantizeGeometry(input.geometry)};
if (input.id != null) output.id = input.id;
if (input.bbox != null) output.bbox = input.bbox;
return output;
}
if (input != null) switch (input.type) {
case "Feature": return quantizeFeature(input);
case "FeatureCollection": {
var output = {type: "FeatureCollection", features: input.features.map(quantizeFeature)};
if (input.bbox != null) output.bbox = input.bbox;
return output;
}
default: return quantizeGeometry(input);
}
return input;
};
function rectangularPolyconicRaw(phi0) {
var sinPhi0 = sin(phi0);
function forward(lambda, phi) {
var A = sinPhi0 ? tan(lambda * sinPhi0 / 2) / sinPhi0 : lambda / 2;
if (!phi) return [2 * A, -phi0];
var E = 2 * atan(A * sin(phi)),
cotPhi = 1 / tan(phi);
return [
sin(E) * cotPhi,
phi + (1 - cos(E)) * cotPhi - phi0
];
}
// TODO return null for points outside outline.
forward.invert = function(x, y) {
if (abs(y += phi0) < epsilon) return [sinPhi0 ? 2 * atan(sinPhi0 * x / 2) / sinPhi0 : x, 0];
var k = x * x + y * y,
phi = 0,
i = 10, delta;
do {
var tanPhi = tan(phi),
secPhi = 1 / cos(phi),
j = k - 2 * y * phi + phi * phi;
phi -= delta = (tanPhi * j + 2 * (phi - y)) / (2 + j * secPhi * secPhi + 2 * (phi - y) * tanPhi);
} while (abs(delta) > epsilon && --i > 0);
var E = x * (tanPhi = tan(phi)),
A = tan(abs(y) < abs(phi + 1 / tanPhi) ? asin(E) * 0.5 : acos(E) * 0.5 + pi / 4) / sin(phi);
return [
sinPhi0 ? 2 * atan(sinPhi0 * A) / sinPhi0 : 2 * A,
phi
];
};
return forward;
}
var rectangularPolyconic = function() {
return parallel1(rectangularPolyconicRaw)
.scale(131.215);
};
var K = [
[0.9986, -0.062],
[1.0000, 0.0000],
[0.9986, 0.0620],
[0.9954, 0.1240],
[0.9900, 0.1860],
[0.9822, 0.2480],
[0.9730, 0.3100],
[0.9600, 0.3720],
[0.9427, 0.4340],
[0.9216, 0.4958],
[0.8962, 0.5571],
[0.8679, 0.6176],
[0.8350, 0.6769],
[0.7986, 0.7346],
[0.7597, 0.7903],
[0.7186, 0.8435],
[0.6732, 0.8936],
[0.6213, 0.9394],
[0.5722, 0.9761],
[0.5322, 1.0000]
];
K.forEach(function(d) {
d[1] *= 1.0144;
});
function robinsonRaw(lambda, phi) {
var i = min(18, abs(phi) * 36 / pi),
i0 = floor(i),
di = i - i0,
ax = (k = K[i0])[0],
ay = k[1],
bx = (k = K[++i0])[0],
by = k[1],
cx = (k = K[min(19, ++i0)])[0],
cy = k[1],
k;
return [
lambda * (bx + di * (cx - ax) / 2 + di * di * (cx - 2 * bx + ax) / 2),
(phi > 0 ? halfPi : -halfPi) * (by + di * (cy - ay) / 2 + di * di * (cy - 2 * by + ay) / 2)
];
}
robinsonRaw.invert = function(x, y) {
var yy = y / halfPi,
phi = yy * 90,
i = min(18, abs(phi / 5)),
i0 = max(0, floor(i));
do {
var ay = K[i0][1],
by = K[i0 + 1][1],
cy = K[min(19, i0 + 2)][1],
u = cy - ay,
v = cy - 2 * by + ay,
t = 2 * (abs(yy) - by) / u,
c = v / u,
di = t * (1 - c * t * (1 - 2 * c * t));
if (di >= 0 || i0 === 1) {
phi = (y >= 0 ? 5 : -5) * (di + i);
var j = 50, delta;
do {
i = min(18, abs(phi) / 5);
i0 = floor(i);
di = i - i0;
ay = K[i0][1];
by = K[i0 + 1][1];
cy = K[min(19, i0 + 2)][1];
phi -= (delta = (y >= 0 ? halfPi : -halfPi) * (by + di * (cy - ay) / 2 + di * di * (cy - 2 * by + ay) / 2) - y) * degrees;
} while (abs(delta) > epsilon2 && --j > 0);
break;
}
} while (--i0 >= 0);
var ax = K[i0][0],
bx = K[i0 + 1][0],
cx = K[min(19, i0 + 2)][0];
return [
x / (bx + di * (cx - ax) / 2 + di * di * (cx - 2 * bx + ax) / 2),
phi * radians
];
};
var robinson = function() {
return d3Geo.geoProjection(robinsonRaw)
.scale(152.63);
};
function satelliteVerticalRaw(P) {
function forward(lambda, phi) {
var cosPhi = cos(phi),
k = (P - 1) / (P - cosPhi * cos(lambda));
return [
k * cosPhi * sin(lambda),
k * sin(phi)
];
}
forward.invert = function(x, y) {
var rho2 = x * x + y * y,
rho = sqrt(rho2),
sinc = (P - sqrt(1 - rho2 * (P + 1) / (P - 1))) / ((P - 1) / rho + rho / (P - 1));
return [
atan2(x * sinc, rho * sqrt(1 - sinc * sinc)),
rho ? asin(y * sinc / rho) : 0
];
};
return forward;
}
function satelliteRaw(P, omega) {
var vertical = satelliteVerticalRaw(P);
if (!omega) return vertical;
var cosOmega = cos(omega),
sinOmega = sin(omega);
function forward(lambda, phi) {
var coordinates = vertical(lambda, phi),
y = coordinates[1],
A = y * sinOmega / (P - 1) + cosOmega;
return [
coordinates[0] * cosOmega / A,
y / A
];
}
forward.invert = function(x, y) {
var k = (P - 1) / (P - 1 - y * sinOmega);
return vertical.invert(k * x, k * y * cosOmega);
};
return forward;
}
var satellite = function() {
var distance = 2,
omega = 0,
m = d3Geo.geoProjectionMutator(satelliteRaw),
p = m(distance, omega);
// As a multiple of radius.
p.distance = function(_) {
if (!arguments.length) return distance;
return m(distance = +_, omega);
};
p.tilt = function(_) {
if (!arguments.length) return omega * degrees;
return m(distance, omega = _ * radians);
};
return p
.scale(432.147)
.clipAngle(acos(1 / distance) * degrees - 1e-6);
};
var epsilon$1 = 1e-4;
var epsilonInverse = 1e4;
var x0 = -180;
var x0e = x0 + epsilon$1;
var x1 = 180;
var x1e = x1 - epsilon$1;
var y0 = -90;
var y0e = y0 + epsilon$1;
var y1 = 90;
var y1e = y1 - epsilon$1;
function nonempty(coordinates) {
return coordinates.length > 0;
}
function quantize$1(x) {
return Math.floor(x * epsilonInverse) / epsilonInverse;
}
function normalizePoint(y) {
return y === y0 || y === y1 ? [0, y] : [x0, quantize$1(y)]; // pole or antimeridian?
}
function clampPoint(p) {
var x = p[0], y = p[1], clamped = false;
if (x <= x0e) x = x0, clamped = true;
else if (x >= x1e) x = x1, clamped = true;
if (y <= y0e) y = y0, clamped = true;
else if (y >= y1e) y = y1, clamped = true;
return clamped ? [x, y] : p;
}
function clampPoints(points) {
return points.map(clampPoint);
}
// For each ring, detect where it crosses the antimeridian or pole.
function extractFragments(rings, polygon, fragments) {
for (var j = 0, m = rings.length; j < m; ++j) {
var ring = rings[j].slice();
// By default, assume that this ring doesn’t need any stitching.
fragments.push({index: -1, polygon: polygon, ring: ring});
for (var i = 0, n = ring.length; i < n; ++i) {
var point = ring[i],
x = point[0],
y = point[1];
// If this is an antimeridian or polar point…
if (x <= x0e || x >= x1e || y <= y0e || y >= y1e) {
ring[i] = clampPoint(point);
// Advance through any antimeridian or polar points…
for (var k = i + 1; k < n; ++k) {
var pointk = ring[k],
xk = pointk[0],
yk = pointk[1];
if (xk > x0e && xk < x1e && yk > y0e && yk < y1e) break;
}
// If this was just a single antimeridian or polar point,
// we don’t need to cut this ring into a fragment;
// we can just leave it as-is.
if (k === i + 1) continue;
// Otherwise, if this is not the first point in the ring,
// cut the current fragment so that it ends at the current point.
// The current point is also normalized for later joining.
if (i) {
var fragmentBefore = {index: -1, polygon: polygon, ring: ring.slice(0, i + 1)};
fragmentBefore.ring[fragmentBefore.ring.length - 1] = normalizePoint(y);
fragments[fragments.length - 1] = fragmentBefore;
}
// If the ring started with an antimeridian fragment,
// we can ignore that fragment entirely.
else fragments.pop();
// If the remainder of the ring is an antimeridian fragment,
// move on to the next ring.
if (k >= n) break;
// Otherwise, add the remaining ring fragment and continue.
fragments.push({index: -1, polygon: polygon, ring: ring = ring.slice(k - 1)});
ring[0] = normalizePoint(ring[0][1]);
i = -1;
n = ring.length;
}
}
}
}
// Now stitch the fragments back together into rings.
function stitchFragments(fragments) {
var i, n = fragments.length;
// To connect the fragments start-to-end, create a simple index by end.
var fragmentByStart = {},
fragmentByEnd = {},
fragment,
start,
startFragment,
end,
endFragment;
// For each fragment…
for (i = 0; i < n; ++i) {
fragment = fragments[i];
start = fragment.ring[0];
end = fragment.ring[fragment.ring.length - 1];
// If this fragment is closed, add it as a standalone ring.
if (start[0] === end[0] && start[1] === end[1]) {
fragment.polygon.push(fragment.ring);
fragments[i] = null;
continue;
}
fragment.index = i;
fragmentByStart[start] = fragmentByEnd[end] = fragment;
}
// For each open fragment…
for (i = 0; i < n; ++i) {
fragment = fragments[i];
if (fragment) {
start = fragment.ring[0];
end = fragment.ring[fragment.ring.length - 1];
startFragment = fragmentByEnd[start];
endFragment = fragmentByStart[end];
delete fragmentByStart[start];
delete fragmentByEnd[end];
// If this fragment is closed, add it as a standalone ring.
if (start[0] === end[0] && start[1] === end[1]) {
fragment.polygon.push(fragment.ring);
continue;
}
if (startFragment) {
delete fragmentByEnd[start];
delete fragmentByStart[startFragment.ring[0]];
startFragment.ring.pop(); // drop the shared coordinate
fragments[startFragment.index] = null;
fragment = {index: -1, polygon: startFragment.polygon, ring: startFragment.ring.concat(fragment.ring)};
if (startFragment === endFragment) {
// Connect both ends to this single fragment to create a ring.
fragment.polygon.push(fragment.ring);
} else {
fragment.index = n++;
fragments.push(fragmentByStart[fragment.ring[0]] = fragmentByEnd[fragment.ring[fragment.ring.length - 1]] = fragment);
}
} else if (endFragment) {
delete fragmentByStart[end];
delete fragmentByEnd[endFragment.ring[endFragment.ring.length - 1]];
fragment.ring.pop(); // drop the shared coordinate
fragment = {index: n++, polygon: endFragment.polygon, ring: fragment.ring.concat(endFragment.ring)};
fragments[endFragment.index] = null;
fragments.push(fragmentByStart[fragment.ring[0]] = fragmentByEnd[fragment.ring[fragment.ring.length - 1]] = fragment);
} else {
fragment.ring.push(fragment.ring[0]); // close ring
fragment.polygon.push(fragment.ring);
}
}
}
}
function stitchFeature(input) {
var output = {type: "Feature", geometry: stitchGeometry(input.geometry)};
if (input.id != null) output.id = input.id;
if (input.bbox != null) output.bbox = input.bbox;
if (input.properties != null) output.properties = input.properties;
return output;
}
function stitchGeometry(input) {
if (input == null) return input;
var output, fragments, i, n;
switch (input.type) {
case "GeometryCollection": output = {type: "GeometryCollection", geometries: input.geometries.map(stitchGeometry)}; break;
case "Point": output = {type: "Point", coordinates: clampPoint(input.coordinates)}; break;
case "MultiPoint": case "LineString": output = {type: input.type, coordinates: clampPoints(input.coordinates)}; break;
case "MultiLineString": output = {type: "MultiLineString", coordinates: input.coordinates.map(clampPoints)}; break;
case "Polygon": {
var polygon = [];
extractFragments(input.coordinates, polygon, fragments = []);
stitchFragments(fragments);
output = {type: "Polygon", coordinates: polygon};
break;
}
case "MultiPolygon": {
fragments = [], i = -1, n = input.coordinates.length;
var polygons = new Array(n);
while (++i < n) extractFragments(input.coordinates[i], polygons[i] = [], fragments);
stitchFragments(fragments);
output = {type: "MultiPolygon", coordinates: polygons.filter(nonempty)};
break;
}
default: return input;
}
if (input.bbox != null) output.bbox = input.bbox;
return output;
}
var stitch = function(input) {
if (input == null) return input;
switch (input.type) {
case "Feature": return stitchFeature(input);
case "FeatureCollection": {
var output = {type: "FeatureCollection", features: input.features.map(stitchFeature)};
if (input.bbox != null) output.bbox = input.bbox;
return output;
}
default: return stitchGeometry(input);
}
};
function timesRaw(lambda, phi) {
var t = tan(phi / 2),
s = sin(quarterPi * t);
return [
lambda * (0.74482 - 0.34588 * s * s),
1.70711 * t
];
}
timesRaw.invert = function(x, y) {
var t = y / 1.70711,
s = sin(quarterPi * t);
return [
x / (0.74482 - 0.34588 * s * s),
2 * atan(t)
];
};
var times = function() {
return d3Geo.geoProjection(timesRaw)
.scale(146.153);
};
// Compute the origin as the midpoint of the two reference points.
// Rotate one of the reference points by the origin.
// Apply the spherical law of sines to compute gamma rotation.
var twoPoint = function(raw, p0, p1) {
var i = d3Geo.geoInterpolate(p0, p1),
o = i(0.5),
a = d3Geo.geoRotation([-o[0], -o[1]])(p0),
b = i.distance / 2,
y = -asin(sin(a[1] * radians) / sin(b)),
R = [-o[0], -o[1], -(a[0] > 0 ? pi - y : y) * degrees],
p = d3Geo.geoProjection(raw(b)).rotate(R),
r = d3Geo.geoRotation(R),
center = p.center;
delete p.rotate;
p.center = function(_) {
return arguments.length ? center(r(_)) : r.invert(center());
};
return p
.clipAngle(90);
};
function twoPointAzimuthalRaw(d) {
var cosd = cos(d);
function forward(lambda, phi) {
var coordinates = d3Geo.geoGnomonicRaw(lambda, phi);
coordinates[0] *= cosd;
return coordinates;
}
forward.invert = function(x, y) {
return d3Geo.geoGnomonicRaw.invert(x / cosd, y);
};
return forward;
}
function twoPointAzimuthalUsa() {
return twoPointAzimuthal([-158, 21.5], [-77, 39])
.clipAngle(60)
.scale(400);
}
function twoPointAzimuthal(p0, p1) {
return twoPoint(twoPointAzimuthalRaw, p0, p1);
}
// TODO clip to ellipse
function twoPointEquidistantRaw(z0) {
if (!(z0 *= 2)) return d3Geo.geoAzimuthalEquidistantRaw;
var lambdaa = -z0 / 2,
lambdab = -lambdaa,
z02 = z0 * z0,
tanLambda0 = tan(lambdab),
S = 0.5 / sin(lambdab);
function forward(lambda, phi) {
var za = acos(cos(phi) * cos(lambda - lambdaa)),
zb = acos(cos(phi) * cos(lambda - lambdab)),
ys = phi < 0 ? -1 : 1;
za *= za, zb *= zb;
return [
(za - zb) / (2 * z0),
ys * sqrt(4 * z02 * zb - (z02 - za + zb) * (z02 - za + zb)) / (2 * z0)
];
}
forward.invert = function(x, y) {
var y2 = y * y,
cosza = cos(sqrt(y2 + (t = x + lambdaa) * t)),
coszb = cos(sqrt(y2 + (t = x + lambdab) * t)),
t,
d;
return [
atan2(d = cosza - coszb, t = (cosza + coszb) * tanLambda0),
(y < 0 ? -1 : 1) * acos(sqrt(t * t + d * d) * S)
];
};
return forward;
}
function twoPointEquidistantUsa() {
return twoPointEquidistant([-158, 21.5], [-77, 39])
.clipAngle(130)
.scale(122.571);
}
function twoPointEquidistant(p0, p1) {
return twoPoint(twoPointEquidistantRaw, p0, p1);
}
function vanDerGrintenRaw(lambda, phi) {
if (abs(phi) < epsilon) return [lambda, 0];
var sinTheta = abs(phi / halfPi),
theta = asin(sinTheta);
if (abs(lambda) < epsilon || abs(abs(phi) - halfPi) < epsilon) return [0, sign(phi) * pi * tan(theta / 2)];
var cosTheta = cos(theta),
A = abs(pi / lambda - lambda / pi) / 2,
A2 = A * A,
G = cosTheta / (sinTheta + cosTheta - 1),
P = G * (2 / sinTheta - 1),
P2 = P * P,
P2_A2 = P2 + A2,
G_P2 = G - P2,
Q = A2 + G;
return [
sign(lambda) * pi * (A * G_P2 + sqrt(A2 * G_P2 * G_P2 - P2_A2 * (G * G - P2))) / P2_A2,
sign(phi) * pi * (P * Q - A * sqrt((A2 + 1) * P2_A2 - Q * Q)) / P2_A2
];
}
vanDerGrintenRaw.invert = function(x, y) {
if (abs(y) < epsilon) return [x, 0];
if (abs(x) < epsilon) return [0, halfPi * sin(2 * atan(y / pi))];
var x2 = (x /= pi) * x,
y2 = (y /= pi) * y,
x2_y2 = x2 + y2,
z = x2_y2 * x2_y2,
c1 = -abs(y) * (1 + x2_y2),
c2 = c1 - 2 * y2 + x2,
c3 = -2 * c1 + 1 + 2 * y2 + z,
d = y2 / c3 + (2 * c2 * c2 * c2 / (c3 * c3 * c3) - 9 * c1 * c2 / (c3 * c3)) / 27,
a1 = (c1 - c2 * c2 / (3 * c3)) / c3,
m1 = 2 * sqrt(-a1 / 3),
theta1 = acos(3 * d / (a1 * m1)) / 3;
return [
pi * (x2_y2 - 1 + sqrt(1 + 2 * (x2 - y2) + z)) / (2 * x),
sign(y) * pi * (-m1 * cos(theta1 + pi / 3) - c2 / (3 * c3))
];
};
var vanDerGrinten = function() {
return d3Geo.geoProjection(vanDerGrintenRaw)
.scale(79.4183);
};
function vanDerGrinten2Raw(lambda, phi) {
if (abs(phi) < epsilon) return [lambda, 0];
var sinTheta = abs(phi / halfPi),
theta = asin(sinTheta);
if (abs(lambda) < epsilon || abs(abs(phi) - halfPi) < epsilon) return [0, sign(phi) * pi * tan(theta / 2)];
var cosTheta = cos(theta),
A = abs(pi / lambda - lambda / pi) / 2,
A2 = A * A,
x1 = cosTheta * (sqrt(1 + A2) - A * cosTheta) / (1 + A2 * sinTheta * sinTheta);
return [
sign(lambda) * pi * x1,
sign(phi) * pi * sqrt(1 - x1 * (2 * A + x1))
];
}
vanDerGrinten2Raw.invert = function(x, y) {
if (!x) return [0, halfPi * sin(2 * atan(y / pi))];
var x1 = abs(x / pi),
A = (1 - x1 * x1 - (y /= pi) * y) / (2 * x1),
A2 = A * A,
B = sqrt(A2 + 1);
return [
sign(x) * pi * (B - A),
sign(y) * halfPi * sin(2 * atan2(sqrt((1 - 2 * A * x1) * (A + B) - x1), sqrt(B + A + x1)))
];
};
var vanDerGrinten2 = function() {
return d3Geo.geoProjection(vanDerGrinten2Raw)
.scale(79.4183);
};
function vanDerGrinten3Raw(lambda, phi) {
if (abs(phi) < epsilon) return [lambda, 0];
var sinTheta = phi / halfPi,
theta = asin(sinTheta);
if (abs(lambda) < epsilon || abs(abs(phi) - halfPi) < epsilon) return [0, pi * tan(theta / 2)];
var A = (pi / lambda - lambda / pi) / 2,
y1 = sinTheta / (1 + cos(theta));
return [
pi * (sign(lambda) * sqrt(A * A + 1 - y1 * y1) - A),
pi * y1
];
}
vanDerGrinten3Raw.invert = function(x, y) {
if (!y) return [x, 0];
var y1 = y / pi,
A = (pi * pi * (1 - y1 * y1) - x * x) / (2 * pi * x);
return [
x ? pi * (sign(x) * sqrt(A * A + 1) - A) : 0,
halfPi * sin(2 * atan(y1))
];
};
var vanDerGrinten3 = function() {
return d3Geo.geoProjection(vanDerGrinten3Raw)
.scale(79.4183);
};
function vanDerGrinten4Raw(lambda, phi) {
if (!phi) return [lambda, 0];
var phi0 = abs(phi);
if (!lambda || phi0 === halfPi) return [0, phi];
var B = phi0 / halfPi,
B2 = B * B,
C = (8 * B - B2 * (B2 + 2) - 5) / (2 * B2 * (B - 1)),
C2 = C * C,
BC = B * C,
B_C2 = B2 + C2 + 2 * BC,
B_3C = B + 3 * C,
lambda0 = lambda / halfPi,
lambda1 = lambda0 + 1 / lambda0,
D = sign(abs(lambda) - halfPi) * sqrt(lambda1 * lambda1 - 4),
D2 = D * D,
F = B_C2 * (B2 + C2 * D2 - 1) + (1 - B2) * (B2 * (B_3C * B_3C + 4 * C2) + 12 * BC * C2 + 4 * C2 * C2),
x1 = (D * (B_C2 + C2 - 1) + 2 * sqrt(F)) / (4 * B_C2 + D2);
return [
sign(lambda) * halfPi * x1,
sign(phi) * halfPi * sqrt(1 + D * abs(x1) - x1 * x1)
];
}
vanDerGrinten4Raw.invert = function(x, y) {
var delta;
if (!x || !y) return [x, y];
y /= pi;
var x1 = sign(x) * x / halfPi,
D = (x1 * x1 - 1 + 4 * y * y) / abs(x1),
D2 = D * D,
B = 2 * y,
i = 50;
do {
var B2 = B * B,
C = (8 * B - B2 * (B2 + 2) - 5) / (2 * B2 * (B - 1)),
C_ = (3 * B - B2 * B - 10) / (2 * B2 * B),
C2 = C * C,
BC = B * C,
B_C = B + C,
B_C2 = B_C * B_C,
B_3C = B + 3 * C,
F = B_C2 * (B2 + C2 * D2 - 1) + (1 - B2) * (B2 * (B_3C * B_3C + 4 * C2) + C2 * (12 * BC + 4 * C2)),
F_ = -2 * B_C * (4 * BC * C2 + (1 - 4 * B2 + 3 * B2 * B2) * (1 + C_) + C2 * (-6 + 14 * B2 - D2 + (-8 + 8 * B2 - 2 * D2) * C_) + BC * (-8 + 12 * B2 + (-10 + 10 * B2 - D2) * C_)),
sqrtF = sqrt(F),
f = D * (B_C2 + C2 - 1) + 2 * sqrtF - x1 * (4 * B_C2 + D2),
f_ = D * (2 * C * C_ + 2 * B_C * (1 + C_)) + F_ / sqrtF - 8 * B_C * (D * (-1 + C2 + B_C2) + 2 * sqrtF) * (1 + C_) / (D2 + 4 * B_C2);
B -= delta = f / f_;
} while (delta > epsilon && --i > 0);
return [
sign(x) * (sqrt(D * D + 4) + D) * pi / 4,
halfPi * B
];
};
var vanDerGrinten4 = function() {
return d3Geo.geoProjection(vanDerGrinten4Raw)
.scale(127.16);
};
var A = 4 * pi + 3 * sqrt(3);
var B = 2 * sqrt(2 * pi * sqrt(3) / A);
var wagner4Raw = mollweideBromleyRaw(B * sqrt(3) / pi, B, A / 6);
var wagner4 = function() {
return d3Geo.geoProjection(wagner4Raw)
.scale(176.84);
};
function wagner6Raw(lambda, phi) {
return [lambda * sqrt(1 - 3 * phi * phi / (pi * pi)), phi];
}
wagner6Raw.invert = function(x, y) {
return [x / sqrt(1 - 3 * y * y / (pi * pi)), y];
};
var wagner6 = function() {
return d3Geo.geoProjection(wagner6Raw)
.scale(152.63);
};
function wagner7Raw(lambda, phi) {
var s = 0.90631 * sin(phi),
c0 = sqrt(1 - s * s),
c1 = sqrt(2 / (1 + c0 * cos(lambda /= 3)));
return [
2.66723 * c0 * c1 * sin(lambda),
1.24104 * s * c1
];
}
wagner7Raw.invert = function(x, y) {
var t1 = x / 2.66723,
t2 = y / 1.24104,
p = sqrt(t1 * t1 + t2 * t2),
c = 2 * asin(p / 2);
return [
3 * atan2(x * tan(c), 2.66723 * p),
p && asin(y * sin(c) / (1.24104 * 0.90631 * p))
];
};
var wagner7 = function() {
return d3Geo.geoProjection(wagner7Raw)
.scale(172.632);
};
function wiechelRaw(lambda, phi) {
var cosPhi = cos(phi),
sinPhi = cos(lambda) * cosPhi,
sin1_Phi = 1 - sinPhi,
cosLambda = cos(lambda = atan2(sin(lambda) * cosPhi, -sin(phi))),
sinLambda = sin(lambda);
cosPhi = sqrt(1 - sinPhi * sinPhi);
return [
sinLambda * cosPhi - cosLambda * sin1_Phi,
-cosLambda * cosPhi - sinLambda * sin1_Phi
];
}
wiechelRaw.invert = function(x, y) {
var w = (x * x + y * y) / -2,
k = sqrt(-w * (2 + w)),
b = y * w + x * k,
a = x * w - y * k,
D = sqrt(a * a + b * b);
return [
atan2(k * b, D * (1 + w)),
D ? -asin(k * a / D) : 0
];
};
var wiechel = function() {
return d3Geo.geoProjection(wiechelRaw)
.rotate([0, -90, 45])
.scale(124.75)
.clipAngle(180 - 1e-3);
};
function winkel3Raw(lambda, phi) {
var coordinates = aitoffRaw(lambda, phi);
return [
(coordinates[0] + lambda / halfPi) / 2,
(coordinates[1] + phi) / 2
];
}
winkel3Raw.invert = function(x, y) {
var lambda = x, phi = y, i = 25;
do {
var cosphi = cos(phi),
sinphi = sin(phi),
sin_2phi = sin(2 * phi),
sin2phi = sinphi * sinphi,
cos2phi = cosphi * cosphi,
sinlambda = sin(lambda),
coslambda_2 = cos(lambda / 2),
sinlambda_2 = sin(lambda / 2),
sin2lambda_2 = sinlambda_2 * sinlambda_2,
C = 1 - cos2phi * coslambda_2 * coslambda_2,
E = C ? acos(cosphi * coslambda_2) * sqrt(F = 1 / C) : F = 0,
F,
fx = 0.5 * (2 * E * cosphi * sinlambda_2 + lambda / halfPi) - x,
fy = 0.5 * (E * sinphi + phi) - y,
dxdlambda = 0.5 * F * (cos2phi * sin2lambda_2 + E * cosphi * coslambda_2 * sin2phi) + 0.5 / halfPi,
dxdphi = F * (sinlambda * sin_2phi / 4 - E * sinphi * sinlambda_2),
dydlambda = 0.125 * F * (sin_2phi * sinlambda_2 - E * sinphi * cos2phi * sinlambda),
dydphi = 0.5 * F * (sin2phi * coslambda_2 + E * sin2lambda_2 * cosphi) + 0.5,
denominator = dxdphi * dydlambda - dydphi * dxdlambda,
dlambda = (fy * dxdphi - fx * dydphi) / denominator,
dphi = (fx * dydlambda - fy * dxdlambda) / denominator;
lambda -= dlambda, phi -= dphi;
} while ((abs(dlambda) > epsilon || abs(dphi) > epsilon) && --i > 0);
return [lambda, phi];
};
var winkel3 = function() {
return d3Geo.geoProjection(winkel3Raw)
.scale(158.837);
};
exports.geoAiry = airy;
exports.geoAiryRaw = airyRaw;
exports.geoAitoff = aitoff;
exports.geoAitoffRaw = aitoffRaw;
exports.geoArmadillo = armadillo;
exports.geoArmadilloRaw = armadilloRaw;
exports.geoAugust = august;
exports.geoAugustRaw = augustRaw;
exports.geoBaker = baker;
exports.geoBakerRaw = bakerRaw;
exports.geoBerghaus = berghaus;
exports.geoBerghausRaw = berghausRaw;
exports.geoBoggs = boggs;
exports.geoBoggsRaw = boggsRaw;
exports.geoBonne = bonne;
exports.geoBonneRaw = bonneRaw;
exports.geoBottomley = bottomley;
exports.geoBottomleyRaw = bottomleyRaw;
exports.geoBromley = bromley;
exports.geoBromleyRaw = bromleyRaw;
exports.geoChamberlin = chamberlin;
exports.geoChamberlinRaw = chamberlinRaw;
exports.geoChamberlinAfrica = chamberlinAfrica;
exports.geoCollignon = collignon;
exports.geoCollignonRaw = collignonRaw;
exports.geoCraig = craig;
exports.geoCraigRaw = craigRaw;
exports.geoCraster = craster;
exports.geoCrasterRaw = crasterRaw;
exports.geoCylindricalEqualArea = cylindricalEqualArea;
exports.geoCylindricalEqualAreaRaw = cylindricalEqualAreaRaw;
exports.geoCylindricalStereographic = cylindricalStereographic;
exports.geoCylindricalStereographicRaw = cylindricalStereographicRaw;
exports.geoEckert1 = eckert1;
exports.geoEckert1Raw = eckert1Raw;
exports.geoEckert2 = eckert2;
exports.geoEckert2Raw = eckert2Raw;
exports.geoEckert3 = eckert3;
exports.geoEckert3Raw = eckert3Raw;
exports.geoEckert4 = eckert4;
exports.geoEckert4Raw = eckert4Raw;
exports.geoEckert5 = eckert5;
exports.geoEckert5Raw = eckert5Raw;
exports.geoEckert6 = eckert6;
exports.geoEckert6Raw = eckert6Raw;
exports.geoEisenlohr = eisenlohr;
exports.geoEisenlohrRaw = eisenlohrRaw;
exports.geoFahey = fahey;
exports.geoFaheyRaw = faheyRaw;
exports.geoFoucaut = foucaut;
exports.geoFoucautRaw = foucautRaw;
exports.geoGilbert = gilbert;
exports.geoGingery = gingery;
exports.geoGingeryRaw = gingeryRaw;
exports.geoGinzburg4 = ginzburg4;
exports.geoGinzburg4Raw = ginzburg4Raw;
exports.geoGinzburg5 = ginzburg5;
exports.geoGinzburg5Raw = ginzburg5Raw;
exports.geoGinzburg6 = ginzburg6;
exports.geoGinzburg6Raw = ginzburg6Raw;
exports.geoGinzburg8 = ginzburg8;
exports.geoGinzburg8Raw = ginzburg8Raw;
exports.geoGinzburg9 = ginzburg9;
exports.geoGinzburg9Raw = ginzburg9Raw;
exports.geoGringorten = gringorten;
exports.geoGringortenRaw = gringortenRaw;
exports.geoGuyou = guyou;
exports.geoGuyouRaw = guyouRaw;
exports.geoHammer = hammer;
exports.geoHammerRaw = hammerRaw;
exports.geoHammerRetroazimuthal = hammerRetroazimuthal;
exports.geoHammerRetroazimuthalRaw = hammerRetroazimuthalRaw;
exports.geoHealpix = healpix;
exports.geoHealpixRaw = healpixRaw;
exports.geoHill = hill;
exports.geoHillRaw = hillRaw;
exports.geoHomolosine = homolosine;
exports.geoHomolosineRaw = homolosineRaw;
exports.geoInterrupt = interrupt;
exports.geoInterruptedBoggs = boggs$1;
exports.geoInterruptedHomolosine = homolosine$1;
exports.geoInterruptedMollweide = mollweide$1;
exports.geoInterruptedMollweideHemispheres = mollweideHemispheres;
exports.geoInterruptedSinuMollweide = sinuMollweide$1;
exports.geoInterruptedSinusoidal = sinusoidal$1;
exports.geoKavrayskiy7 = kavrayskiy7;
exports.geoKavrayskiy7Raw = kavrayskiy7Raw;
exports.geoLagrange = lagrange;
exports.geoLagrangeRaw = lagrangeRaw;
exports.geoLarrivee = larrivee;
exports.geoLarriveeRaw = larriveeRaw;
exports.geoLaskowski = laskowski;
exports.geoLaskowskiRaw = laskowskiRaw;
exports.geoLittrow = littrow;
exports.geoLittrowRaw = littrowRaw;
exports.geoLoximuthal = loximuthal;
exports.geoLoximuthalRaw = loximuthalRaw;
exports.geoMiller = miller;
exports.geoMillerRaw = millerRaw;
exports.geoModifiedStereographic = modifiedStereographic;
exports.geoModifiedStereographicRaw = modifiedStereographicRaw;
exports.geoModifiedStereographicAlaska = modifiedStereographicAlaska;
exports.geoModifiedStereographicGs48 = modifiedStereographicGs48;
exports.geoModifiedStereographicGs50 = modifiedStereographicGs50;
exports.geoModifiedStereographicMiller = modifiedStereographicMiller;
exports.geoModifiedStereographicLee = modifiedStereographicLee;
exports.geoMollweide = mollweide;
exports.geoMollweideRaw = mollweideRaw;
exports.geoMtFlatPolarParabolic = mtFlatPolarParabolic;
exports.geoMtFlatPolarParabolicRaw = mtFlatPolarParabolicRaw;
exports.geoMtFlatPolarQuartic = mtFlatPolarQuartic;
exports.geoMtFlatPolarQuarticRaw = mtFlatPolarQuarticRaw;
exports.geoMtFlatPolarSinusoidal = mtFlatPolarSinusoidal;
exports.geoMtFlatPolarSinusoidalRaw = mtFlatPolarSinusoidalRaw;
exports.geoNaturalEarth = naturalEarth;
exports.geoNaturalEarthRaw = naturalEarthRaw;
exports.geoNaturalEarth2 = naturalEarth2;
exports.geoNaturalEarth2Raw = naturalEarth2Raw;
exports.geoNellHammer = nellHammer;
exports.geoNellHammerRaw = nellHammerRaw;
exports.geoPatterson = patterson;
exports.geoPattersonRaw = pattersonRaw;
exports.geoPolyconic = polyconic;
exports.geoPolyconicRaw = polyconicRaw;
exports.geoPolyhedral = polyhedral;
exports.geoPolyhedralButterfly = butterfly;
exports.geoPolyhedralCollignon = collignon$1;
exports.geoPolyhedralWaterman = waterman;
exports.geoProject = index;
exports.geoGringortenQuincuncial = gringorten$1;
exports.geoPeirceQuincuncial = peirce;
exports.geoPierceQuincuncial = peirce;
exports.geoQuantize = quantize;
exports.geoQuincuncial = quincuncial;
exports.geoRectangularPolyconic = rectangularPolyconic;
exports.geoRectangularPolyconicRaw = rectangularPolyconicRaw;
exports.geoRobinson = robinson;
exports.geoRobinsonRaw = robinsonRaw;
exports.geoSatellite = satellite;
exports.geoSatelliteRaw = satelliteRaw;
exports.geoSinuMollweide = sinuMollweide;
exports.geoSinuMollweideRaw = sinuMollweideRaw;
exports.geoSinusoidal = sinusoidal;
exports.geoSinusoidalRaw = sinusoidalRaw;
exports.geoStitch = stitch;
exports.geoTimes = times;
exports.geoTimesRaw = timesRaw;
exports.geoTwoPointAzimuthal = twoPointAzimuthal;
exports.geoTwoPointAzimuthalRaw = twoPointAzimuthalRaw;
exports.geoTwoPointAzimuthalUsa = twoPointAzimuthalUsa;
exports.geoTwoPointEquidistant = twoPointEquidistant;
exports.geoTwoPointEquidistantRaw = twoPointEquidistantRaw;
exports.geoTwoPointEquidistantUsa = twoPointEquidistantUsa;
exports.geoVanDerGrinten = vanDerGrinten;
exports.geoVanDerGrintenRaw = vanDerGrintenRaw;
exports.geoVanDerGrinten2 = vanDerGrinten2;
exports.geoVanDerGrinten2Raw = vanDerGrinten2Raw;
exports.geoVanDerGrinten3 = vanDerGrinten3;
exports.geoVanDerGrinten3Raw = vanDerGrinten3Raw;
exports.geoVanDerGrinten4 = vanDerGrinten4;
exports.geoVanDerGrinten4Raw = vanDerGrinten4Raw;
exports.geoWagner4 = wagner4;
exports.geoWagner4Raw = wagner4Raw;
exports.geoWagner6 = wagner6;
exports.geoWagner6Raw = wagner6Raw;
exports.geoWagner7 = wagner7;
exports.geoWagner7Raw = wagner7Raw;
exports.geoWiechel = wiechel;
exports.geoWiechelRaw = wiechelRaw;
exports.geoWinkel3 = winkel3;
exports.geoWinkel3Raw = winkel3Raw;
Object.defineProperty(exports, '__esModule', { value: true });
})));
Raw
d3-geo-voronoi.js
// https://github.com/Fil/d3-geo-voronoi Version 0.0.3. Copyright 2017 Philippe Riviere.
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-array'), require('d3-collection'), require('d3-geo'), require('d3-voronoi')) :
typeof define === 'function' && define.amd ? define(['exports', 'd3-array', 'd3-collection', 'd3-geo', 'd3-voronoi'], factory) :
(factory((global.d3 = global.d3 || {}),global.d3,global.d3,global.d3,global.d3));
}(this, (function (exports,d3Array,d3Collection,d3Geo,d3Voronoi) { 'use strict';
//
// Copyright (c) 2016, by Loren Petrich
//
// Distributed under the terms of the MIT License
//
// http://lpetrich.org/
//
// Calculates the spherical Delaunay triangulation of a set of points
// These points are entered as an array of arrays of coordinates: 0, 1, 2
// Any extra members are ignored
// FindDelaunayTriangulation(Positions) and
// FindDelaunayTriangulationIndexed(Positions, Indices)
// work from an array of points as specified above,
// the second one working from a set of indices into the array,
// and return an object with these members:
// positions -- vectors on a unit sphere
// indices -- of all the vertices
// triangles -- array of TriangleObject
// edges -- array of EdgeObject
// hull -- array of vertex indices -- the convex hull
// vor_positions -- positions of triangles' circumcircle centers (Voronoi vertices)
// vor_edges -- pair of indices in vor_positions (empty one: [-1,-1])
// vor_polygons -- object indexed by vertex index,
// and containing edges (EdgeObject), triangles (TriangleObject),
// and boundary (vertices in vor_positions)
// Open ones have a -1 at each end.
// Warning: ImproveTriangulation() is mysteriously buggy
// and is effectively disabled for now
function dotprd(x,y)
{
var sum = 0.0;
for (var ic=0; ic<3; ic++)
sum += x[ic]*y[ic];
return sum;
}
function crsprd(x,y)
{
var prod = new Array(3);
for (var ic=0; ic<3; ic++)
{
var ic1 = ic + 1;
if (ic1 >= 3) ic1 -= 3;
var ic2 = ic + 2;
if (ic2 >= 3) ic2 -= 3;
prod[ic] = x[ic1]*y[ic2] - x[ic2]*y[ic1];
}
return prod;
}
function triple_prd(x,y,z)
{
return dotprd(crsprd(x,y),z);
}
// This distance formula has some nice properties:
// distance and not square of distance;
// the square roots give better numerical resolution
// distance of antipode(p) to p' = - (distance of p to p')
// Range: -2 to +2
function ptdist(x,y)
{
var dst1 = 0.0;
var dst2 = 0.0;
for (var ic=0; ic<3; ic++)
{
var diff1 = y[ic] - x[ic];
dst1 += diff1*diff1;
var diff2 = y[ic] + x[ic];
dst2 += diff2*diff2;
}
return Math.sqrt(dst1) - Math.sqrt(dst2);
}
function Normalize(vec)
{
var vecres = new Array(3);
var sum = 0.0, nrmult;
for (var ic=0; ic<3; ic++)
{
var val = vec[ic];
sum += val*val;
}
if (sum > 0)
nrmult = 1/Math.sqrt(sum);
else
nrmult = 0;
for (var ic=0; ic<3; ic++)
{
vecres[ic] = nrmult*vec[ic];
}
return vecres;
}
function crsprd_ix(Positions,ix,iy)
{
return crsprd(Positions[ix],Positions[iy]);
}
function triple_prd_ix(Positions,ix,iy,iz)
{
return triple_prd(Positions[ix],Positions[iy],Positions[iz]);
}
function ptdist_ix(Positions,ix,iy)
{
return ptdist(Positions[ix],Positions[iy]);
}
// Returns a zero 3-vector
function zerovec()
{
var vec = new Array(3);
for (var ic=0; ic<3; ic++)
vec[ic] = 0.0;
return vec;
}
// Implements copying
function vec_copy(x)
{
var vec = new Array(3);
for (var ic=0; ic<3; ic++)
vec[ic] = x[ic];
return vec;
}
// Implements x += y
function vec_add_to(x,y)
{
for (var ic=0; ic<3; ic++)
x[ic] += y[ic];
}
// Implements x *= y
function vec_mult_scalar_to(x,y)
{
for (var ic=0; ic<3; ic++)
x[ic] *= y;
}
// Implements x - y
function vec_difference(x,y)
{
var diff = zerovec();
for (var ic=0; ic<3; ic++)
diff[ic] = x[ic] - y[ic];
return diff;
}
// JavaScript's counterpart of "null" / "None":
function IsNull(x)
{
return (typeof(x) == 'undefined');
}
function TrianglesEqual(tr1, tr2)
{
if (IsNull(tr1)) return false;
if (IsNull(tr2)) return false;
for (var iv=0; iv<3; iv++)
if (tr1.verts[iv] != tr2.verts[iv])
return false;
return true;
}
function EdgesEqual(ed1, ed2)
{
if (IsNull(ed1)) return false;
if (IsNull(ed2)) return false;
for (var iv=0; iv<2; iv++)
if (ed1.verts[iv] != ed2.verts[iv])
return false;
return true;
}
function max(x,y)
{
return (y > x) ? y : x;
}
function TriangleObject(Positions, verts)
{
this.verts = verts;
this.edges = new Array(3);
// Find directions for testing whether a point is inside
this.dirs = new Array(3);
for (var ic=0; ic<3; ic++)
{
var ic1 = ic + 1;
if (ic1 >= 3) ic1 -= 3;
var ic2 = ic + 2;
if (ic2 >= 3) ic2 -= 3;
this.dirs[ic] = crsprd_ix(Positions,verts[ic1],verts[ic2]);
}
// Tetrahedral volume factor
this.vol = triple_prd_ix(Positions,verts[0],verts[1],verts[2]);
// Adjust to get the signs correct for the point-inside test;
// the vertex opposite the edges' vertices ought to give a dot product of 1
for (var ic=0; ic<3; ic++)
vec_mult_scalar_to(this.dirs[ic],1/this.vol);
// Circumcircle test
var ccdir = zerovec();
for (var ic=0; ic<3; ic++)
vec_add_to(ccdir,this.dirs[ic]);
this.ccdir = Normalize(ccdir);
var ccdsq = 0;
for (var ic=0; ic<3; ic++)
ccdsq += ptdist(this.ccdir,Positions[verts[ic]]);
ccdsq /= 3;
this.ccdsq = ccdsq;
}
// For copying in vertex info from another triangle
TriangleObject.prototype.copy_vert_info = function(src)
{
this.verts = src.verts;
this.dirs = src.dirs;
this.vol = src.vol;
this.ccdir = src.ccdir;
this.ccdsq = src.ccdsq;
};
TriangleObject.prototype.IsVertexOrderCorrect = function()
{
return this.vol >= 0;
};
TriangleObject.prototype.IsPointInside = function(p)
{
for (var ic=0; ic<3; ic++)
if (dotprd(p,this.dirs[ic]) < 0) return false;
return true;
};
TriangleObject.prototype.IsPointInCircumcircle = function(p)
{
return (ptdist(this.ccdir,p) < this.ccdsq);
};
TriangleObject.prototype.IsVertex = function(ix)
{
for (var ic=0; ic<3; ic++)
if (ix == this.verts[ic]) return true;
return false;
};
TriangleObject.prototype.VertexIndexIn = function(ix)
{
for (var ic=0; ic<3; ic++)
if (ix == this.verts[ic]) return ic;
return -1;
};
TriangleObject.prototype.EdgeIndexIn = function(ed)
{
for (var ic=0; ic<3; ic++)
if (EdgesEqual(this.edges[ic], ed)) return ic;
return -1;
};
function EdgeObject(verts)
{
this.verts = verts;
this.polys = new Array(2);
}
EdgeObject.prototype.IsVertex = function(ix)
{
for (var ic=0; ic<2; ic++)
if (ix == this.verts[ic]) return true;
return false;
};
EdgeObject.prototype.VertexIndexIn = function(ix)
{
for (var ic=0; ic<2; ic++)
if (ix == this.verts[ic]) return ic;
return -1;
};
EdgeObject.prototype.PolyIndexIn = function(pl)
{
for (var ic=0; ic<2; ic++)
if (TrianglesEqual(this.polys[ic],pl)) return ic;
return -1;
};
function EdgeCheckObject(Positions,verts)
{
this.verts = verts;
this.pdst = ptdist_ix(Positions,verts[0],verts[1]);
this.direc = Normalize(crsprd_ix(Positions,verts[0],verts[1]));
var midpnt = zerovec();
vec_add_to(midpnt,Positions[verts[0]]);
vec_add_to(midpnt,Positions[verts[1]]);
this.midpnt = Normalize(midpnt);
}
// Check on the possible intersection with another edge-check object
// return a boolean of whether or not it does
EdgeCheckObject.prototype.intersects = function(Positions,other)
{
// Assume that sharing a vertex means non-intersecting
for (var ic=0; ic<2; ic++)
for (var ict=0; ict<2; ict++)
if (this.verts[ic] == other.verts[ict]) return false;
// Find intersection point; will test it and its antipode
var itsc = Normalize(crsprd(this.direc, other.direc));
// Find dot product with midpoints to test if the intersection
// is in the near hemispheres of the lines' midpoints.
// If it is in both near hemispheres or both far hemispheres, it's OK
// In both far hemispheres: antipode is in both near hemispheres
var near0 = dotprd(itsc,this.midpnt) > 0;
var near1 = dotprd(itsc,other.midpnt) > 0;
if (near1 != near0) return false;
var pd0 = [];
for (var ic=0; ic<2; ic++)
{
var pd = ptdist(itsc, Positions[this.verts[ic]]);
pd0.push(pd);
}
var pd1 = [];
for (var ic=0; ic<2; ic++)
{
var pd = ptdist(itsc, Positions[other.verts[ic]]);
pd1.push(pd);
}
var mxpd0 = max(pd0[0],pd0[1]);
var mxpd1 = max(pd1[0],pd1[1]);
if ((mxpd0 <= this.pdst) && (mxpd1 <= other.pdst) && near0) return true;
// Handle its antipode; use antipode-related shortcuts
// like reversing the distance value and the hemisphere-presence value
vec_mult_scalar_to(itsc, -1);
near0 = !near0;
for (var ic=0; ic<2; ic++)
{
pd0[ic] = - pd0[ic];
pd1[ic] = - pd1[ic];
}
mxpd0 = max(pd0[0],pd0[1]);
mxpd1 = max(pd1[0],pd1[1]);
if ((mxpd0 <= this.pdst) && (mxpd1 <= other.pdst) && near0) return true;
return false;
};
// Adds to an array if it was not already present;
// Must resort to this kludge because JavaScript doesn't have sets
function AddUnique(arr, x)
{
for (var i=0; i<arr.length; i++)
if (x == arr[i]) return;
arr.push(x);
}
// Version for edges, since testing equality of objects
// doesn't work that well in JavaScript
function AddUniqueEdge(arr, ed)
{
for (var i=0; i<arr.length; i++)
if (EdgesEqual(arr[i],ed)) return;
arr.push(ed);
}
// Find the set intersection
function FindShared(arr1, arr2)
{
var resarr = [];
for (var i1=0; i1<arr1.length; i1++)
{
var x1 = arr1[i1];
for (var i2=0; i2<arr2.length; i2++)
{
var x2 = arr2[i2];
if (x1 == x2)
{
resarr.push(x1);
break;
}
}
}
return resarr;
}
// Version for edges
function FindSharedEdges(arr1, arr2)
{
var resarr = [];
for (var i1=0; i1<arr1.length; i1++)
{
var x1 = arr1[i1];
for (var i2=0; i2<arr2.length; i2++)
{
var x2 = arr2[i2];
if (EdgesEqual(x1,x2))
{
resarr.push(x1);
break;
}
}
}
return resarr;
}
// Takes all the members of of arr2 out of arr1
// and ignores the arr2 members not present in arr1
function FindSetDifference(arr1, arr2)
{
if (arr2.length == 0) return;
var diffarr = [];
for (var i1=0; i1<arr1.length; i1++)
{
var x1 = arr1[i1];
var AddThisOne = true;
for (var i2=0; i2<arr2.length; i2++)
{
var x2 = arr2[i2];
if (x2 == x1)
{
AddThisOne = false;
break;
}
}
if (AddThisOne) diffarr.push(x1);
}
// Clear the array
arr1.splice(0,arr1.length);
for (var i=0; i<diffarr.length; i++)
arr1.push(diffarr[i]);
}
// Version for edges
function FindSetDifferenceEdges(arr1, arr2)
{
if (arr2.length == 0) return;
var diffarr = [];
for (var i1=0; i1<arr1.length; i1++)
{
var x1 = arr1[i1];
var AddThisOne = true;
for (var i2=0; i2<arr2.length; i2++)
{
var x2 = arr2[i2];
if (EdgesEqual(x1,x2))
{
AddThisOne = false;
break;
}
}
if (AddThisOne) diffarr.push(x1);
}
// Clear the array
arr1.splice(0,arr1.length);
for (var i=0; i<diffarr.length; i++)
arr1.push(diffarr[i]);
}
// Specified by index ix; returns whether it was possible to do so
function AddPointInside(TriSet, ix)
{
var Positions = TriSet.positions;
var p = Positions[ix];
var NumTris = TriSet.triangles.length;
for (var j=0; j<NumTris; j++)
{
var tri = TriSet.triangles[j];
if (tri.IsPointInside(p))
{
// Create three new triangles and their edges
var eds = tri.edges;
var trixs = [];
for (var ic=0; ic<3; ic++)
trixs.push(eds[ic].PolyIndexIn(tri));
var newtris = Array(3);
var neweds = Array(3);
for (var ic=0; ic<3; ic++)
{
var ic1 = ic + 1;
if (ic1 >= 3) ic1 -= 3;
newtris[ic] = new TriangleObject(Positions,[tri.verts[ic],tri.verts[ic1],ix]);
neweds[ic] = new EdgeObject([tri.verts[ic],ix]);
}
// Connect those triangles and edges
for (var ic=0; ic<3; ic++)
{
var ic1 = ic + 1;
if (ic1 >= 3) ic1 -= 3;
newtris[ic].edges[0] = neweds[ic1];
newtris[ic].edges[1] = neweds[ic];
neweds[ic].polys[0] = newtris[ic];
neweds[ic1].polys[1] = newtris[ic];
}
// Find which external edges go with which triangles
for (var ic=0; ic<3; ic++)
{
var ed = eds[ic];
var trix = trixs[ic];
for (var ict=0; ict<3; ict++)
{
var newtri = newtris[ict];
var numverts = 0;
for (var iv=0; iv<2; iv++)
{
if (newtri.IsVertex(ed.verts[iv]))
numverts++;
if (numverts == 2)
{
ed.polys[trix] = newtri;
newtri.edges[2] = ed;
break;
}
}
}
}
// Insert those triangles and edges into the lists
TriSet.triangles[j] = newtris[0];
for (var ic=1; ic<3; ic++)
TriSet.triangles.push(newtris[ic]);
for (var ic=0; ic<3; ic++)
TriSet.edges.push(neweds[ic]);
// All done; indicate that the point was added
return true;
}
}
// The point was inside no triangle, and thus was not added
return false;
}
function ImproveTriangulation(TriSet)
{
var Positions = TriSet.positions;
var quad_verts = new Array(4);
for (var itr=0; itr<100; itr++)
{
var numflips = 0;
for (var i=0; i<TriSet.edges.length; i++)
{
var edge = TriSet.edges[i];
var tris = edge.polys;
// Skip over external edges
if (IsNull(tris[0])) continue;
if (IsNull(tris[1])) continue;
// Find the containing quadrangle's vertices
for (var ic=0; ic<3; ic++)
{
var ix = tris[0].verts[ic];
if (!edge.IsVertex(ix)) break;
}
var ic1 = ic + 1;
if (ic1 >= 3) ic1 -= 3;
var ic2 = ic + 2;
if (ic2 >= 3) ic2 -= 3;
quad_verts[0] = ix;
quad_verts[1] = tris[0].verts[ic1];
quad_verts[3] = tris[0].verts[ic2];
for (var ic=0; ic<3; ic++)
{
var ix = tris[1].verts[ic];
if (!edge.IsVertex(ix)) break;
}
quad_verts[2] = ix;
// Are the non-edge points in the other triangles' circumcircles?
var incc0 = tris[0].IsPointInCircumcircle(Positions[quad_verts[2]]);
var incc1 = tris[1].IsPointInCircumcircle(Positions[quad_verts[0]]);
if ((!incc0) && (!incc1)) continue;
// Are the would-be triangles properly oriented?
var newtri0 = new TriangleObject(Positions, [quad_verts[0],quad_verts[1],quad_verts[2]]);
if (!newtri0.IsVertexOrderCorrect()) continue;
var newtri1 = new TriangleObject(Positions, [quad_verts[0],quad_verts[2],quad_verts[3]]);
if (!newtri1.IsVertexOrderCorrect()) continue;
// If so, then flip
numflips++;
// Adjust the edge and triangle memberships:
// 0-3 goes from 0 to 1, 1-2 goes from 1 to 0
for (var ic=0; ic<3; ic++)
{
var ed = tris[0].edges[ic];
if (EdgesEqual(ed,edge)) continue;
else if (ed.IsVertex(quad_verts[3]))
{
var ed03 = ed;
var edix03 = ic;
break;
}
}
for (var ic=0; ic<3; ic++)
{
var ed = tris[1].edges[ic];
if (EdgesEqual(ed,edge)) continue;
else if (ed.IsVertex(quad_verts[1]))
{
var ed12 = ed;
var edix12 = ic;
break;
}
}
var trix0 = ed03.PolyIndexIn(tris[0]);
var trix1 = ed12.PolyIndexIn(tris[1]);
ed03.polys[trix0] = tris[1];
ed12.polys[trix1] = tris[0];
tris[0].edges[edix03] = ed12;
tris[1].edges[edix12] = ed03;
// Add the vertices
tris[0].copy_vert_info(newtri0);
tris[1].copy_vert_info(newtri1);
edge.verts = [quad_verts[0], quad_verts[2]];
}
if (numflips == 0) break;
}
}
function FindConvexHull(TriSet)
{
// var Positions = TriSet.positions;
// Find boundary loop -- use as convex hull
var NextVertex = new Object;
var VtxStart = -1;
for (var i=0; i<TriSet.edges.length; i++)
{
var edge = TriSet.edges[i];
// Find a boundary one -- look for the triangle that it contains
if (IsNull(edge.polys[0]))
{
if (IsNull(edge.polys[1]))
continue;
else
var tri = edge.polys[1];
}
else
{
if (IsNull(edge.polys[1]))
var tri = edge.polys[0];
else
continue;
}
// Ensure that the hull is in the same direction as the triangles
var ix0 = edge.verts[0];
var ix1 = edge.verts[1];
var posdiff = tri.VertexIndexIn(ix1) - tri.VertexIndexIn(ix0);
if (posdiff < 0) posdiff += 3;
if (posdiff != 1)
{
var ixs = ix0;
ix0 = ix1;
ix1 = ixs;
}
NextVertex[ix0] = ix1;
VtxStart = ix0;
}
if (VtxStart >= 0)
{
var ix = VtxStart;
var hull = [ix];
while(true)
{
var ixnext = NextVertex[ix];
if (ixnext == VtxStart) break;
hull.push(ixnext);
ix = ixnext;
}
TriSet.hull = hull;
}
}
// Finds the dual of the Delaunay triangulation
// Won't bother to do the sort of connectivity
// that was necessary for the Delaunay triangulation
function FindVoronoiDiagram(TriSet)
{
// Special cases: 3 or fewer points
if (TriSet.triangles.length == 1)
{
// A single triangle
if (TriSet.hull.length == 3)
{
var tri = TriSet.triangles[0];
TriSet.vor_positions.push(tri.ccdir);
for (var k=0; k<3; k++)
{
var kx = k + 1;
if (kx >= 3) kx = 0;
var ky = k - 1;
if (ky < 0) ky = 2;
var v1 = TriSet.positions[TriSet.hull[k]];
var v2 = TriSet.positions[TriSet.hull[kx]];
var posdiff = vec_difference(v2,v1);
TriSet.vor_positions.push(Normalize(crsprd(posdiff,tri.ccdir)));
TriSet.vor_edges.push([0, k+1, 4]);
var ix = TriSet.hull[k];
TriSet.vor_polygons[ix] = new Object;
var vor_poly = TriSet.vor_polygons[ix];
var iy = TriSet.hull[ky];
for (var l=0; l<3; l++)
{
var edge = TriSet.edges[l];
var shrd = FindShared([iy,ix],edge.verts);
if (shrd.length == 2) break;
}
vor_poly.edges = [edge];
vor_poly.triangles = [tri];
vor_poly.boundary = [0, ky+1, 4, k+1];
}
var ept = vec_copy(tri.ccdir);
vec_mult_scalar_to(ept,-1);
TriSet.vor_positions.push(ept);
}
return;
}
else if (TriSet.triangles.length == 0)
{
// A biangle
if (TriSet.hull.length == 2)
{
var v0 = TriSet.positions[TriSet.hull[0]];
var v1 = TriSet.positions[TriSet.hull[1]];
var vt0 = zerovec();
vec_add_to(vt0,v0);
vec_add_to(vt0,v1);
vt0 = Normalize(vt0);
TriSet.vor_positions.push(vt0);
var vt1 = Normalize(crsprd(v0,v1));
TriSet.vor_positions.push(vt1);
var vt2 = vec_copy(vt0);
vec_mult_scalar_to(vt2,-1);
TriSet.vor_positions.push(vt2);
var vt3 = vec_copy(vt1);
vec_mult_scalar_to(vt3,-1);
TriSet.vor_positions.push(vt3);
TriSet.vor_edges.push([0, 1, 2, 3, 0]);
edge = TriSet.edges[0];
for (var k=0; k<2; k++)
{
var ix = TriSet.hull[k];
TriSet.vor_polygons[ix] = new Object;
var vor_poly = TriSet.vor_polygons[ix];
vor_poly.edges = [edge];
vor_poly.triangles = [0];
if (k == 0)
vor_poly.boundary = [0, 1, 2, 3];
else if (k == 1)
vor_poly.boundary = [0, 3, 2, 1];
}
}
return;
}
// Create the array of Voronoi-vertex positions:
// Add indices to the triangle objects for convenience
for (var i=0; i<TriSet.triangles.length; i++)
{
var tri = TriSet.triangles[i];
tri.index = i;
TriSet.vor_positions.push(tri.ccdir);
}
// Voronoi edges: a cinch
// Voronoi edges parallel original edges
for (var i=0; i<TriSet.edges.length; i++)
{
var edge = TriSet.edges[i];
if (!IsNull(edge.polys[0]) && !IsNull(edge.polys[1]))
{
var vor_edge = [edge.polys[0].index, edge.polys[1].index];
TriSet.vor_edges.push(vor_edge);
}
}
// Voronoi polygons: -1 at ends means an open one
// First, collect the edges and triangles at each vertex
// Put them into vor_polygons, because each one
// is for each original vertex
for (var i=0; i<TriSet.indices.length; i++)
{
var ix = TriSet.indices[i];
TriSet.vor_polygons[ix] = new Object;
var vor_poly = TriSet.vor_polygons[ix];
vor_poly.edges = [];
vor_poly.triangles = [];
vor_poly.boundary = [];
}
for (var i=0; i<TriSet.edges.length; i++)
{
var edge = TriSet.edges[i];
for (var ic=0; ic<2; ic++)
TriSet.vor_polygons[edge.verts[ic]].edges.push(edge);
}
for (var i=0; i<TriSet.triangles.length; i++)
{
var tri = TriSet.triangles[i];
for (var ic=0; ic<3; ic++)
TriSet.vor_polygons[tri.verts[ic]].triangles.push(tri);
}
for (var i=0; i<TriSet.indices.length; i++)
{
var ix = TriSet.indices[i];
var vor_poly = TriSet.vor_polygons[ix];
// First triangle
var init_tri = vor_poly.triangles[0];
var tri = init_tri;
vor_poly.boundary.push(tri.index);
// First edge
for (var ic=0; ic<3; ic++)
{
var edge = tri.edges[ic];
if (edge.IsVertex(ix))
break;
}
var init_edge = edge;
// The next triangle and edge
var IsInside = false;
while(true)
{
var iv = edge.PolyIndexIn(tri);
tri = edge.polys[1-iv];
if (IsNull(tri)) break;
if (TrianglesEqual(tri,init_tri))
{
IsInside = true;
break;
}
vor_poly.boundary.push(tri.index);
for (var ic=0; ic<3; ic++)
{
var next_edge = tri.edges[ic];
if (EdgesEqual(next_edge,edge)) continue;
if (next_edge.IsVertex(ix))
{
edge = next_edge;
break;
}
}
}
if (!IsInside)
{
vor_poly.boundary.reverse();
tri = init_tri;
// First edge the other way
for (var ic=0; ic<3; ic++)
{
edge = tri.edges[ic];
if (EdgesEqual(edge,init_edge)) continue;
if (edge.IsVertex(ix))
break;
}
while(true)
{
var iv = edge.PolyIndexIn(tri);
tri = edge.polys[1-iv];
if (IsNull(tri)) break;
vor_poly.boundary.push(tri.index);
for (var ic=0; ic<3; ic++)
{
var next_edge = tri.edges[ic];
if (EdgesEqual(next_edge,edge)) continue;
if (next_edge.IsVertex(ix))
{
edge = next_edge;
break;
}
}
}
}
// Add -1 on ends for open polygon:
if (!IsInside)
{
vor_poly.boundary.reverse();
vor_poly.boundary.push(-1);
vor_poly.boundary.reverse();
vor_poly.boundary.push(-1);
}
}
// Handle the area outside of the convex hull
if (TriSet.hull.length >= 3)
{
// Set up the initial boundary lines
// The boundary lines contain:
// Index of Voronoi vertex / triangle center / intersection (in VorPos)
// Indices of original vertices on each side of the line
var VorBdLns = new Array();
var Positions = TriSet.positions;
var hlen = TriSet.hull.length;
for (var ic=0; ic<hlen; ic++)
{
var ix = TriSet.hull[ic];
var icx = ic + 1;
if (icx >= hlen) icx = 0;
var ixa = TriSet.hull[icx];
var edset1 = TriSet.vor_polygons[ix].edges;
var edset2 = TriSet.vor_polygons[ixa].edges;
var edsetshr = FindSharedEdges(edset1,edset2);
var edge = edsetshr[0];
var tvrt = edge.polys[0].index;
var vt0 = Positions[ix];
var vt1 = Positions[ixa];
var vtdf = vec_difference(vt1,vt0);
// Contains: triangle index (Voronoi vertex),
// vertex index 1 (Voronoi region), position
// vertex index 2 (Voronoi region), position,
// great-circle normal
var VorBdLn = [tvrt, TriSet.vor_positions[tvrt], ix, vt0, ixa, vt1, vtdf];
VorBdLns.push(VorBdLn);
}
// Find intersections
while (VorBdLns.length > 3)
{
// Check all combinations of neighbors
var n = VorBdLns.length;
var itscpts = new Array();
var ptitscs = new Array();
for (var k=0; k<n; k++)
ptitscs.push(new Array());
for (var k=0; k<n; k++)
{
// Find the intersection point; use the convex hull's direction
var kx = k + 1;
if (kx >= n) kx = 0;
var itscpt = Normalize(crsprd(VorBdLns[k][6],VorBdLns[kx][6]));
vec_mult_scalar_to(itscpt,-1);
ptitscs[k].push(itscpts.length);
ptitscs[kx].push(itscpts.length);
itscpts.push(itscpt);
}
// Find the intersection points that are the closest to their parent points
for (var k=0; k<n; k++)
{
var ptitsc = ptitscs[k];
if (ptitsc.length >= 2)
{
var dists = new Array();
for (var kp=0; kp<ptitsc.length; kp++)
dists.push(ptdist(itscpts[ptitsc[kp]],VorBdLns[k][1]));
var dx = 0;
var dmin = dists[dx];
for (var dxi=0; dxi<dists.length; dxi++)
{
var dst = dists[dxi];
if (dst < dmin)
{
dx = dxi; dmin = dst;
}
}
var ptitscrd = ptitsc[dx];
}
else if (ptitsc.length == 1)
var ptitscrd = ptitsc[0];
else
var ptitscrd = -1;
ptitscs[k] = ptitscrd;
}
var NewVorBdLns = new Array();
for (var k=0; k<n; k++)
{
// Find all matched intersection points and add them
var kx = k + 1;
if (kx >= n) kx = 0;
var ky = k - 1;
if (ky < 0) ky = n - 1;
// 0 is lone, 1 is leading, 2 is trailing
// vorvtidx is the index of the Voronoi vertex
var pstat = 0;
var ptitsc = ptitscs[k], ptitsc_next;
if (ptitsc != -1)
{
var ptitsc_prev = ptitscs[ky];
if (ptitsc == ptitsc_prev)
pstat = 2;
else
{
ptitsc_next = ptitscs[kx];
if (ptitsc == ptitsc_next)
pstat = 1;
}
}
if (pstat == 0)
{
// Keep the Voronoi line without merging
NewVorBdLns.push(VorBdLns[k]);
}
else if (pstat == 1)
{
// Merge the Voronoi lines and create a new one
var VorBdLn0 = VorBdLns[k];
var VorBdLn1 = VorBdLns[kx];
var itscpt = itscpts[k];
var tvrt0 = VorBdLn0[0];
var tvrt1 = VorBdLn1[0];
var PointOK = (tvrt1 != tvrt0);
if (PointOK)
{
var nitx = TriSet.vor_positions.length;
var ix0 = VorBdLn0[2];
var vt0 = VorBdLn0[3];
var ix1 = VorBdLn1[4];
var vt1 = VorBdLn1[5];
var dst_in = undefined;
var dst_out = undefined;
for (var m=0; m<n; m++)
{
var ptstm = ptdist(VorBdLns[m][3],itscpt);
var mrl = m - k;
while (mrl < 0) mrl += n;
while (mrl >= n) mrl -= n;
if (mrl <= 2)
{
if (dst_in == undefined)
dst_in = ptstm;
else if (ptstm < dst_in)
dst_in = ptstm;
}
else
{
if (dst_out == undefined)
dst_out = ptstm;
else if (ptstm < dst_out)
dst_out = ptstm;
}
}
PointOK = (dst_in < dst_out);
}
if (PointOK)
{
var vtdf = vec_difference(vt1,vt0);
var VorBdLn = [nitx, itscpt, ix0, vt0, ix1, vt1, vtdf];
NewVorBdLns.push(VorBdLn);
TriSet.vor_positions.push(itscpt);
var ixi = VorBdLn0[4];
// Should be equal:
// ixi = VorBdLn2[2];
TriSet.vor_edges.push([tvrt0, nitx]);
TriSet.vor_edges.push([tvrt1, nitx]);
// Add the point to the center Voronoi region and close it
TriSet.vor_polygons[ixi].boundary.shift();
var vpln = TriSet.vor_polygons[ixi].boundary.length;
TriSet.vor_polygons[ixi].boundary[vpln-1] = nitx;
// Add the point to the left Voronoi region
if (TriSet.vor_polygons[ix0].boundary[1] == tvrt0)
{
TriSet.vor_polygons[ix0].boundary.unshift(-1);
TriSet.vor_polygons[ix0].boundary[1] = nitx;
}
else
{
vpln = TriSet.vor_polygons[ix0].boundary.length;
if (TriSet.vor_polygons[ix0].boundary[vpln-2] == tvrt0)
{
TriSet.vor_polygons[ix0].boundary.push(-1);
vpln = TriSet.vor_polygons[ix0].boundary.length;
TriSet.vor_polygons[ix0].boundary[vpln-2] = nitx;
}
}
// Add the point to the right Voronoi region
if (TriSet.vor_polygons[ix1].boundary[1] == tvrt1)
{
TriSet.vor_polygons[ix1].boundary.unshift(-1);
TriSet.vor_polygons[ix1].boundary[1] = nitx;
}
else
{
vpln = TriSet.vor_polygons[ix1].boundary.length;
if (TriSet.vor_polygons[ix1].boundary[vpln-2] == tvrt1)
{
TriSet.vor_polygons[ix1].boundary.push(-1);
vpln = TriSet.vor_polygons[ix1].boundary.length;
TriSet.vor_polygons[ix1].boundary[vpln-2] = nitx;
}
}
}
else
{
NewVorBdLns.push(VorBdLn0);
NewVorBdLns.push(VorBdLn1);
}
}
/*
else if (pstat == 2)
{
// Do nothing
}
*/
}
if (NewVorBdLns.length == VorBdLns.length) break;
VorBdLns = NewVorBdLns;
}
// Special cases: only two or three points left
if (VorBdLns.length == 2)
{
if (VorBdLns[0][0] != VorBdLns[1][0])
{
var VorLn = [];
for (var k=0; k<2; k++)
{
// Connecting line
var kx = VorBdLns[k][0];
VorLn.push(kx);
// Close the Voronoi region by deleting the end -1's
kx = VorBdLns[k][2];
TriSet.vor_polygons[kx].boundary.shift();
TriSet.vor_polygons[kx].boundary.pop();
}
TriSet.vor_edges.push(VorLn);
}
}
else if (VorBdLns.length == 3)
{
var ic0 = VorBdLns[0][0];
var ic1 = VorBdLns[1][0];
var ic2 = VorBdLns[2][0];
if (ic0 != ic1 && ic0 != ic2 && ic1 != ic2)
{
var nitx = TriSet.vor_positions.length;
var v0 = VorBdLns[0][3];
var v1 = VorBdLns[1][3];
var v2 = VorBdLns[2][3];
var itscpt = zerovec();
vec_add_to(itscpt,crsprd(v0,v1));
vec_add_to(itscpt,crsprd(v1,v2));
vec_add_to(itscpt,crsprd(v2,v0));
itscpt = Normalize(itscpt);
vec_mult_scalar_to(itscpt,-1);
TriSet.vor_positions.push(itscpt);
for (var k=0; k<3; k++)
{
// Connecting line
var VorBdLn = VorBdLns[k];
TriSet.vor_edges.push([VorBdLn[0], nitx]);
// Add the point to the Voronoi region and close it
var ix = VorBdLn[2];
TriSet.vor_polygons[ix].boundary.shift();
var vpln = TriSet.vor_polygons[ix].boundary.length;
TriSet.vor_polygons[ix].boundary[vpln-1] = nitx;
}
}
}
}
// Adjust the orientations
for (var k=0; k<TriSet.vor_polygons.length; k++)
{
vor_poly = TriSet.vor_polygons[k];
if (vor_poly.boundary.length >= 3 && vor_poly.boundary[0] >= 0)
{
tri = new TriangleObject(TriSet.vor_positions,vor_poly.boundary.slice(0,3));
if (!tri.IsVertexOrderCorrect())
vor_poly.boundary.reverse();
}
}
}
function FindDelaunayTriangulationIndexed(Positions, Indices)
{
// Create the triangle-set object
var TriSet = new Object;
TriSet.positions = Positions;
TriSet.indices = Indices;
TriSet.triangles = [];
TriSet.edges = [];
TriSet.hull = [];
TriSet.vor_positions = [];
TriSet.vor_edges = [];
TriSet.vor_polygons = new Object;
// Create the first triangle, if it is possible to create any
if (Indices.length < 3)
{
if (Indices.length == 2)
{
TriSet.edges.push(new EdgeObject(Indices));
TriSet.hull = Indices;
}
FindVoronoiDiagram(TriSet);
return TriSet;
}
var tri = new TriangleObject(Positions, Indices.slice(0,3));
if (!tri.IsVertexOrderCorrect())
tri = new TriangleObject(Positions, [Indices[0],Indices[2],Indices[1]]);
TriSet.triangles.push(tri);
var echs = new Array(3);
for (var ic=0; ic<3; ic++)
{
var ic1 = ic + 1;
if (ic1 >= 3) ic1 -= 3;
var ix = Indices[ic];
var ix1 = Indices[ic1];
var vts = [ix, ix1];
var edge = new EdgeObject(vts);
var echeck = new EdgeCheckObject(Positions,vts);
echeck.edge = edge;
echs[ic] = echeck;
tri.edges[ic] = edge;
edge.polys[0] = tri;
TriSet.edges.push(edge);
}
// Place those crossing checkers in a boundary object;
// will have to use various kludges since JavaScript doesn't have sets
var BoundaryVerts = Indices.slice(0,3);
var BoundaryEdges = echs;
var Verts = Object;
for (var ic=0; ic<3; ic++)
{
var ic1 = ic + 2;
if (ic1 >= 3) ic1 -= 3;
var ix = Indices[ic];
Verts[ix] = [echs[ic],echs[ic+1]];
}
// Add points until it is no longer possible
for (var i=3; i<Indices.length; i++)
{
var ix = Indices[i];
// If possible, add the point inside
if (AddPointInside(TriSet, ix)) continue;
// Point was not inside
Verts[ix] = [];
var NewEdges = [];
var VertsAddedTo = [];
var EdgesToDelete = [];
// Find all the non-intersecting edges
for (var j=0; j<BoundaryVerts.length; j++)
{
var ix1 = BoundaryVerts[j];
var echk = new EdgeCheckObject(Positions, [ix, ix1]);
var DoesIntersect = false;
for (var k=0; k<BoundaryEdges.length; k++)
{
var echk1 = BoundaryEdges[k];
DoesIntersect = echk.intersects(Positions,echk1);
if (DoesIntersect) break;
}
if (DoesIntersect) continue;
var edge = new EdgeObject(echk.verts);
echk.edge = edge;
AddUniqueEdge(NewEdges,echk);
AddUniqueEdge(Verts[ix],echk);
AddUnique(VertsAddedTo,ix);
AddUniqueEdge(Verts[ix1],echk);
AddUnique(VertsAddedTo,ix1);
}
// Add the new vertex itself
AddUnique(BoundaryVerts,ix);
// Find all the triangles
for (var j=0; j<BoundaryEdges.length; j++)
{
var echk = BoundaryEdges[j];
var echks = [];
for (var iv=0; iv<2; iv++)
{
var vset = FindSharedEdges(Verts[ix],Verts[echk.verts[iv]]);
if (vset.length == 0) continue;
echks.push(vset[0]);
}
if (echks.length < 2) continue;
var empt_indx = -1;
for (var iv=0; iv<2; iv++)
{
if (IsNull(echk.edge.polys[iv]))
{
empt_indx = iv;
break;
}
}
if (empt_indx < 0) continue;
var oldtri = echk.edge.polys[1-empt_indx];
var v0 = echk.verts[0];
var i0 = oldtri.VertexIndexIn(v0);
var v1 = echk.verts[1];
var i1 = oldtri.VertexIndexIn(v1);
var i01 = i1 - i0;
if (i01 < 0) i01 += 3;
if (i01 == 1)
{
// Order in original: other, v0, v1
var NewTriVerts = [ix, v1, v0];
}
else if (i01 == 2)
{
// Order in original: other, v1, v0
var NewTriVerts = [ix, v0, v1];
}
var tri = new TriangleObject(Positions, NewTriVerts);
if (!tri.IsVertexOrderCorrect()) continue;
// Add the new triangle
// Also, add the new edges,
// or remove them from the lists if necessary
TriSet.triangles.push(tri);
echk.edge.polys[empt_indx] = tri;
tri.edges[0] = echk.edge;
tri.edges[1] = echks[0].edge;
tri.edges[2] = echks[1].edge;
AddUniqueEdge(EdgesToDelete, echk);
for (var iv=0; iv<2; iv++)
{
var echki = echks[iv];
if (IsNull(echki.edge.polys[0]))
{
echki.edge.polys[0] = tri;
TriSet.edges.push(echki.edge);
}
else
{
echki.edge.polys[1] = tri;
AddUniqueEdge(EdgesToDelete,echki);
}
}
}
// Add the new edges and remove the edges and vertices
// that are now in the interior
for (var j=0; j<NewEdges.length; j++)
AddUniqueEdge(BoundaryEdges,NewEdges[j]);
FindSetDifferenceEdges(BoundaryEdges, EdgesToDelete);
var BoundaryVertsToRemove = [];
for (var j=0; j<VertsAddedTo.length; j++)
{
var ixa = VertsAddedTo[j];
FindSetDifferenceEdges(Verts[ixa],EdgesToDelete);
if (Verts[ixa].length == 0)
BoundaryVertsToRemove.push(ixa);
}
FindSetDifference(BoundaryVerts, BoundaryVertsToRemove);
}
// Improve it iteratively
ImproveTriangulation(TriSet);
// Find the boundary line of this region
FindConvexHull(TriSet);
// Find the regions around each point:
FindVoronoiDiagram(TriSet);
return TriSet;
}
function FindDelaunayTriangulation(Positions)
{
var Indices = new Array(Positions.length);
for (var i=0; i<Indices.length; i++)
Indices[i] = i;
return FindDelaunayTriangulationIndexed(Positions, Indices);
}
//
// (c) 2016 Philippe Riviere
//
// https://github.com/Fil/
//
// This software is distributed under the terms of the MIT License
var geoVoronoi = function () {
var radians = Math.PI / 180;
var cartesian = function (spherical) {
var lambda = spherical[0] * radians,
phi = spherical[1] * radians,
cosphi = Math.cos(phi);
return [
cosphi * Math.cos(lambda),
cosphi * Math.sin(lambda),
Math.sin(phi)
];
};
var spherical = function (cartesian) {
var r = Math.sqrt(cartesian[0] * cartesian[0] + cartesian[1] * cartesian[1]),
lat = Math.atan2(cartesian[2], r),
lng = Math.atan2(cartesian[1], cartesian[0]);
return [lng / radians, lat / radians];
};
var mapline = function (Positions, Verts) {
return Verts
.map(function (v) {
return spherical(Positions[v]);
});
};
var diagram = d3Voronoi.voronoi()([]);
var DT = diagram.DT = null,
sites = diagram.sites = [],
pos = diagram.pos = [],
x = function (d) {
if (typeof d == 'object' && 'type' in d) {
return d3Geo.geoCentroid(d)[0];
}
if (0 in d) return d[0];
},
y = function (d) {
if (typeof d == 'object' && 'type' in d) {
return d3Geo.geoCentroid(d)[1];
}
if (0 in d) return d[1];
};
var voro = function (data) {
diagram._hull = diagram._polygons = diagram._links = diagram._triangles = null;
if (typeof data == 'object' && data.type == 'FeatureCollection') {
data = data.features;
}
sites = data.map(function (site, i) {
site.index = i;
return site;
});
pos = data.map(function (site) {
return [x(site), y(site)];
});
DT = FindDelaunayTriangulation(pos.map(cartesian));
return diagram;
};
diagram.links = voro.links = function (s) {
if (s) voro(s);
if (diagram._links) return diagram._links;
var _index = d3Collection.map();
var features = DT.edges.map(function (i, n) {
_index.set(d3Array.extent(i.verts), n);
var properties = {
source: sites[i.verts[0]],
target: sites[i.verts[1]],
urquhart: true, // will be changed to false later
length: d3Geo.geoDistance(pos[i.verts[0]], pos[i.verts[1]])
};
// add left and right sites (?)
// make GeoJSON
return {
type: "Feature",
geometry: {
type: 'LineString',
coordinates: [spherical(DT.positions[i.verts[0]]), spherical(DT.positions[i.verts[1]])]
},
properties: properties
};
});
// Urquhart Graph? tag longer link from each triangle
DT.triangles.forEach(function (t) {
var l = 0,
length = 0,
remove, v;
for (var j = 0; j < 3; j++) {
v = d3Array.extent([t.verts[j], t.verts[(j + 1) % 3]]);
var n = _index.get(v);
length = features[n].properties.length;
if (length > l) {
l = length;
remove = n;
}
}
features[remove].properties.urquhart = false;
});
return diagram._links = {
type: "FeatureCollection",
features: features
};
};
diagram.triangles = voro.triangles = function (s) {
if (s) voro(s);
if (diagram._triangles) return diagram._triangles;
var features = DT.triangles
.map(function (t) {
t.spherical = t.verts.map(function (v) {
return DT.positions[v];
})
.map(spherical);
// correct winding order
if (t.ccdsq < 0) {
t.spherical = t.spherical.reverse();
t.ccdsq *= -1;
}
return t;
})
// make geojson
.map(function (t) {
return {
type: "Feature",
geometry: {
type: "Polygon",
coordinates: [t.spherical.concat([t.spherical[0]])]
},
properties: {
sites: t.verts.map(function (i) {
return sites[i];
}),
area: t.vol, // steradians
circumcenter: spherical(t.ccdir),
// ccdsq is *not* the geodesic distance
/* circumradius: (2-t.ccdsq) * 53 */
}
}
});
return diagram._triangles = {
type: "FeatureCollection",
features: features
};
};
diagram.polygons = voro.polygons = function (s) {
if (s) voro(s);
if (diagram._polygons) return diagram._polygons;
var features = DT.indices.map(function (i, n) {
var vor_poly = DT.vor_polygons[DT.indices[i]];
var geometry = {};
if (vor_poly == undefined) {
geometry.type = "Sphere";
} else {
var line = mapline(DT.vor_positions,
vor_poly.boundary.concat([vor_poly.boundary[0]])
);
// correct winding order
var b = {
type: "Polygon",
coordinates: [[pos[i], line[0], line[1], pos[i]]]
};
if (d3Geo.geoArea(b) > 2 * Math.PI + 1e-10) {
line = line.reverse();
}
geometry.type = "Polygon";
geometry.coordinates = [line];
}
return {
type: "Feature",
geometry: geometry,
properties: {
site: sites[i],
sitecoordinates: pos[i],
neighbours: vor_poly.edges.map(function (e) {
return e.verts.filter(function (j) {
return j !== i;
})[0];
})
}
};
});
return diagram._polygons = {
type: "FeatureCollection",
features: features
};
};
diagram.hull = voro.hull = function (s) {
if (s) voro(s);
if (diagram._hull) return diagram._hull;
if (!DT.hull.length) {
return null; // What is a null GeoJSON?
}
// seems that DT.hull is always clockwise
var hull = DT.hull.reverse();
// make GeoJSON
return diagram._hull = {
type: "Feature",
geometry: {
type: "Polygon",
coordinates: [hull.concat([hull[0]]).map(function (i) {
return pos[i];
})]
},
properties: {
sites: hull.map(function (i) {
return sites[i];
})
}
};
};
diagram.find = function (x, y, radius) {
var features = diagram.polygons().features;
// optimization: start from most recent result
var i, next = diagram.find.found || 0;
var cell = features[next] || features[next = 0];
var dist = d3Geo.geoDistance([x, y], cell.properties.sitecoordinates);
do {
cell = features[i = next];
next = null;
cell.properties.neighbours.forEach(function (e) {
var ndist = d3Geo.geoDistance([x, y], features[e].properties.sitecoordinates);
if (ndist < dist) {
dist = ndist;
next = e;
return;
}
});
} while (next !== null);
diagram.find.found = i;
if (!radius || dist < radius * radius) return cell.properties.site;
};
voro.x = function (f) {
if (!f) return x;
x = f;
return voro;
};
voro.y = function (f) {
if (!f) return y;
y = f;
return voro;
};
voro.extent = function (f) {
if (!f) return null;
return voro;
};
voro.size = function (f) {
if (!f) return null;
return voro;
};
return voro;
};
exports.geoVoronoi = geoVoronoi;
Object.defineProperty(exports, '__esModule', { value: true });
})));
Raw
d3-geo.js
This file has been truncated, but you can view the full file.
// https://d3js.org/d3-geo/ Version 1.6.4. Copyright 2017 Mike Bostock.
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-array')) :
typeof define === 'function' && define.amd ? define(['exports', 'd3-array'], factory) :
(factory((global.d3 = global.d3 || {}),global.d3));
}(this, function (exports,d3Array) { 'use strict';
// Adds floating point numbers with twice the normal precision.
// Reference: J. R. Shewchuk, Adaptive Precision Floating-Point Arithmetic and
// Fast Robust Geometric Predicates, Discrete & Computational Geometry 18(3)
// 305–363 (1997).
// Code adapted from GeographicLib by Charles F. F. Karney,
// http://geographiclib.sourceforge.net/
function adder() {
return new Adder;
}
function Adder() {
this.reset();
}
Adder.prototype = {
constructor: Adder,
reset: function() {
this.s = // rounded value
this.t = 0; // exact error
},
add: function(y) {
add(temp, y, this.t);
add(this, temp.s, this.s);
if (this.s) this.t += temp.t;
else this.s = temp.t;
},
valueOf: function() {
return this.s;
}
};
var temp = new Adder;
function add(adder, a, b) {
var x = adder.s = a + b,
bv = x - a,
av = x - bv;
adder.t = (a - av) + (b - bv);
}
var epsilon = 1e-6;
var epsilon2 = 1e-12;
var pi = Math.PI;
var halfPi = pi / 2;
var quarterPi = pi / 4;
var tau = pi * 2;
var degrees = 180 / pi;
var radians = pi / 180;
var abs = Math.abs;
var atan = Math.atan;
var atan2 = Math.atan2;
var cos = Math.cos;
var ceil = Math.ceil;
var exp = Math.exp;
var log = Math.log;
var pow = Math.pow;
var sin = Math.sin;
var sign = Math.sign || function(x) { return x > 0 ? 1 : x < 0 ? -1 : 0; };
var sqrt = Math.sqrt;
var tan = Math.tan;
function acos(x) {
return x > 1 ? 0 : x < -1 ? pi : Math.acos(x);
}
function asin(x) {
return x > 1 ? halfPi : x < -1 ? -halfPi : Math.asin(x);
}
function haversin(x) {
return (x = sin(x / 2)) * x;
}
function noop() {}
function streamGeometry(geometry, stream) {
if (geometry && streamGeometryType.hasOwnProperty(geometry.type)) {
streamGeometryType[geometry.type](geometry, stream);
}
}
var streamObjectType = {
Feature: function(object, stream) {
streamGeometry(object.geometry, stream);
},
FeatureCollection: function(object, stream) {
var features = object.features, i = -1, n = features.length;
while (++i < n) streamGeometry(features[i].geometry, stream);
}
};
var streamGeometryType = {
Sphere: function(object, stream) {
stream.sphere();
},
Point: function(object, stream) {
object = object.coordinates;
stream.point(object[0], object[1], object[2]);
},
MultiPoint: function(object, stream) {
var coordinates = object.coordinates, i = -1, n = coordinates.length;
while (++i < n) object = coordinates[i], stream.point(object[0], object[1], object[2]);
},
LineString: function(object, stream) {
streamLine(object.coordinates, stream, 0);
},
MultiLineString: function(object, stream) {
var coordinates = object.coordinates, i = -1, n = coordinates.length;
while (++i < n) streamLine(coordinates[i], stream, 0);
},
Polygon: function(object, stream) {
streamPolygon(object.coordinates, stream);
},
MultiPolygon: function(object, stream) {
var coordinates = object.coordinates, i = -1, n = coordinates.length;
while (++i < n) streamPolygon(coordinates[i], stream);
},
GeometryCollection: function(object, stream) {
var geometries = object.geometries, i = -1, n = geometries.length;
while (++i < n) streamGeometry(geometries[i], stream);
}
};
function streamLine(coordinates, stream, closed) {
var i = -1, n = coordinates.length - closed, coordinate;
stream.lineStart();
while (++i < n) coordinate = coordinates[i], stream.point(coordinate[0], coordinate[1], coordinate[2]);
stream.lineEnd();
}
function streamPolygon(coordinates, stream) {
var i = -1, n = coordinates.length;
stream.polygonStart();
while (++i < n) streamLine(coordinates[i], stream, 1);
stream.polygonEnd();
}
function geoStream(object, stream) {
if (object && streamObjectType.hasOwnProperty(object.type)) {
streamObjectType[object.type](object, stream);
} else {
streamGeometry(object, stream);
}
}
var areaRingSum = adder();
var areaSum = adder();
var lambda00;
var phi00;
var lambda0;
var cosPhi0;
var sinPhi0;
var areaStream = {
point: noop,
lineStart: noop,
lineEnd: noop,
polygonStart: function() {
areaRingSum.reset();
areaStream.lineStart = areaRingStart;
areaStream.lineEnd = areaRingEnd;
},
polygonEnd: function() {
var areaRing = +areaRingSum;
areaSum.add(areaRing < 0 ? tau + areaRing : areaRing);
this.lineStart = this.lineEnd = this.point = noop;
},
sphere: function() {
areaSum.add(tau);
}
};
function areaRingStart() {
areaStream.point = areaPointFirst;
}
function areaRingEnd() {
areaPoint(lambda00, phi00);
}
function areaPointFirst(lambda, phi) {
areaStream.point = areaPoint;
lambda00 = lambda, phi00 = phi;
lambda *= radians, phi *= radians;
lambda0 = lambda, cosPhi0 = cos(phi = phi / 2 + quarterPi), sinPhi0 = sin(phi);
}
function areaPoint(lambda, phi) {
lambda *= radians, phi *= radians;
phi = phi / 2 + quarterPi; // half the angular distance from south pole
// Spherical excess E for a spherical triangle with vertices: south pole,
// previous point, current point. Uses a formula derived from Cagnoli’s
// theorem. See Todhunter, Spherical Trig. (1871), Sec. 103, Eq. (2).
var dLambda = lambda - lambda0,
sdLambda = dLambda >= 0 ? 1 : -1,
adLambda = sdLambda * dLambda,
cosPhi = cos(phi),
sinPhi = sin(phi),
k = sinPhi0 * sinPhi,
u = cosPhi0 * cosPhi + k * cos(adLambda),
v = k * sdLambda * sin(adLambda);
areaRingSum.add(atan2(v, u));
// Advance the previous points.
lambda0 = lambda, cosPhi0 = cosPhi, sinPhi0 = sinPhi;
}
function area(object) {
areaSum.reset();
geoStream(object, areaStream);
return areaSum * 2;
}
function spherical(cartesian) {
return [atan2(cartesian[1], cartesian[0]), asin(cartesian[2])];
}
function cartesian(spherical) {
var lambda = spherical[0], phi = spherical[1], cosPhi = cos(phi);
return [cosPhi * cos(lambda), cosPhi * sin(lambda), sin(phi)];
}
function cartesianDot(a, b) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
function cartesianCross(a, b) {
return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
}
// TODO return a
function cartesianAddInPlace(a, b) {
a[0] += b[0], a[1] += b[1], a[2] += b[2];
}
function cartesianScale(vector, k) {
return [vector[0] * k, vector[1] * k, vector[2] * k];
}
// TODO return d
function cartesianNormalizeInPlace(d) {
var l = sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
d[0] /= l, d[1] /= l, d[2] /= l;
}
function cartesianEqual(a, b) {
var dx = b[0] - a[0],
dy = b[1] - a[1],
dz = b[2] - a[2];
return dx * dx + dy * dy + dz * dz < epsilon2 * epsilon2;
}
var lambda0$1;
var phi0;
var lambda1;
var phi1;
var lambda2;
var lambda00$1;
var phi00$1;
var p0;
var deltaSum = adder();
var ranges;
var range$1;
var boundsStream = {
point: boundsPoint,
lineStart: boundsLineStart,
lineEnd: boundsLineEnd,
polygonStart: function() {
boundsStream.point = boundsRingPoint;
boundsStream.lineStart = boundsRingStart;
boundsStream.lineEnd = boundsRingEnd;
deltaSum.reset();
areaStream.polygonStart();
},
polygonEnd: function() {
areaStream.polygonEnd();
boundsStream.point = boundsPoint;
boundsStream.lineStart = boundsLineStart;
boundsStream.lineEnd = boundsLineEnd;
if (areaRingSum < 0) lambda0$1 = -(lambda1 = 180), phi0 = -(phi1 = 90);
else if (deltaSum > epsilon) phi1 = 90;
else if (deltaSum < -epsilon) phi0 = -90;
range$1[0] = lambda0$1, range$1[1] = lambda1;
}
};
function boundsPoint(lambda, phi) {
ranges.push(range$1 = [lambda0$1 = lambda, lambda1 = lambda]);
if (phi < phi0) phi0 = phi;
if (phi > phi1) phi1 = phi;
}
function linePoint(lambda, phi) {
var p = cartesian([lambda * radians, phi * radians]);
if (p0) {
var normal = cartesianCross(p0, p),
equatorial = [normal[1], -normal[0], 0],
inflection = cartesianCross(equatorial, normal);
cartesianNormalizeInPlace(inflection);
inflection = spherical(inflection);
var delta = lambda - lambda2,
sign = delta > 0 ? 1 : -1,
lambdai = inflection[0] * degrees * sign,
phii,
antimeridian = abs(delta) > 180;
if (antimeridian ^ (sign * lambda2 < lambdai && lambdai < sign * lambda)) {
phii = inflection[1] * degrees;
if (phii > phi1) phi1 = phii;
} else if (lambdai = (lambdai + 360) % 360 - 180, antimeridian ^ (sign * lambda2 < lambdai && lambdai < sign * lambda)) {
phii = -inflection[1] * degrees;
if (phii < phi0) phi0 = phii;
} else {
if (phi < phi0) phi0 = phi;
if (phi > phi1) phi1 = phi;
}
if (antimeridian) {
if (lambda < lambda2) {
if (angle(lambda0$1, lambda) > angle(lambda0$1, lambda1)) lambda1 = lambda;
} else {
if (angle(lambda, lambda1) > angle(lambda0$1, lambda1)) lambda0$1 = lambda;
}
} else {
if (lambda1 >= lambda0$1) {
if (lambda < lambda0$1) lambda0$1 = lambda;
if (lambda > lambda1) lambda1 = lambda;
} else {
if (lambda > lambda2) {
if (angle(lambda0$1, lambda) > angle(lambda0$1, lambda1)) lambda1 = lambda;
} else {
if (angle(lambda, lambda1) > angle(lambda0$1, lambda1)) lambda0$1 = lambda;
}
}
}
} else {
ranges.push(range$1 = [lambda0$1 = lambda, lambda1 = lambda]);
}
if (phi < phi0) phi0 = phi;
if (phi > phi1) phi1 = phi;
p0 = p, lambda2 = lambda;
}
function boundsLineStart() {
boundsStream.point = linePoint;
}
function boundsLineEnd() {
range$1[0] = lambda0$1, range$1[1] = lambda1;
boundsStream.point = boundsPoint;
p0 = null;
}
function boundsRingPoint(lambda, phi) {
if (p0) {
var delta = lambda - lambda2;
deltaSum.add(abs(delta) > 180 ? delta + (delta > 0 ? 360 : -360) : delta);
} else {
lambda00$1 = lambda, phi00$1 = phi;
}
areaStream.point(lambda, phi);
linePoint(lambda, phi);
}
function boundsRingStart() {
areaStream.lineStart();
}
function boundsRingEnd() {
boundsRingPoint(lambda00$1, phi00$1);
areaStream.lineEnd();
if (abs(deltaSum) > epsilon) lambda0$1 = -(lambda1 = 180);
range$1[0] = lambda0$1, range$1[1] = lambda1;
p0 = null;
}
// Finds the left-right distance between two longitudes.
// This is almost the same as (lambda1 - lambda0 + 360°) % 360°, except that we want
// the distance between ±180° to be 360°.
function angle(lambda0, lambda1) {
return (lambda1 -= lambda0) < 0 ? lambda1 + 360 : lambda1;
}
function rangeCompare(a, b) {
return a[0] - b[0];
}
function rangeContains(range, x) {
return range[0] <= range[1] ? range[0] <= x && x <= range[1] : x < range[0] || range[1] < x;
}
function bounds(feature) {
var i, n, a, b, merged, deltaMax, delta;
phi1 = lambda1 = -(lambda0$1 = phi0 = Infinity);
ranges = [];
geoStream(feature, boundsStream);
// First, sort ranges by their minimum longitudes.
if (n = ranges.length) {
ranges.sort(rangeCompare);
// Then, merge any ranges that overlap.
for (i = 1, a = ranges[0], merged = [a]; i < n; ++i) {
b = ranges[i];
if (rangeContains(a, b[0]) || rangeContains(a, b[1])) {
if (angle(a[0], b[1]) > angle(a[0], a[1])) a[1] = b[1];
if (angle(b[0], a[1]) > angle(a[0], a[1])) a[0] = b[0];
} else {
merged.push(a = b);
}
}
// Finally, find the largest gap between the merged ranges.
// The final bounding box will be the inverse of this gap.
for (deltaMax = -Infinity, n = merged.length - 1, i = 0, a = merged[n]; i <= n; a = b, ++i) {
b = merged[i];
if ((delta = angle(a[1], b[0])) > deltaMax) deltaMax = delta, lambda0$1 = b[0], lambda1 = a[1];
}
}
ranges = range$1 = null;
return lambda0$1 === Infinity || phi0 === Infinity
? [[NaN, NaN], [NaN, NaN]]
: [[lambda0$1, phi0], [lambda1, phi1]];
}
var W0;
var W1;
var X0;
var Y0;
var Z0;
var X1;
var Y1;
var Z1;
var X2;
var Y2;
var Z2;
var lambda00$2;
var phi00$2;
var x0;
var y0;
var z0;
// previous point
var centroidStream = {
sphere: noop,
point: centroidPoint,
lineStart: centroidLineStart,
lineEnd: centroidLineEnd,
polygonStart: function() {
centroidStream.lineStart = centroidRingStart;
centroidStream.lineEnd = centroidRingEnd;
},
polygonEnd: function() {
centroidStream.lineStart = centroidLineStart;
centroidStream.lineEnd = centroidLineEnd;
}
};
// Arithmetic mean of Cartesian vectors.
function centroidPoint(lambda, phi) {
lambda *= radians, phi *= radians;
var cosPhi = cos(phi);
centroidPointCartesian(cosPhi * cos(lambda), cosPhi * sin(lambda), sin(phi));
}
function centroidPointCartesian(x, y, z) {
++W0;
X0 += (x - X0) / W0;
Y0 += (y - Y0) / W0;
Z0 += (z - Z0) / W0;
}
function centroidLineStart() {
centroidStream.point = centroidLinePointFirst;
}
function centroidLinePointFirst(lambda, phi) {
lambda *= radians, phi *= radians;
var cosPhi = cos(phi);
x0 = cosPhi * cos(lambda);
y0 = cosPhi * sin(lambda);
z0 = sin(phi);
centroidStream.point = centroidLinePoint;
centroidPointCartesian(x0, y0, z0);
}
function centroidLinePoint(lambda, phi) {
lambda *= radians, phi *= radians;
var cosPhi = cos(phi),
x = cosPhi * cos(lambda),
y = cosPhi * sin(lambda),
z = sin(phi),
w = atan2(sqrt((w = y0 * z - z0 * y) * w + (w = z0 * x - x0 * z) * w + (w = x0 * y - y0 * x) * w), x0 * x + y0 * y + z0 * z);
W1 += w;
X1 += w * (x0 + (x0 = x));
Y1 += w * (y0 + (y0 = y));
Z1 += w * (z0 + (z0 = z));
centroidPointCartesian(x0, y0, z0);
}
function centroidLineEnd() {
centroidStream.point = centroidPoint;
}
// See J. E. Brock, The Inertia Tensor for a Spherical Triangle,
// J. Applied Mechanics 42, 239 (1975).
function centroidRingStart() {
centroidStream.point = centroidRingPointFirst;
}
function centroidRingEnd() {
centroidRingPoint(lambda00$2, phi00$2);
centroidStream.point = centroidPoint;
}
function centroidRingPointFirst(lambda, phi) {
lambda00$2 = lambda, phi00$2 = phi;
lambda *= radians, phi *= radians;
centroidStream.point = centroidRingPoint;
var cosPhi = cos(phi);
x0 = cosPhi * cos(lambda);
y0 = cosPhi * sin(lambda);
z0 = sin(phi);
centroidPointCartesian(x0, y0, z0);
}
function centroidRingPoint(lambda, phi) {
lambda *= radians, phi *= radians;
var cosPhi = cos(phi),
x = cosPhi * cos(lambda),
y = cosPhi * sin(lambda),
z = sin(phi),
cx = y0 * z - z0 * y,
cy = z0 * x - x0 * z,
cz = x0 * y - y0 * x,
m = sqrt(cx * cx + cy * cy + cz * cz),
w = asin(m), // line weight = angle
v = m && -w / m; // area weight multiplier
X2 += v * cx;
Y2 += v * cy;
Z2 += v * cz;
W1 += w;
X1 += w * (x0 + (x0 = x));
Y1 += w * (y0 + (y0 = y));
Z1 += w * (z0 + (z0 = z));
centroidPointCartesian(x0, y0, z0);
}
function centroid(object) {
W0 = W1 =
X0 = Y0 = Z0 =
X1 = Y1 = Z1 =
X2 = Y2 = Z2 = 0;
geoStream(object, centroidStream);
var x = X2,
y = Y2,
z = Z2,
m = x * x + y * y + z * z;
// If the area-weighted ccentroid is undefined, fall back to length-weighted ccentroid.
if (m < epsilon2) {
x = X1, y = Y1, z = Z1;
// If the feature has zero length, fall back to arithmetic mean of point vectors.
if (W1 < epsilon) x = X0, y = Y0, z = Z0;
m = x * x + y * y + z * z;
// If the feature still has an undefined ccentroid, then return.
if (m < epsilon2) return [NaN, NaN];
}
return [atan2(y, x) * degrees, asin(z / sqrt(m)) * degrees];
}
function constant(x) {
return function() {
return x;
};
}
function compose(a, b) {
function compose(x, y) {
return x = a(x, y), b(x[0], x[1]);
}
if (a.invert && b.invert) compose.invert = function(x, y) {
return x = b.invert(x, y), x && a.invert(x[0], x[1]);
};
return compose;
}
function rotationIdentity(lambda, phi) {
return [lambda > pi ? lambda - tau : lambda < -pi ? lambda + tau : lambda, phi];
}
rotationIdentity.invert = rotationIdentity;
function rotateRadians(deltaLambda, deltaPhi, deltaGamma) {
return (deltaLambda %= tau) ? (deltaPhi || deltaGamma ? compose(rotationLambda(deltaLambda), rotationPhiGamma(deltaPhi, deltaGamma))
: rotationLambda(deltaLambda))
: (deltaPhi || deltaGamma ? rotationPhiGamma(deltaPhi, deltaGamma)
: rotationIdentity);
}
function forwardRotationLambda(deltaLambda) {
return function(lambda, phi) {
return lambda += deltaLambda, [lambda > pi ? lambda - tau : lambda < -pi ? lambda + tau : lambda, phi];
};
}
function rotationLambda(deltaLambda) {
var rotation = forwardRotationLambda(deltaLambda);
rotation.invert = forwardRotationLambda(-deltaLambda);
return rotation;
}
function rotationPhiGamma(deltaPhi, deltaGamma) {
var cosDeltaPhi = cos(deltaPhi),
sinDeltaPhi = sin(deltaPhi),
cosDeltaGamma = cos(deltaGamma),
sinDeltaGamma = sin(deltaGamma);
function rotation(lambda, phi) {
var cosPhi = cos(phi),
x = cos(lambda) * cosPhi,
y = sin(lambda) * cosPhi,
z = sin(phi),
k = z * cosDeltaPhi + x * sinDeltaPhi;
return [
atan2(y * cosDeltaGamma - k * sinDeltaGamma, x * cosDeltaPhi - z * sinDeltaPhi),
asin(k * cosDeltaGamma + y * sinDeltaGamma)
];
}
rotation.invert = function(lambda, phi) {
var cosPhi = cos(phi),
x = cos(lambda) * cosPhi,
y = sin(lambda) * cosPhi,
z = sin(phi),
k = z * cosDeltaGamma - y * sinDeltaGamma;
return [
atan2(y * cosDeltaGamma + z * sinDeltaGamma, x * cosDeltaPhi + k * sinDeltaPhi),
asin(k * cosDeltaPhi - x * sinDeltaPhi)
];
};
return rotation;
}
function rotation(rotate) {
rotate = rotateRadians(rotate[0] * radians, rotate[1] * radians, rotate.length > 2 ? rotate[2] * radians : 0);
function forward(coordinates) {
coordinates = rotate(coordinates[0] * radians, coordinates[1] * radians);
return coordinates[0] *= degrees, coordinates[1] *= degrees, coordinates;
}
forward.invert = function(coordinates) {
coordinates = rotate.invert(coordinates[0] * radians, coordinates[1] * radians);
return coordinates[0] *= degrees, coordinates[1] *= degrees, coordinates;
};
return forward;
}
// Generates a circle centered at [0°, 0°], with a given radius and precision.
function circleStream(stream, radius, delta, direction, t0, t1) {
if (!delta) return;
var cosRadius = cos(radius),
sinRadius = sin(radius),
step = direction * delta;
if (t0 == null) {
t0 = radius + direction * tau;
t1 = radius - step / 2;
} else {
t0 = circleRadius(cosRadius, t0);
t1 = circleRadius(cosRadius, t1);
if (direction > 0 ? t0 < t1 : t0 > t1) t0 += direction * tau;
}
for (var point, t = t0; direction > 0 ? t > t1 : t < t1; t -= step) {
point = spherical([cosRadius, -sinRadius * cos(t), -sinRadius * sin(t)]);
stream.point(point[0], point[1]);
}
}
// Returns the signed angle of a cartesian point relative to [cosRadius, 0, 0].
function circleRadius(cosRadius, point) {
point = cartesian(point), point[0] -= cosRadius;
cartesianNormalizeInPlace(point);
var radius = acos(-point[1]);
return ((-point[2] < 0 ? -radius : radius) + tau - epsilon) % tau;
}
function circle() {
var center = constant([0, 0]),
radius = constant(90),
precision = constant(6),
ring,
rotate,
stream = {point: point};
function point(x, y) {
ring.push(x = rotate(x, y));
x[0] *= degrees, x[1] *= degrees;
}
function circle() {
var c = center.apply(this, arguments),
r = radius.apply(this, arguments) * radians,
p = precision.apply(this, arguments) * radians;
ring = [];
rotate = rotateRadians(-c[0] * radians, -c[1] * radians, 0).invert;
circleStream(stream, r, p, 1);
c = {type: "Polygon", coordinates: [ring]};
ring = rotate = null;
return c;
}
circle.center = function(_) {
return arguments.length ? (center = typeof _ === "function" ? _ : constant([+_[0], +_[1]]), circle) : center;
};
circle.radius = function(_) {
return arguments.length ? (radius = typeof _ === "function" ? _ : constant(+_), circle) : radius;
};
circle.precision = function(_) {
return arguments.length ? (precision = typeof _ === "function" ? _ : constant(+_), circle) : precision;
};
return circle;
}
function clipBuffer() {
var lines = [],
line;
return {
point: function(x, y, i, t) {
var point = [x, y];
// when called by clipPolygon, store index and t
if (arguments.length > 2) { point.index = i; point.t = t; }
line.push(point);
},
lineStart: function() {
lines.push(line = []);
},
lineEnd: noop,
rejoin: function() {
if (lines.length > 1) lines.push(lines.pop().concat(lines.shift()));
},
result: function() {
var result = lines;
lines = [];
line = null;
return result;
}
};
}
function clipLine(a, b, x0, y0, x1, y1) {
var ax = a[0],
ay = a[1],
bx = b[0],
by = b[1],
t0 = 0,
t1 = 1,
dx = bx - ax,
dy = by - ay,
r;
r = x0 - ax;
if (!dx && r > 0) return;
r /= dx;
if (dx < 0) {
if (r < t0) return;
if (r < t1) t1 = r;
} else if (dx > 0) {
if (r > t1) return;
if (r > t0) t0 = r;
}
r = x1 - ax;
if (!dx && r < 0) return;
r /= dx;
if (dx < 0) {
if (r > t1) return;
if (r > t0) t0 = r;
} else if (dx > 0) {
if (r < t0) return;
if (r < t1) t1 = r;
}
r = y0 - ay;
if (!dy && r > 0) return;
r /= dy;
if (dy < 0) {
if (r < t0) return;
if (r < t1) t1 = r;
} else if (dy > 0) {
if (r > t1) return;
if (r > t0) t0 = r;
}
r = y1 - ay;
if (!dy && r < 0) return;
r /= dy;
if (dy < 0) {
if (r > t1) return;
if (r > t0) t0 = r;
} else if (dy > 0) {
if (r < t0) return;
if (r < t1) t1 = r;
}
if (t0 > 0) a[0] = ax + t0 * dx, a[1] = ay + t0 * dy;
if (t1 < 1) b[0] = ax + t1 * dx, b[1] = ay + t1 * dy;
return true;
}
function pointEqual(a, b) {
return abs(a[0] - b[0]) < epsilon && abs(a[1] - b[1]) < epsilon;
}
function Intersection(point, points, other, entry) {
this.x = point;
this.z = points;
this.o = other; // another intersection
this.e = entry; // is an entry?
this.v = false; // visited
this.n = this.p = null; // next & previous
}
// A generalized polygon clipping algorithm: given a polygon that has been cut
// into its visible line segments, and rejoins the segments by interpolating
// along the clip edge.
function clipRejoin(segments, compareIntersection, startInside, interpolate, stream) {
var subject = [],
clip = [],
i,
n;
segments.forEach(function(segment) {
if ((n = segment.length - 1) <= 0) return;
var n, p0 = segment[0], p1 = segment[n], x;
// If the first and last points of a segment are coincident, then treat as a
// closed ring. TODO if all rings are closed, then the winding order of the
// exterior ring should be checked.
if (pointEqual(p0, p1)) {
stream.lineStart();
for (i = 0; i < n; ++i) stream.point((p0 = segment[i])[0], p0[1]);
stream.lineEnd();
return;
}
subject.push(x = new Intersection(p0, segment, null, true));
clip.push(x.o = new Intersection(p0, null, x, false));
subject.push(x = new Intersection(p1, segment, null, false));
clip.push(x.o = new Intersection(p1, null, x, true));
});
if (!subject.length) return;
clip.sort(compareIntersection);
link(subject);
link(clip);
for (i = 0, n = clip.length; i < n; ++i) {
clip[i].e = startInside = !startInside;
}
var start = subject[0],
points,
point;
while (1) {
// Find first unvisited intersection.
var current = start,
isSubject = true;
while (current.v) if ((current = current.n) === start) return;
points = current.z;
stream.lineStart();
do {
current.v = current.o.v = true;
if (current.e) {
if (isSubject) {
for (i = 0, n = points.length; i < n; ++i) stream.point((point = points[i])[0], point[1]);
} else {
interpolate(current.x, current.n.x, 1, stream);
}
current = current.n;
} else {
if (isSubject) {
points = current.p.z;
for (i = points.length - 1; i >= 0; --i) stream.point((point = points[i])[0], point[1]);
} else {
interpolate(current.x, current.p.x, -1, stream);
}
current = current.p;
}
current = current.o;
points = current.z;
isSubject = !isSubject;
} while (!current.v);
stream.lineEnd();
}
}
function link(array) {
if (!(n = array.length)) return;
var n,
i = 0,
a = array[0],
b;
while (++i < n) {
a.n = b = array[i];
b.p = a;
a = b;
}
a.n = b = array[0];
b.p = a;
}
var clipMax = 1e9;
var clipMin = -clipMax;
// TODO Use d3-polygon’s polygonContains here for the ring check?
// TODO Eliminate duplicate buffering in clipBuffer and polygon.push?
function clipExtent(x0, y0, x1, y1) {
function visible(x, y) {
return x0 <= x && x <= x1 && y0 <= y && y <= y1;
}
function interpolate(from, to, direction, stream) {
var a = 0, a1 = 0;
if (from == null
|| (a = corner(from, direction)) !== (a1 = corner(to, direction))
|| comparePoint(from, to) < 0 ^ direction > 0) {
do stream.point(a === 0 || a === 3 ? x0 : x1, a > 1 ? y1 : y0);
while ((a = (a + direction + 4) % 4) !== a1);
} else {
stream.point(to[0], to[1]);
}
}
function corner(p, direction) {
return abs(p[0] - x0) < epsilon ? direction > 0 ? 0 : 3
: abs(p[0] - x1) < epsilon ? direction > 0 ? 2 : 1
: abs(p[1] - y0) < epsilon ? direction > 0 ? 1 : 0
: direction > 0 ? 3 : 2; // abs(p[1] - y1) < epsilon
}
function compareIntersection(a, b) {
return comparePoint(a.x, b.x);
}
function comparePoint(a, b) {
var ca = corner(a, 1),
cb = corner(b, 1);
return ca !== cb ? ca - cb
: ca === 0 ? b[1] - a[1]
: ca === 1 ? a[0] - b[0]
: ca === 2 ? a[1] - b[1]
: b[0] - a[0];
}
return function(stream) {
var activeStream = stream,
bufferStream = clipBuffer(),
segments,
polygon,
ring,
x__, y__, v__, // first point
x_, y_, v_, // previous point
first,
clean;
var clipStream = {
point: point,
lineStart: lineStart,
lineEnd: lineEnd,
polygonStart: polygonStart,
polygonEnd: polygonEnd
};
function point(x, y) {
if (visible(x, y)) activeStream.point(x, y);
}
function polygonInside() {
var winding = 0;
for (var i = 0, n = polygon.length; i < n; ++i) {
for (var ring = polygon[i], j = 1, m = ring.length, point = ring[0], a0, a1, b0 = point[0], b1 = point[1]; j < m; ++j) {
a0 = b0, a1 = b1, point = ring[j], b0 = point[0], b1 = point[1];
if (a1 <= y1) { if (b1 > y1 && (b0 - a0) * (y1 - a1) > (b1 - a1) * (x0 - a0)) ++winding; }
else { if (b1 <= y1 && (b0 - a0) * (y1 - a1) < (b1 - a1) * (x0 - a0)) --winding; }
}
}
return winding;
}
// Buffer geometry within a polygon and then clip it en masse.
function polygonStart() {
activeStream = bufferStream, segments = [], polygon = [], clean = true;
}
function polygonEnd() {
var startInside = polygonInside(),
cleanInside = clean && startInside,
visible = (segments = d3Array.merge(segments)).length;
if (cleanInside || visible) {
stream.polygonStart();
if (cleanInside) {
stream.lineStart();
interpolate(null, null, 1, stream);
stream.lineEnd();
}
if (visible) {
clipRejoin(segments, compareIntersection, startInside, interpolate, stream);
}
stream.polygonEnd();
}
activeStream = stream, segments = polygon = ring = null;
}
function lineStart() {
clipStream.point = linePoint;
if (polygon) polygon.push(ring = []);
first = true;
v_ = false;
x_ = y_ = NaN;
}
// TODO rather than special-case polygons, simply handle them separately.
// Ideally, coincident intersection points should be jittered to avoid
// clipping issues.
function lineEnd() {
if (segments) {
linePoint(x__, y__);
if (v__ && v_) bufferStream.rejoin();
segments.push(bufferStream.result());
}
clipStream.point = point;
if (v_) activeStream.lineEnd();
}
function linePoint(x, y) {
var v = visible(x, y);
if (polygon) ring.push([x, y]);
if (first) {
x__ = x, y__ = y, v__ = v;
first = false;
if (v) {
activeStream.lineStart();
activeStream.point(x, y);
}
} else {
if (v && v_) activeStream.point(x, y);
else {
var a = [x_ = Math.max(clipMin, Math.min(clipMax, x_)), y_ = Math.max(clipMin, Math.min(clipMax, y_))],
b = [x = Math.max(clipMin, Math.min(clipMax, x)), y = Math.max(clipMin, Math.min(clipMax, y))];
if (clipLine(a, b, x0, y0, x1, y1)) {
if (!v_) {
activeStream.lineStart();
activeStream.point(a[0], a[1]);
}
activeStream.point(b[0], b[1]);
if (!v) activeStream.lineEnd();
clean = false;
} else if (v) {
activeStream.lineStart();
activeStream.point(x, y);
clean = false;
}
}
}
x_ = x, y_ = y, v_ = v;
}
return clipStream;
};
}
function extent() {
var x0 = 0,
y0 = 0,
x1 = 960,
y1 = 500,
cache,
cacheStream,
clip;
return clip = {
stream: function(stream) {
return cache && cacheStream === stream ? cache : cache = clipExtent(x0, y0, x1, y1)(cacheStream = stream);
},
extent: function(_) {
return arguments.length ? (x0 = +_[0][0], y0 = +_[0][1], x1 = +_[1][0], y1 = +_[1][1], cache = cacheStream = null, clip) : [[x0, y0], [x1, y1]];
}
};
}
var sum = adder();
function polygonContains(polygon, point) {
var lambda = point[0],
phi = point[1],
normal = [sin(lambda), -cos(lambda), 0],
angle = 0,
winding = 0;
sum.reset();
for (var i = 0, n = polygon.length; i < n; ++i) {
if (!(m = (ring = polygon[i]).length)) continue;
var ring,
m,
point0 = ring[m - 1],
lambda0 = point0[0],
phi0 = point0[1] / 2 + quarterPi,
sinPhi0 = sin(phi0),
cosPhi0 = cos(phi0);
for (var j = 0; j < m; ++j, lambda0 = lambda1, sinPhi0 = sinPhi1, cosPhi0 = cosPhi1, point0 = point1) {
var point1 = ring[j],
lambda1 = point1[0],
phi1 = point1[1] / 2 + quarterPi,
sinPhi1 = sin(phi1),
cosPhi1 = cos(phi1),
delta = lambda1 - lambda0,
sign = delta >= 0 ? 1 : -1,
absDelta = sign * delta,
antimeridian = absDelta > pi,
k = sinPhi0 * sinPhi1;
sum.add(atan2(k * sign * sin(absDelta), cosPhi0 * cosPhi1 + k * cos(absDelta)));
angle += antimeridian ? delta + sign * tau : delta;
// Are the longitudes either side of the point’s meridian (lambda),
// and are the latitudes smaller than the parallel (phi)?
if (antimeridian ^ lambda0 >= lambda ^ lambda1 >= lambda) {
var arc = cartesianCross(cartesian(point0), cartesian(point1));
cartesianNormalizeInPlace(arc);
var intersection = cartesianCross(normal, arc);
cartesianNormalizeInPlace(intersection);
var phiArc = (antimeridian ^ delta >= 0 ? -1 : 1) * asin(intersection[2]);
if (phi > phiArc || phi === phiArc && (arc[0] || arc[1])) {
winding += antimeridian ^ delta >= 0 ? 1 : -1;
}
}
}
}
// First, determine whether the South pole is inside or outside:
//
// It is inside if:
// * the polygon winds around it in a clockwise direction.
// * the polygon does not (cumulatively) wind around it, but has a negative
// (counter-clockwise) area.
//
// Second, count the (signed) number of times a segment crosses a lambda
// from the point to the South pole. If it is zero, then the point is the
// same side as the South pole.
return (angle < -epsilon || angle < epsilon && sum < -epsilon) ^ (winding & 1);
}
var lengthSum = adder();
var lambda0$2;
var sinPhi0$1;
var cosPhi0$1;
var lengthStream = {
sphere: noop,
point: noop,
lineStart: lengthLineStart,
lineEnd: noop,
polygonStart: noop,
polygonEnd: noop
};
function lengthLineStart() {
lengthStream.point = lengthPointFirst;
lengthStream.lineEnd = lengthLineEnd;
}
function lengthLineEnd() {
lengthStream.point = lengthStream.lineEnd = noop;
}
function lengthPointFirst(lambda, phi) {
lambda *= radians, phi *= radians;
lambda0$2 = lambda, sinPhi0$1 = sin(phi), cosPhi0$1 = cos(phi);
lengthStream.point = lengthPoint;
}
function lengthPoint(lambda, phi) {
lambda *= radians, phi *= radians;
var sinPhi = sin(phi),
cosPhi = cos(phi),
delta = abs(lambda - lambda0$2),
cosDelta = cos(delta),
sinDelta = sin(delta),
x = cosPhi * sinDelta,
y = cosPhi0$1 * sinPhi - sinPhi0$1 * cosPhi * cosDelta,
z = sinPhi0$1 * sinPhi + cosPhi0$1 * cosPhi * cosDelta;
lengthSum.add(atan2(sqrt(x * x + y * y), z));
lambda0$2 = lambda, sinPhi0$1 = sinPhi, cosPhi0$1 = cosPhi;
}
function length(object) {
lengthSum.reset();
geoStream(object, lengthStream);
return +lengthSum;
}
var coordinates = [null, null];
var object = {type: "LineString", coordinates: coordinates};
function distance(a, b) {
coordinates[0] = a;
coordinates[1] = b;
return length(object);
}
var containsObjectType = {
Feature: function(object, point) {
return containsGeometry(object.geometry, point);
},
FeatureCollection: function(object, point) {
Raw
index.html
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Raw
kruskal.js
View raw

(Sorry about that, but we can’t show files that are this big right now.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment