Skip to content

Instantly share code, notes, and snippets.

@Fil
Last active May 3, 2018 13:33
Show Gist options
  • Save Fil/ce9d55cc070126b66fb86bc03b84b3f4 to your computer and use it in GitHub Desktop.
Save Fil/ce9d55cc070126b66fb86bc03b84b3f4 to your computer and use it in GitHub Desktop.
Cox Projection (North polar aspect)
license: mit

Conformal projection of the Sphere within an equilateral triangle.

This is the North polar aspect, see the original Cox Projection (South polar aspect).

Original research by Philippe Rivière for d3-geo-projection, issue #118.

Read a write-up at visionscarto.net.

Reference papers:

  • J. F. Cox, Représentation de la surface entière de la terre dans un triangle équilatéral. Bulletin de la Classe des Sciences, Académie Royale de Belgique, 5, 21: 66-71. Bruxelles, 1935.

  • J. Magis, Calcul du canevas de la représentation conforme de la sphère entière dans un triangle équilatéral. Bulletin Geodésique, 247-256. Paris, 1938.

  • L.P. Lee, Conformal Projections Based On Dixon Elliptic Functions, Cartographica: The International Journal for Geographic Information and Geovisualization > Vol. 13, # 1, June 1976.

  • M. Douglas McIlroy, Wallpaper maps, 2011.


[

](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>
<canvas width="960" height="500"></canvas>
<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>
// Approximate \int _0 ^sm(z) dt / (1 - t^3)^(2/3)
// sm maps a triangle to a disc, sm^-1 does the opposite
function sm_1(s){
var z = Complex(s),
w = Complex([-1/2, Math.sqrt(3)/2]),
k = Complex(0);
// 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 y = rot.clone().mul(z).mul(-1).add(1);
// var n = 3
// w1 = gamma(1/n) * gamma(1 - 2/n) / n / gamma(1 - 1/n)
// https://bl.ocks.org/Fil/852557838117687bbd985e4b38ff77d4
var w1 = 1.7666387502854533;
// McIlroy formula 5 p6 and table for F3 page 16
var F0 = [
1.44224957030741,
0.240374928384568,
0.0686785509670194,
0.0178055502507087,
0.00228276285265497,
-1.48379585422573e-3,
-1.64287728109203e-3,
-1.02583417082273e-3,
-4.83607537673571e-4,
-1.67030822094781e-4,
-2.45024395166263e-5,
2.14092375450951e-5,
2.55897270486771e-5,
1.73086854400834e-5,
8.72756299984649e-6,
3.18304486798473e-6,
4.79323894565283e-7
-4.58968389565456e-7,
-5.62970586787826e-7,
-3.92135372833465e-7
];
var F = Complex(0);
for (var i = F0.length; i--;) {
F = Complex(F0[i]).add(F.mul(y));
}
// /!\ mutates y
var k = Complex(w1).add(y.pow(1-2/3).mul(-1).mul(F)).mul(rot.pow(2))
// when we are close to [0,0] we switch to another approximation:
// https://www.wolframalpha.com/input/?i=(-2%2F3+choose+k)++*+(-1)%5Ek++%2F+(k%2B1)+with+k%3D0,1,2,3,4
// the difference is _very_ tiny but necessary
// if we want projection(0,0) === [0,0]
var n = z.abs(), m = 0.3;
if (n < m) {
var H0 = [
1, 1/3, 5/27, 10/81, 22/243 //…
];
var z3 = z.clone().pow(3);
var h = Complex(0);
for (i = H0.length; i--;) {
h = Complex(H0[i]).add(h.mul(z3));
}
h = h.mul(z);
k.mul(n / m).add(h.mul(1 - n / m));
}
return k.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);
function coxRawN(lambda, phi){
var s = d3.geoLagrangeRaw(0.5)(lambda, phi);
var t = sm_1([-s[1] / 2, -s[0] / 2]);
return [-t[1], -t[0]]
}
projection = d3.geoProjection(coxRawN)
.rotate([-10.8,0])
.scale(188.68216)
.translate([480, 500 * 1/3])
.precision(0.01);
// the Sphere should go *exactly* to the vertices of the triangles
// because they're singular points
function sphere() {
var degrees = 180 / Math.PI,
c = 2 * Math.asin(1/Math.sqrt(5))*degrees;
return {
type: "Polygon",
coordinates: [
[ [ 0,-90 ], [ -180, c ], [ 0, 90 ], [ 180, c ], [ 0,-90 ] ]
]
};
}
var stream_ = projection.stream;
projection.stream = function(stream) {
var rotate = projection.rotate(),
rotateStream = stream_(stream),
sphereStream = (projection.rotate([0, 0]), stream_(stream));
projection.rotate(rotate);
rotateStream.sphere = function() { d3.geoStream(sphere(), sphereStream); };
return rotateStream;
};
var 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() {
context.fillStyle = "#fff";
context.fillRect(0, 0, width, height);
context.beginPath();
path({type:"Sphere"});
context.strokeStyle = "black";
context.lineWidth = 1.5;
context.stroke(), context.closePath();
context.beginPath();
path(land);
context.lineWidth = 1;
context.strokeStyle = "#000";
context.stroke();
context.fillStyle = "#000";
context.fill();
context.closePath();
context.beginPath();
path(d3.geoGraticule()());
context.strokeStyle = "#777";
context.lineWidth = 0.5;
context.stroke(), context.closePath();
};
render();
});
</script>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment