Computing the inverse of the ChamberlinAfrica projection for issue #85.
When a projection does not have an invert
function, we apply the Newton-Raphson method (see gist 8903368) on a function of the cartesian coordinates (x,y,z). This avoids lots of problems at the poles.
This approach seems to work well with:
It doesn't work with
The main code is the following:
function solve(project, precision) {
var n = 100,
step = (halfPi - epsilon) / n,
i,
j,
grid = [];
for (i = 0; i <= 4 * n; i++) {
grid[i] = [];
for (j = 0; j <= 2 * n; j++) {
grid[i][j] = project((i - 2 * n) * step, (j - n) * step);
}
}
function invert (x, y) {
// find a start point c "close enough" to x,y
var i, j,
c,
m, min = +Infinity;
// d3.scan
for (i = 0; i <= 4 * n; i++) {
for (j = 0; j <= 2 * n; j++) {
m = abs(x - grid[i][j][0]) + abs(y - grid[i][j][1]);
if (m < min) {
c = [i,j];
min = m;
}
}
}
c = [ step * (c[0] - 2 * n), step * (c[1] - n) ];
c = cartesian(c);
// solve for x,y
var solution = Newton.Solve(g => {
var norm = g[0] * g[0] + g[1] * g[1] + g[2] * g[2];
cartesianScale(g, 1 / sqrt(norm));
var s = spherical(g),
p = project(s[0], s[1]);
return [p[0], p[1], norm];
}, [x, y, 1], { start: c, acc: precision });
var norm = solution[0] * solution[0] + solution[1] * solution[1] + solution[2] * solution[2];
cartesianScale(solution, 1 / sqrt(norm));
return spherical(solution);
}
return invert;
}
/*
* Newton's method for finding roots
*
* code adapted from D.V. Fedorov,
* “Introduction to Numerical Methods with examples in Javascript”
* http://owww.phys.au.dk/~fedorov/nucltheo/Numeric/11/book.pdf
* (licensed under the GPL)
* by Philippe Riviere <philippe.riviere@illisible.net> March 2014
* modified for compatibility with Chrome/Safari
* added a max iterations parameter
*
* Usage: Newton.Solve(Math.exp, 2)); => ~ log(2)
* Newton.Solve(d3.geo.chamberlin(), [200,240])
*/
var Newton = {version: "1.0.0"}; // semver
Newton.Norm = function (v) {
return Math.sqrt(v.reduce(function (s, e) {
return s + e * e
}, 0));
}
Newton.Dot = function (a, b) {
var s = 0;
for (var i in a) s += a[i] * b[i];
return s;
}
// QR - decomposition A=QR of matrix A
Newton.QRDec = function(A){
var m = A.length, R = [], i, j, k;
for (j in A) {
R[j] = [];
for (i in A) R[j][i] = 0;
}
var Q = [];
for (i in A) {
Q[i] = [];
for (j in A[0]) Q[i][j] = A[i][j];
}
// Q is a copy of A
for (i = 0; i < m; i++) {
var e = Q[i],
r = Math.sqrt(Newton.Dot(e, e));
if (r == 0) throw "Newton.QRDec: singular matrix"
R[i][i] = r;
for (k in e) e[k] /= r;
// normalization
for (j = i + 1; j < m; j++) {
var q = Q[j],
s = Newton.Dot(e, q);
for (k in q) q[k] -= s * e[k];
// orthogonalization
R[j][i] = s;
}
}
return [Q, R];
};
// QR - backsubstitution
// input: matrices Q,R, array b; output: array x such that QRx=b
Newton.QRBack = function(Q, R, b) {
var m = Q.length,
c = new Array(m),
x = new Array(m),
i, k, s;
for (i in Q) {
// c = QˆT b
c[i] = 0;
for (k in b) c[i] += Q[i][k] * b[k];
}
for (i = m - 1; i >= 0; i--) {
// back substitution
for (s = 0, k = i + 1; k < m; k++) s += R[k][i] * x[k];
x[i] = (c[i] - s) / R[i][i];
}
return x;
}
// calculates inverse of matrix A
Newton.Inverse = function(A){
var t = Newton.QRDec(A),
Q = t[0],
R = t[1];
var m = [], i, k, n;
for (i in A) {
n = [];
for (k in A) {
n[k] = (k == i ? 1 : 0)
}
m[i] = Newton.QRBack(Q, R, n);
}
return m;
};
Newton.Zero = function (fs, x, opts={} /* acc, dx, max */) {
// Newton's root-finding method
var i, j, k;
if (opts.acc == undefined) opts.acc = 1e-6
if (opts.dx == undefined) opts.dx = 1e-3
if (opts.max == undefined) opts.max = 40 // max iterations
var J = [];
for (i in x) {
J[i] = [];
for (j in x) J[i][j] = 0;
}
var minusfx = [];
var v = fs(x);
if (v == null) throw "unable to compute fs at "+JSON.stringify(x)
for (i in x) minusfx[i] = -v[i];
do {
if (opts.max-- < 0) return null;
for (i in x)
for (k in x) {
// calculate Jacobian
x[k] += opts.dx
v = fs(x);
if (v == null) throw "unable to compute fs at "+JSON.stringify(x)
J[k][i] = (v[i] + minusfx[i]) / opts.dx
x[k] -= opts.dx
}
var t = Newton.QRDec(J),
Q = t[0],
R = t[1],
Dx = Newton.QRBack(Q, R, minusfx)
// Newton's step
var s = 2
do {
// simple backtracking line search
s = s / 2;
var z = [];
for (i in x) {
z[i] = x[i] + s * Dx[i];
}
var minusfz = [];
v = fs(z);
if (v == null) throw "unable to compute fs at "+JSON.stringify(z)
for (i in x) {
minusfz[i] = -v[i];
}
}
while (Newton.Norm(minusfz) > (1 - s / 2) * Newton.Norm(minusfx) && s > 1. / 128)
minusfx = minusfz;
x = z;
// step done
}
while (Newton.Norm(minusfx) > opts.acc)
return x;
}
Newton.Solve = function(fs, res, opts={}){
if (typeof res != 'object') res = [ typeof res == 'number'
? + res
: 0
];
var _fs = fs;
fs = function(x) {
var r = _fs(x);
if (typeof r == 'number') r = [ r ];
for (var i in r) r[i] -= res[i];
return r;
}
var start = [];
if (opts.start) {
start = opts.start;
} else {
for (var i in res) start[i]=0;
}
var n = Newton.Zero(fs, start, opts);
if (n && n.length == 1) n = n[0];
return n;
};
forked from Fil's block: Interrupted Homolosine
forked from Fil's block: ChamberlinAfrica inverse projection
xxxxxxxxxx
<meta charset="utf-8">
<style>
.stroke {
fill: none;
stroke: #000;
stroke-width: 3px;
}
.fill {
fill: #fff;
}
.graticule {
fill: none;
stroke: #777;
stroke-width: 0.5px;
stroke-opacity: 0.5;
}
.land {
fill: #222;
}
.boundary {
fill: none;
stroke: #fff;
stroke-width: 0.5px;
}
</style>
<svg width="960" height="484"></svg>
<script src="https://d3js.org/d3.v4.min.js"></script>
<script src="d3-geo.js"></script>
<script src="geoNewtonInverse.js"></script>
<script src="https://d3js.org/d3-geo-projection.v2.js"></script>
<script src="https://d3js.org/topojson.v1.min.js"></script>
<script>
var svg = d3.select("svg"),
width = +svg.attr("width"),
height = +svg.attr("height");
// import from math
var sin = Math.sin, cos = Math.cos, pi = Math.PI;
// import from d3-geo-projection
var hammer = d3.geoHammerRaw(1.68, 2);
var visionscarto_bertin_1953_alpha3 = function (lambda, phi) {
var fu = 1.4;
if (lambda + phi < -fu) {
var u = (lambda - phi + 1.6) * (lambda + phi + fu) / 8;
lambda += u;
phi -= 0.8 * u * sin(phi + pi/2);
}
var r = hammer(lambda, phi);
var k = 12,
d = (1 - cos(lambda * phi)) / k;
if (r[1] < 0) {
r[0] *= 1 + d;
}
if (r[1] > 0) {
r[1] *= 1 + d / 1.5 * r[0] * r[0];
}
return r;
};
(d3.geoVisionscarto = (function () {
visionscarto_bertin_1953_alpha3.invert = geoInverse(visionscarto_bertin_1953_alpha3, 0.01, d3.geoRobinsonRaw.invert);
return d3.geoProjection(visionscarto_bertin_1953_alpha3)
.rotate([-16.5, -42])
.scale(176.5);
}));
var k = 2 * Math.PI / 180; // min geoDistance in radians to be happy with the result
var farscale = d3.scaleSqrt()
.domain([0,k])
.range(['green', 'red'])
var projection = d3.geoVisionscarto();
var graticule = d3.geoGraticule();
var path = d3.geoPath()
.projection(projection);
var defs = svg.append("defs");
defs.append("path")
.datum({type: "Sphere"})
.attr("id", "sphere")
.attr("d", path);
defs.append("clipPath")
.attr("id", "clip")
.append("use")
.attr("xlink:href", "#sphere");
svg.append("use")
.attr("class", "stroke")
.attr("xlink:href", "#sphere");
svg.append("use")
.attr("class", "fill")
.attr("xlink:href", "#sphere");
svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("clip-path", "url(#clip)")
.attr("d", path);
var step = 5;
d3.range(-89, 89, step)
.forEach(lat => d3.range(-180, 180, step / Math.cos(lat * Math.PI / 180))
.forEach(lon => {
var o = [lon, lat],
m = projection(o);
var i = null;
try {
i = projection.invert(m);
} catch(e){
console.log('m', m, e)
}
svg.append("path")
.attr("fill", !i ? 'blue' : farscale(d3.geoDistance(o,i)))
.attr("d", path({type:"Point", coordinates: o}))
}))
var marker = svg.append("path")
.attr("fill", "orange")
svg.on('mousemove', function() {
var mouse = d3.mouse(this);
marker.attr("d", path({type:"Point", coordinates: projection.invert(mouse)}));
});
d3.json("https://gist.githubusercontent.com/mbostock/4090846/raw/d534aba169207548a8a3d670c9c2cc719ff05c47/world-50m.json", function(error, world) {
if (error) throw error;
svg.insert("path", ".graticule")
.datum(topojson.feature(world, world.objects.land))
.attr("class", "land")
.attr("clip-path", "url(#clip)")
.attr("d", path);
svg.insert("path", ".graticule")
.datum(topojson.mesh(world, world.objects.countries, function(a, b) { return a !== b; }))
.attr("class", "boundary")
.attr("clip-path", "url(#clip)")
.attr("d", path);
});
</script>
https://d3js.org/d3.v4.min.js
https://d3js.org/d3-geo-projection.v2.js
https://d3js.org/topojson.v1.min.js