block by fil 8783c504d2ebb736ce12bcfd1606e840

L’Équateur de Marseille

Full Screen

L’équateur de Marseille s’appuie sur la Canebière et la prolonge dans les deux directions, jusqu’à faire le tour du monde.

Cette approche a été conçue et développée par Jean-Luc Arnaud dans « Marseille et la Canebière - Echelles du monde », 2017. <halshs-01484193>

Visualisation interactive par Philippe Rivière – visionscarto.net

Forked from mbostock‘s block: Versor Dragging II

index.html

<!DOCTYPE html>
<html>
<head>
  <style>
    svg { display: none }
    #download { display: none; font-size: small; }
    #download { text-decoration: none; color: black; }
    body.svg svg { display: block }
    body.svg #download { display: inline-block }
    body.svg canvas { display: none }
  </style>
</head>
<body class="svg">
  <div id="svgwrap">    <svg version="1.1"
    	xmlns="//www.w3.org/2000/svg"
    	xmlns:xlink="//www.w3.org/1999/xlink"
    	xmlns:inkscape="//www.inkscape.org/namespaces/inkscape"
 width="960" height="600">
</svg></div>
  <canvas width="960" height="600"></canvas>
  <select style="position: absolute; top: 0; left: 0"></select>
  
  <label style="position: absolute; top: 1.4em; left: 0"><input name="svgcanvas" type="radio" checked="checked">SVG <a id="download" href="#">[download]</a></label>
  <label style="position: absolute; top: 2.4em; left: 0"><input  name="svgcanvas" type="radio">canvas</label>

<script src="https://d3js.org/d3.v4.min.js"></script>
<script src="https://d3js.org/d3-geo-projection.v2.js"></script>
<script src="https://d3js.org/topojson.v2.min.js"></script>
<script src="versor.js"></script>

<!--
<script src="d3.v4.min.js"></script>
<script src="d3-geo-projection.js"></script>
<script src="topojson.v2.min.js"></script>
<script src="versor.js"></script>
<script src="geoinverseNewton.js"></script>
-->
<script>
d3.selectAll('input').on('change', () => {
  d3.select('body').classed('svg', !d3.select('body').classed('svg'));
  render();
})

d3.select('#download').on('click', function(){
  var b64      = btoa(d3.select("#svgwrap").html());
  d3.select('#download')
    .attr('href', "data:image/svg+xml;base64,\n" + b64)
    .attr('download', "map.svg")
  	.dispatch('click')
})


var menu = d3.select('select')
.on('change', () => {
  var sel = menu.node().parentElement;
  projection = d3['geo' + sel.options[sel.selectedIndex].value]()
    .precision(0.1)
	.rotate([-canebiere[0], -canebiere[1], -angle]);
  render();
})
.selectAll('option')
.data([ 'Orthographic', 'Homolosine', 'Mollweide', 'Fahey', 'Sinusoidal', 'Miller', 'InterruptedMollweideHemispheres', 'Equirectangular', 'Mercator', 'Lagrange', 'Larrivee',  ])
.enter()
.append('option')
.attr('value', d => d)
.text(d => d);

// import from math
var sin = Math.sin, cos = Math.cos, pi = Math.PI;
var degrees = 180 / pi;

// retina display
var canvas = d3.select("canvas"),
    svg = d3.select("svg"),
  width = canvas.property("width"),
  height = canvas.property("height"),
  context = canvas.node().getContext("2d");

var devicePixelRatio = window.devicePixelRatio || 1;
canvas.attr('width', width * devicePixelRatio);
canvas.attr('height', height * devicePixelRatio);
canvas.style('width', width + 'px');
canvas.style('height', height + 'px');
context.setTransform(devicePixelRatio,0,0,devicePixelRatio,0,0)

canebiere = [ 5.3796269, 43.2971309 ];
angle = 25 + 44/60;  // 25°44'

antipode = [canebiere[0]+180, -canebiere[1]]

projection = d3.geoOrthographic()
	.precision(0.1)
	.rotate([-canebiere[0], -canebiere[1], -angle]);

  
// rotator = envoie Marseille vers le [0,0] de notre référentiel
// rotator.invert = envoie [0,0] vers Marseille
rotator = d3.geoRotation([-canebiere[0], -canebiere[1], -angle]);

  
equateur = d3.geoGraticule()
.extentMajor([[-180,-10], [180,10]])
.extentMinor([[-180,-5], [180,5]])
()

equateur.coordinates = equateur.coordinates.slice(4,5).map(d => d.map(rotator.invert))

var graticule = d3.geoGraticule()();

