Lee’s Tetrahedral Conformal Projection, South by Fil

index.html

<!DOCTYPE html>
<!-- 578 ~= 500 * 2 / sqrt(3) -->
<div style="text-align:center"><canvas width="578" height="500"></canvas></div>
<script src="https://d3js.org/d3.v4.js"></script>
<!-- <script src="https://d3js.org/d3-geo-projection.v2.js"></script> -->
<script src="d3-geo.js"></script>
<script src="d3-geo-projection.js"></script>
<script src="https://d3js.org/topojson.v2.min.js"></script>
<script src="https://unpkg.com/complex.js"></script>
  

<script>

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

  // retina display
  var devicePixelRatio = window.devicePixelRatio || 1;
  canvas.style('width', canvas.attr('width')+'px');
  canvas.style('height', canvas.attr('height')+'px');
  canvas.attr('width', canvas.attr('width') * devicePixelRatio);
  canvas.attr('height', canvas.attr('height') * devicePixelRatio);
  context.scale(devicePixelRatio,devicePixelRatio);


  
  
// import for d3.geo
var geoStream = d3.geoStream;

// import from math
var epsilon = 1e-6, epsilon2 = epsilon * epsilon, asin = Math.asin;
var pi = Math.PI, degrees = 180 / pi, asin1_3 = Math.asin(1 / 3);



var leeRaw = function(lambda, phi) {
  // return d3.geoGnomonicRaw(...arguments);
  var w = Complex([-1/2, Math.sqrt(3)/2]),
      k = Complex(0),
      h = Complex(0),
      z = Complex(d3.geoStereographicRaw(lambda, phi)).mul(Math.sqrt(2));

  // rotate to have s ~= 1
  var rot = w.clone().pow(d3.scan([0,1,2].map(
    i => -(z.clone().mul(w.clone().pow(i))).re
  )));

  var n = z.abs();

  if (n > 0.3) {


    // if |z| > 0.5, use the approx based on y = (1-z)
    // McIlroy formula 6 p6 and table for G page 16
    var y = rot.clone().mul(z).mul(-1).add(1);

    // w1 = gamma(1/3) * gamma(1/2) / 3 / gamma(5/6);
    // https://bl.ocks.org/Fil/1aeff1cfda7188e9fbf037d8e466c95c
    var w1 = 1.4021821053254548;

    var G0 = [
      1.15470053837925,
      0.192450089729875,
      0.0481125224324687,
      0.010309826235529,
      3.34114739114366e-4,
      -1.50351632601465e-3,
      -1.23044177962310e-3,
      -6.75190201960282e-4,
      -2.84084537293856e-4,
      -8.21205120500051e-5,
      -1.59257630018706e-6,
      1.91691805888369e-5,
      1.73095888028726e-5,
      1.03865580818367e-5,
      4.70614523937179e-6,
      1.4413500104181e-6,
      1.92757960170179e-8,
      -3.82869799649063e-7,
      -3.57526015225576e-7,
      -2.2175964844211e-7
    ];

    var G = Complex(0);
    for (var i = G0.length; i--;) {
      G = Complex(G0[i]).add(G.mul(y));
    }

    k = Complex(w1).add(y.sqrt().mul(-1).mul(G)).mul(rot).mul(rot)

  }

  if (n < 0.5) {

    // if |z| < 0.3
    // https://www.wolframalpha.com/input/?i=series+of+((1-z%5E3))+%5E+(-1%2F2)+at+z%3D0 (and ask for "more terms")
    // 1 + z^3/2 + (3 z^6)/8 + (5 z^9)/16 + (35 z^12)/128 + (63 z^15)/256 + (231 z^18)/1024 + O(z^21)
    // https://www.wolframalpha.com/input/?i=integral+of+1+%2B+z%5E3%2F2+%2B+(3+z%5E6)%2F8+%2B+(5+z%5E9)%2F16+%2B+(35+z%5E12)%2F128+%2B+(63+z%5E15)%2F256+%2B+(231+z%5E18)%2F1024
    // (231 z^19)/19456 + (63 z^16)/4096 + (35 z^13)/1664 + z^10/32 + (3 z^7)/56 + z^4/8 + z + constant
    var H0 = [
      1, 1/8, 3/56, 1/32, 35/1664, 63/4096, 231/19456
    ]
    var z3 = z.clone().pow(3);
    for (var i = H0.length; i--;) {
      h = Complex(H0[i]).add(h.mul(z3));
    }
    h = h.mul(z);
  }


  if (n < 0.3) return h.toVector();
  if (n > 0.5) return k.toVector();

  // in between 0.3 and 0.5, interpolate
  var t = (n - 0.3) / (0.5 - 0.3);
  return k.mul(t).add(h.mul(1 - t)).toVector();

}
  

