Bertin1953 inverse projection [UNLISTED]

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:

ChamberlinAfrica
ModifiedStereographicAlaska
ModifiedStereographicGs48
ModifiedStereographicMiller
ModifiedStereographicLee
ModifiedStereographicGs50 (the part over Europe is broken)
TwoPointAzimuthalUsa

It doesn't work with

InterruptedSinuMollweide (problem of definition of the invert in interrupted.js)

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

<!DOCTYPE html>
<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);

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