d3.json("110m.json", function(error, world) {
  if (error) throw error;
  var land = topojson.feature(world, world.objects["land"]);

  var defs = svg.append('defs');
  defs.append("clipPath")
    .attr("id", "clip")
    .append("use")
    .attr("xlink:href", "#sphere");
  defs
  .append('path')
  .datum({type: "Sphere"})
  .attr('id', 'sphere');

  svg.append('path')
  .datum(graticule)
  .attr('fill', 'none')
  .attr('stroke-width', 0.5)
  .attr('stroke', '#999')
  .style('clip-path', 'url(#clip)');

  svg.append('path')
  .datum(land)
  .style('clip-path', 'url(#clip)');

  svg.append('path')
  .datum(equateur)
  .attr('fill', 'none')
  .attr('stroke', 'green')
  .attr('stroke-width', 3)
  .style('clip-path', 'url(#clip)');

  svg.append('path')
  .datum({type: "Sphere"})
  .attr('fill', 'none')
  .attr('stroke', 'black');

  console.log([[0,0], [90,0], [180,0], [-90,0],[0,90],[0,-90]].map(rotator.invert))
  svg.append('path')
  .datum({
    type:"MultiPoint",
    coordinates: [[0,0], [90,0], [180,0], [-90,0],[0,90],[0,-90]].map(rotator.invert)
  })
  .attr('fill', 'orange')
  

  if (0) d3.timer(e => {
    projection.rotate([e/150, -canebiere[1]])
    svg.selectAll('path').attr('d', path)
  }, 100);

  
  
  render = function() {
      var path = d3.geoPath().projection(projection);

    if (d3.select('body').classed('svg')) {
      svg.selectAll('path,clipPath').attr('d', path);
      return;
    }
    
    path = path.context(context)
    context.save();

    (context.fillStyle = "#fff"), context.fillRect(0, 0, width, height);

    context.beginPath(), path({type: "Sphere"}), context.closePath(), context.clip();

    context.beginPath(), path(graticule), (context.lineWidth = 0.5, context.strokeStyle =
      "#999"), context.stroke(), context.lineWidth = 1, context.closePath();


    context.beginPath(), path(land), (context.strokeStyle = context.fillStyle =
      "#000"), context.fill(), context.lineWidth = 0.1, context.stroke(), context.lineWidth = 1, context.closePath();

    context.beginPath(), path(equateur), (context.lineWidth = 3, context.strokeStyle = "green", context.fillStyle =
      "rgba(0,0,0,0)"), context.fill(), context.stroke(), context.closePath(), context.lineWidth = 1;

    context.beginPath(), context.lineWidth = 2, path({type: "Sphere"}), (context.strokeStyle =
      "black"), context.stroke(), context.closePath();

    context.restore();
  };

  render();
});

canvas.call(d3.drag().on("start", dragstarted).on("drag", dragged));
svg.call(d3.drag().on("start", dragstarted).on("drag", dragged));

var render = function() {},
  v0, // Mouse position in Cartesian coordinates at start of drag gesture.
  r0, // Projection rotation as Euler angles at start.
  q0; // Projection rotation as versor at start.

function dragstarted() {
  v0 = versor.cartesian(projection.invert(d3.mouse(this)));
  r0 = projection.rotate();
  q0 = versor(r0);
}

function dragged() {
  var inv = projection.rotate(r0).invert(d3.mouse(this));
  if (isNaN(inv[0])) return;
  var v1 = versor.cartesian(inv),
    q1 = versor.multiply(q0, versor.delta(v0, v1)),
    r1 = versor.rotation(q1);
  projection.rotate(r1);
  render();
}

  

</script>
</body>
</html>

Makefile

topo-sans-antarctica:
	ndjson-cat ~/Sites/maps/some-geo-data/ne_110m_admin_0_countries/countries.geojson | ndjson-split 'd.features' | ndjson-filter 'd.properties.iso_a3 !== "ATA"' | ndjson-reduce 'p.features.push(d), p' '{type: "FeatureCollection", features: []}' | geo2topo > ~/Sites/d3/269920a4dd223bc48620f284e886cc6c/110m-sans-antarctica.json

geoInverseNewton.js

// import {abs, epsilon, halfPi, sqrt} from 'math'
var abs = Math.abs, epsilon = 1e-6, halfPi = Math.PI/2, sqrt = Math.sqrt;

// import {cartesian, cartesianScale, spherical} from 'cartesian'
var asin  =Math.asin, atan2 = Math.atan2, cos = Math.cos, sin = Math.sin;

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 cartesianScale(vector, k) {
  return [vector[0] * k, vector[1] * k, vector[2] * k];
}

