Skip to content

Instantly share code, notes, and snippets.

@Fil

Fil/.block Secret

Last active July 25, 2017 05:26
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Fil/5ad65a8075e634ab7e4de8ace95fb3c5 to your computer and use it in GitHub Desktop.
Save Fil/5ad65a8075e634ab7e4de8ace95fb3c5 to your computer and use it in GitHub Desktop.
Lee’s Tetrahedric Conformal Projection (work in progress)
license: mit

The goal is to reproduce this map LeeTetrahedral @ progonos with d3, see d3-geo-projection issue #107.

Reference paper Laurence Patrick Lee, Cartographica: The International Journal for Geographic Information and Geovisualization > Vol. 13, # 1, June 1976 (100+ pages!)

See North Polar aspect.


[

](https://github.com/Fil/) Questions and comments welcome on [gitter.im/d3](https://gitter.im/d3/d3), [twitter](https://twitter.com/@recifs) or [slack](https://d3js.slack.com). <script> (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){ (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o), m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m) })(window,document,'script','https://www.google-analytics.com/analytics.js','ga'); ga('create', 'UA-58621-8', 'auto'); ga('send', 'pageview'); </script>
<!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.min.js"></script>
<script src="https://d3js.org/topojson.v2.min.js"></script>
<script src="https://unpkg.com/complex.js"></script>
<script>
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 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);
var pi = Math.PI, degrees = 180 / pi, asin1_3 = Math.asin(1 / 3);
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) // should be clipPolygon!
.rotate([-30, 0])
.precision(0.05)
//.rotate([30, 180]) // for North Pole aspect, needs clipPolygon
.fitExtent([[1, 1], [width-1, height-1]], { type: "Sphere" });
};
projection = d3.geoTetrahedralLee();
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({type:"Sphere"});
context.strokeStyle = "black";
context.lineWidth = 1.5;
context.stroke(), context.closePath();
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.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>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment