Search code examples
d3.jsgisprojectiontopojsonmercator

D3.js : How to calculate accurately the Mercator transformation ratio `r`?


Given a square ABCD on Earth's surface (Equirectangular), with A & B on the Greenwitch meridian, B & C on meridian longitude = 10⁰ :

A( 0.0; 50.0)          C(10.0; 50.0)

B( 0.0; 40.0)          B(10.0; 40.0)

Given my D3js dataviz works in d3.geo.mercator() projection, my square is vertically transformed by a ratio r= mercator_height in px/width in px about 1.5.

How to calculate accurately this Mercator transformation ratio r ?

Note: this is non-linear since it imply some 1/cos() [2].


Edit: I'am tempted to think we should first reproject each point using d3.geo.mercator() on the screen (HOW? Which syntax ?), so D3 do all the hard maths. We could then GET the point's pixels coordinates, so we can calculate the length AB and the length AC in pixels, and finally r=AC/AB. Also, it's a bit how to convert decimal degrees coordinates into projected pixels coordinates function of the chosen d3.geo.<PROJECTIONNAME>() ?

[2]: Mercator: scale factor is changed along the meridians as a function of latitude?


Solution

  • I will assume that the points are A: (0, 50), B: (0, 40), C: (10, 50) and D: (10, 40). The feature enclosed by the points (A, C, D, B) will look as a square using the equirectangular projection. Now, the points are longitude, latitude pairs, you can compute the great-arc distance between the points using d3.geo.distance. This will give you the angular distance between the points. For instance:

    // Points (lon, lat)
    var A = [ 0, 50],
        B = [ 0, 40],
        C = [10, 50],
        D = [10, 40];
    
    // Geographic distance between AB and AC
    var distAB = d3.geo.distance(A, B),  // 0.17453292519943306 radians
        distAC = d3.geo.distance(A, C);  // 0.11210395570214343 radians
    

    Now, these distances are the angles between the points, as you can see, the feature wasn't a square. If we project the points using the D3 Mercator projection:

    // The map will fit in 800 px horizontally
    var width = 800;
    var mercator = d3.geo.mercator()
        .scale(width / (2 * Math.PI));
    
    // Project (lon, lat) points using the projection, to get pixel coordinates.
    var pA = mercator(A),  // [480, 121] (rounded)
        pB = mercator(B),  // [480, 152] (rounded)
        pC = mercator(C);  // [502, 121] (rounded)
    

    And now use the euclidean distance to compute the distance between the projected points pA, pB and pC.

    function dist(p, q) {
        return Math.sqrt(Math.pow(p[0] - q[0], 2) + Math.pow(p[1] - q[1], 2));
    }
    
    var pDistAB = dist(pA, pB),  // 31.54750649588999 pixels
        pDistAC = dist(pA, pC);  // 22.22222222222223 pixels
    

    If you use angular distances as reference, you will get two ratios, one for AB and other for AC:

    var ratioAB = distAB / pDistAB, // 0.005532384159178197 radians/pixels
        ratioAC = distAC / pDistAC; // 0.005044678006596453 radians/pixels
    

    If you use the equirectangular projection as reference, you can use the euclidean distance between the points (as if they were in a plane surface):

    var ratioAB = dist(A, B) / pDistAB, // 0.3169822629659431  degrees/pixels
        ratioAC = dist(A, C) / pDistAC; // 0.44999999999999984 degrees/pixels