function initial_general(project) {
      var n = 50,
      step = (halfPi - epsilon) / n,
      i,
      j,
      grid = [];
    for (i = 0; i <= 4 * n; i++) {
      grid[i] = [];
      for (j = 0; j <= 2 * n; j++) {
        var p = [(i - 2 * n) * step, (j - n) * step]
        grid[i][j] = project(p[0], p[1]);
      }
    }

  return function(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) ];
      return c;
  }
}

function geoInverse(project, precision, initial) {
  
  _initial = initial_general(project)
  _____initial = function(x,y) {
    if (y > 0.3) return ___initial(x+0.16,y-0.035);
    if (x > 0) return ___initial(x+0.16,y+0.015);
    return ___initial(x-0.1,y);
  }
  function invert (x,y) {
    // find a start point c
    var point = [x,y],
        c = _initial(x,y);
        
    c[0] *= 0.999;
    
    c = cartesian(c);
    // solve for x,y
    try {
      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], 1/norm];
        },
        [point[0], point[1], 1],
        { start: c, acc: precision, dx: 1e-5, max: 100 }
      );
    } catch(e) {
    	console.log(e);
    }
    if (solution) 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 x /* bad approximation better that nothing ?? */;
        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),
            no = Newton.Norm(Dx)
            // 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];// * Math.min(no, 1e-2);
            }
//console.log(z, no)
            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;
  };

versor.js

// Version 0.0.0. Copyright 2017 Mike Bostock.
(function(global, factory) {
  typeof exports === 'object' && typeof module !== 'undefined' ? module.exports = factory() :
  typeof define === 'function' && define.amd ? define(factory) :
  (global.versor = factory());
}(this, (function() {'use strict';

var acos = Math.acos,
    asin = Math.asin,
    atan2 = Math.atan2,
    cos = Math.cos,
    max = Math.max,
    min = Math.min,
    PI = Math.PI,
    sin = Math.sin,
    sqrt = Math.sqrt,
    radians = PI / 180,
    degrees = 180 / PI;

// Returns the unit quaternion for the given Euler rotation angles [λ, φ, γ].
function versor(e) {
  var l = e[0] / 2 * radians, sl = sin(l), cl = cos(l), // λ / 2
      p = e[1] / 2 * radians, sp = sin(p), cp = cos(p), // φ / 2
      g = e[2] / 2 * radians, sg = sin(g), cg = cos(g); // γ / 2
  return [
    cl * cp * cg + sl * sp * sg,
    sl * cp * cg - cl * sp * sg,
    cl * sp * cg + sl * cp * sg,
    cl * cp * sg - sl * sp * cg
  ];
}

// Returns Cartesian coordinates [x, y, z] given spherical coordinates [λ, φ].
versor.cartesian = function(e) {
  var l = e[0] * radians, p = e[1] * radians, cp = cos(p);
  return [cp * cos(l), cp * sin(l), sin(p)];
};

// Returns the Euler rotation angles [λ, φ, γ] for the given quaternion.
versor.rotation = function(q) {
  return [
    atan2(2 * (q[0] * q[1] + q[2] * q[3]), 1 - 2 * (q[1] * q[1] + q[2] * q[2])) * degrees,
    asin(max(-1, min(1, 2 * (q[0] * q[2] - q[3] * q[1])))) * degrees,
    atan2(2 * (q[0] * q[3] + q[1] * q[2]), 1 - 2 * (q[2] * q[2] + q[3] * q[3])) * degrees
  ];
};

// Returns the quaternion to rotate between two cartesian points on the sphere.
versor.delta = function(v0, v1) {
  var w = cross(v0, v1), l = sqrt(dot(w, w));
  if (!l) return [1, 0, 0, 0];
  var t = acos(max(-1, min(1, dot(v0, v1)))) / 2, s = sin(t); // t = θ / 2
  return [cos(t), w[2] / l * s, -w[1] / l * s, w[0] / l * s];
};

// Returns the quaternion that represents q0 * q1.
versor.multiply = function(q0, q1) {
  return [
    q0[0] * q1[0] - q0[1] * q1[1] - q0[2] * q1[2] - q0[3] * q1[3],
    q0[1] * q1[0] + q0[0] * q1[1] + q0[2] * q1[3] - q0[3] * q1[2],
    q0[0] * q1[2] - q0[1] * q1[3] + q0[2] * q1[0] + q0[3] * q1[1],
    q0[0] * q1[3] + q0[1] * q1[2] - q0[2] * q1[1] + q0[3] * q1[0]
  ];
};

function cross(v0, v1) {
  return [
    v0[1] * v1[2] - v0[2] * v1[1],
    v0[2] * v1[0] - v0[0] * v1[2],
    v0[0] * v1[1] - v0[1] * v1[0]
  ];
}

function dot(v0, v1) {
  return v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2];
}

return versor;
})));