var centers = [
  [0, 90],
  [-180, -asin1_3 * degrees],
  [-60, -asin1_3 * degrees],
  [60, -asin1_3 * degrees]
];
  var tetrahedron = [[1, 2, 3], [0, 2, 1], [0, 3, 2], [0, 1, 3]].map(function(
    face
  ) {
    return face.map(function(i) {
      return centers[i];
    });
  });

 
d3.geoTetrahedralLee = function(faceProjection) {

  faceProjection =
    faceProjection ||
    function(face) {
      var c = d3.geoCentroid({ type: "MultiPoint", coordinates: face }),
          rotate = [ -c[0], -c[1], 30 ];
      if (Math.abs(c[1]) == 90) {
        rotate = [ 0, -c[1], -30 ];
      }
      return d3
        .geoProjection(leeRaw)
        .scale(1)
        .translate([0, 0])
        .rotate(rotate);
    };

  var faces = tetrahedron.map(function(face) {
    return { face: face, project: faceProjection(face) };
  });

  [-1, 0, 0, 0].forEach(function(d, i) {
    var node = faces[d];
    node && (node.children || (node.children = [])).push(faces[i]);
  });

  return d3
    .geoPolyhedral(
      faces[0],
      function(lambda, phi) {
        lambda *= degrees;
        phi *= degrees;
        for (var i = 0; i < faces.length; i++) {
          if (
            d3.geoContains(
              {
                type: "Polygon",
                coordinates: [[...tetrahedron[i], tetrahedron[i][0]]]
              },
              [lambda, phi]
            )
          ) {
            return faces[i];
          }
        }
      },
      pi / 6
    )
    .clipAngle(180) // this is only to avoid antimeridian clipping on the Sphere
    .precision(0.05)
    //.rotate([-30, 0])
    .rotate([30, 180]) // for North Pole aspect, needs clipPolygon
    .fitExtent([[1, 1], [width-1, height-1]], {type: "Sphere"})
  
};

projection = d3.geoTetrahedralLee();

  
// OR, from the `complex+clipPolygon`  branch of d3-geo-projection:
if (0) projection = d3.geoPolyhedralLee(30)
  .rotate([30, 180]) // for North Pole aspect, needs clipPolygon
  .fitExtent([[1, 1], [width-1, height-1]], {type: "Sphere"});

  
  

  
  
var init_scale = projection.scale(),
  path = d3.geoPath().projection(projection).context(context);

d3.json("https://unpkg.com/world-atlas@1/world/110m.json", function(
  error,
  world
) {
  if (error) throw error;

  var land = topojson.merge(world, world.objects.countries.geometries);

  render = function() {
    var tiling = false;
    context.fillStyle = "#fff";
    context.fillRect(0, 0, width, height);

    if (!tiling) {

    context.beginPath();
    path(d3.geoGraticule()());
    context.strokeStyle = "#777";
    context.lineWidth = 0.5;
    context.stroke(), context.closePath();

    // equator
    context.beginPath();
    path(d3.geoCircle().center([0,90]).radius(90)());
    context.strokeStyle = "#000";
    context.lineWidth = 1;
    context.stroke(), context.closePath();

    // inner triangle
    context.beginPath();
    {
      let rotate = projection.rotate();
      var inner = centers.map(projection.rotate([0,0]));
      projection.rotate(rotate);
    }
      context.moveTo(inner[1][0], inner[1][1]);
      context.lineTo(inner[2][0], inner[2][1]);
      context.lineTo(inner[3][0], inner[3][1]);
      context.lineTo(inner[1][0], inner[1][1]);
    context.strokeStyle = "#777";
    context.lineWidth = 0.5;
    context.setLineDash([5,3]);
    context.stroke(), context.closePath();
    context.setLineDash([]);

    }

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

    context.beginPath();
    var now = performance.now();
    path(land);
    console.log('time', Math.round(performance.now()-now)+'ms');
    context.lineWidth = 1;
    context.strokeStyle = "#000";
    context.stroke();
    context.fillStyle = "#000";
    context.fill();
    context.closePath();

    console.log(projection([0,-90]))
    if (tiling) {
    context.beginPath();
    context.rotate(pi);
    context.translate(-1247,-500);
    path(land);
    context.translate(575,0);
    path(land);
    context.lineWidth = 1;
    //context.strokeStyle = "red";
    context.stroke();
    //context.fillStyle = "pink";
    context.fill();
    context.closePath();
    }
  };

  render();
});


</script>