Skip to content

Instantly share code, notes, and snippets.

@mpetroff

mpetroff/index.html

Created October 30, 2021 20:41
Show Gist options
  • Save mpetroff/dcf7b090eabda85d081d47e8a0c71d4a to your computer and use it in GitHub Desktop.
Save mpetroff/dcf7b090eabda85d081d47e8a0c71d4a to your computer and use it in GitHub Desktop.
Download ZIP
Implementation of "A Square Equal-Area Map Projection with Low Angular Distortion, Minimal Cusps, and Closed-Form Solutions": https://doi.org/10.1145/3460521
Raw
index.html
<!DOCTYPE html>
<html lang="en">
<meta charset="utf-8">
<title>New Projection</title>
<style>
.stroke {
fill: none;
stroke: #000;
stroke-width: 2px;
}
.fill {
fill: #fff;
}
.graticule {
fill: none;
stroke: #666;
stroke-width: 2px;
mix-blend-mode: exclusion;
}
.land {
fill: #222;
}
</style>
<script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/d3.min.js" integrity="sha256-Xb6SSzhH3wEPC4Vy3W70Lqh9Y3Du/3KxPqI2JHQSpTw=" crossorigin="anonymous"></script>
<script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/d3-geo-projection.min.js" integrity="sha256-LY1qVgxfyQnPl8CVKGRKGIfi5BUjq9aMVgxEWC26gT0=" crossorigin="anonymous"></script>
<script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/topojson.min.js" integrity="sha256-tHoAPGoNdhIR28YHl9DWLzeRfdwigkH7OCBXMrHXhoM=" crossorigin="anonymous"></script>
<script src="newprojection.js"></script>
<svg width="678" height="678" id="small"></svg>
<script>
function projectionRaw(lambda, phi) {
const sLambda = Math.sign(lambda),
sPhi = Math.sign(phi),
cosPhi = Math.cos(phi),
x = Math.cos(lambda) * cosPhi,
y = Math.sin(lambda) * cosPhi;
let z = Math.sin(sPhi * phi);
lambda = Math.abs(Math.atan2(y, z));
phi = Math.asin(x);
if (Math.abs(lambda - Math.PI / 2) > 1e-6) lambda %= Math.PI / 2;
const point = hexadecant(lambda > Math.PI / 4 ? Math.PI / 2 - lambda : lambda, phi);
if (lambda > Math.PI / 4) z = point[0], point[0] = -point[1], point[1] = -z;
return (point[0] *= sLambda, point[1] *= -sPhi, point);
};
function projectionRawInvert(x, y) {
if (Math.abs(x) > 1)
x = Math.sign(x) * 2 - x;
if (Math.abs(y) > 1)
y = Math.sign(y) * 2 - y;
const sx = Math.sign(x),
sy = Math.sign(y),
x0 = -sx * x,
y0 = -sy * y,
t = y0 / x0 < 1,
p = hexadecantInvert(t ? y0 : x0, t ? x0 : y0),
lambda = t ? -Math.PI / 2 - p[0] : p[0],
phi = p[1],
cosPhi = Math.cos(phi);
return [sx * (Math.atan2(Math.sin(lambda) * cosPhi, -Math.sin(phi)) + Math.PI), sy * Math.asin(Math.cos(lambda) * cosPhi)];
};
function projectionRaw2(lambda, phi) {
// Go through inverse to make sure it works
var tmp = projectionRaw(lambda, phi);
tmp = projectionRawInvert(tmp[0], tmp[1]);
return projectionRaw(tmp[0], tmp[1]);
}
projection = d3.geoQuincuncial(projectionRaw2).scale(176.423);
var svg = d3.select("#small"),
width = +svg.attr("width"),
height = +svg.attr("height");
var projection = projection
.scale(239)
.translate([width / 2, height / 2])
.precision(0.1);
var path = d3.geoPath()
.projection(projection);
var graticule = d3.geoGraticule();
svg.append("defs").append("path")
.datum({type: "Sphere"})
.attr("id", "sphere")
.attr("d", path);
svg.append("use")
.attr("class", "fill")
.attr("xlink:href", "#sphere");
svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("d", path);
svg.append("use")
.attr("class", "stroke")
.attr("xlink:href", "#sphere");
d3.json("https://cdn.jsdelivr.net/npm/[email protected]/land-110m.json").then(world => {
svg.insert("path", ".graticule")
.datum(topojson.feature(world, world.objects.land))
.attr("class", "land")
.attr("d", path);
});
</script>
Raw
newprojection.js
/*
* Matthew Petroff, October 2021
*
* Implementation of "A Square Equal-Area Map Projection with Low Angular
* Distortion, Minimal Cusps, and Closed-Form Solutions":
* https://doi.org/10.1145/3460521
*
* This work is placed into the public domain via the CC0 1.0 public domain
* dedication: https://creativecommons.org/publicdomain/zero/1.0/
*
*/
function hexadecant(lambda, phi) {
// Handle edge cases
if (phi === Math.PI / 2) return [0, 0];
if (lambda > Math.PI / 4 - 1e-12)
lambda = Math.PI / 4 - 1e-12;
const phi0 = 3 * Math.PI / 8;
const cos_phi0 = Math.cos(phi0);
const sin_phi0 = Math.sin(phi0);
const psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
const psi1 = Math.PI - 2 * psi0;
const rho = Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0)));
const hprime = 12 / Math.PI * (psi0 + rho - Math.PI / 2);
const xiprime = Math.atan((Math.PI * (hprime - 3) * (hprime - 3)) / (Math.sqrt(3) *
(Math.PI * (hprime * hprime - 2 * hprime + 45) - 96 * psi0 - 48 * rho)));
const psi0prime = Math.atan(Math.sqrt(3) / hprime);
const psi1prime = 7 * Math.PI / 6 - psi0prime - xiprime;
const psi2prime = xiprime - Math.PI / 6;
const rhoprime = Math.atan(hprime / Math.sqrt(3));
const lambda0 = Math.floor((lambda + Math.PI / 4) / (Math.PI / 2)) * Math.PI / 2;
const theta = Math.abs(Math.atan2(Math.cos(phi) * Math.sin(lambda - lambda0), sin_phi0 *
Math.cos(phi) * Math.cos(lambda - lambda0) - cos_phi0 * Math.sin(phi)));
const r = Math.acos(sin_phi0 * Math.sin(phi) + cos_phi0 * Math.cos(phi) *
Math.cos(lambda - lambda0));
let beta;
if (theta <= psi0)
beta = psi0 - theta;
else if (theta <= psi0 + psi1)
beta = theta - psi0;
else
beta = Math.PI - theta;
const c = theta <= psi0 + psi1 ? Math.acos(cos_phi0 / Math.SQRT2) : Math.PI / 2 - phi0;
let G, Gprime, F;
if (theta <= psi0) {
G = psi0;
Gprime = psi0prime;
F = rho;
} else if (theta <= psi0 + psi1) {
G = psi1;
Gprime = psi1prime;
F = Math.PI / 2 - rho;
} else {
G = psi0;
Gprime = psi2prime;
F = Math.PI / 4;
}
const aprime = theta <= psi0 ? hprime : Math.sqrt(hprime * hprime + 3) *
Math.sin(Math.PI / 3 - rhoprime) / Math.sin(xiprime);
const cprime = theta <= psi0 + psi1 ? Math.sqrt(hprime * hprime + 3) : 3 - hprime;
const x = Math.acos(Math.cos(r) * Math.cos(c) + Math.sin(r) * Math.sin(c) * Math.cos(beta));
const gamma = Math.asin(Math.sin(beta) * Math.sin(r) / Math.sin(x));
const epsilon = Math.acos(Math.sin(G) * Math.sin(gamma) * Math.cos(c) - Math.cos(G) * Math.cos(gamma));
const upupvp = (gamma + G + epsilon - Math.PI) / (F + G - Math.PI / 2);
const cos_xy = Math.sqrt(1 - Math.pow(Math.sin(G) * Math.sin(c) / Math.sin(epsilon), 2));
const xpxpyp = Math.sqrt((1 - Math.cos(x)) / (1 - cos_xy));
const uprime = aprime * upupvp;
const xpyp = Math.sqrt(uprime * uprime + cprime * cprime - 2 * uprime * cprime * Math.cos(Gprime));
const cos_gammaprime = Math.sqrt(1 - Math.pow((uprime * Math.sin(Gprime) / xpyp), 2));
const xprime = xpyp * xpxpyp;
const yprime = xpyp - xprime;
let rprime = Math.sqrt(xprime * xprime + cprime * cprime - 2 * xprime * cprime * cos_gammaprime);
const alphaprime = Math.acos((yprime * yprime - uprime * uprime - rprime * rprime) / (-2 * uprime * rprime));
// Put theta back in the correct section
let thetaprime;
if (theta <= psi0)
thetaprime = alphaprime;
else if (theta <= psi0 + psi1)
thetaprime = 7 * Math.PI / 6 - xiprime - alphaprime;
else
thetaprime = 7 * Math.PI / 6 - xiprime + alphaprime;
// Handle edge cases
if (lambda < 1e-12)
thetaprime = phi < phi0 ? 0 : Math.PI;
if (x === 0) {
thetaprime = Gprime;
rprime = cprime;
}
const xm = rprime * Math.sin(thetaprime);
const ym = -rprime * Math.cos(thetaprime) - (3 - hprime);
const xy0 = xm * Math.sqrt(3) / 3;
const xy1 = ym / 3;
return [xy0, xy1];
}
function hexadecantInvert(x_m, y_m) {
const phi0 = 3 * Math.PI / 8,
cos_phi0 = Math.cos(phi0),
sin_phi0 = Math.sin(phi0);
const x_c = 3 * x_m / Math.sqrt(3),
y_h = 3 * y_m;
const hprime = 12 / Math.PI * (Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0)) +
Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0))) - Math.PI / 2);
const rprime = Math.sqrt(x_c * x_c + (hprime - 3 - y_h) * (hprime - 3 - y_h));
const thetaprime = Math.abs(Math.atan2(x_c, hprime - 3 - y_h));
const xiprime = Math.atan((Math.PI * (hprime - 3) * (hprime - 3)) / (Math.sqrt(3) *
(Math.PI * (hprime * hprime - 2 * hprime + 45) -
96 * Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0)) -
48 * Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0))))));
const psi0prime = Math.atan(Math.sqrt(3) / hprime);
const psi1prime = 7 * Math.PI / 6 - psi0prime - xiprime;
const psi2prime = xiprime - Math.PI / 6;
let alphaprime;
if (thetaprime <= psi0prime)
alphaprime = thetaprime;
else if (thetaprime <= psi0prime + psi1prime)
alphaprime = Math.PI - psi2prime - thetaprime;
else
alphaprime = thetaprime + psi2prime - Math.PI;
let c = thetaprime <= psi0prime + psi1prime ? Math.acos(cos_phi0 / Math.SQRT2) : Math.PI / 2 - phi0;
const psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
const psi1 = Math.PI - 2 * psi0;
const rho = Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0)));
// The original publication neglected to mention that the non-prime to
// prime replacement should only be done in the condition statements, so G
// is still set to psi0 / psi1.
let G, Gprime, F;
if (thetaprime <= psi0prime) {
G = psi0;
Gprime = psi0prime;
F = rho;
} else if (thetaprime <= psi0prime + psi1prime) {
G = psi1;
Gprime = psi1prime;
F = Math.PI / 2 - rho;
} else {
G = psi0;
Gprime = psi2prime;
F = Math.PI / 4;
}
const aprime = thetaprime <= psi0prime ? hprime : Math.sqrt(hprime * hprime + 3) *
Math.sin(Math.PI / 3 - Math.atan(hprime / Math.sqrt(3))) / Math.sin(xiprime);
const cprime = thetaprime <= psi0prime + psi1prime ? Math.sqrt(hprime * hprime + 3) : 3 - hprime;
let b;
if (thetaprime <= psi0prime)
b = Math.PI / 4;
else if (thetaprime <= psi0prime + psi1prime)
b = Math.atan(Math.SQRT2 * Math.tan(phi0));
else
b = Math.PI / 2 - Math.atan(Math.SQRT2 * Math.tan(phi0));
const betaprime = Gprime - alphaprime;
const xprime = Math.sqrt(rprime * rprime + cprime * cprime - 2 * rprime *
cprime * Math.cos(betaprime));
const gammaprime = xprime > 0 ? Math.asin(rprime * Math.sin(betaprime) / xprime) : 0; // Avoid divide by zero
const epsilonprime = Math.PI - Gprime - gammaprime;
const yprime = epsilonprime > 0 ? rprime * Math.sin(alphaprime) / Math.sin(epsilonprime) : 0; // Avoid divide by zero
const uprime = Math.sqrt(Math.abs(cprime * cprime +
(xprime + yprime) * (xprime + yprime) -
2 * cprime * (xprime + yprime) * Math.cos(gammaprime))); // Take absolute value to avoid negative numbers due numerical precision issues
const vprime = aprime - uprime;
const delta = Math.atan(
(-Math.sin(vprime * (F + G - Math.PI / 2) / aprime))
/ (Math.cos(b) - Math.cos(vprime * (F + G - Math.PI / 2) / aprime))
);
const gamma = F - delta;
const cos_xy = 1 / Math.sqrt(1 + (Math.tan(b) / Math.cos(delta)) * (Math.tan(b) / Math.cos(delta)));
const x = Math.acos(1 - (xprime / (xprime + yprime)) * (xprime / (xprime + yprime)) * (1 - cos_xy));
const r = Math.acos(Math.cos(x) * Math.cos(c) + Math.sin(x) * Math.sin(c) * Math.cos(gamma));
const beta = Math.asin(Math.sin(x) * Math.sin(gamma) / Math.sin(r));
let alpha;
if (thetaprime <= psi0prime)
alpha = psi0 - beta;
else if (thetaprime <= psi0prime + psi1prime)
alpha = beta + psi0; // There was a sign error here in the original publication
else
alpha = Math.PI - beta;
const phi = Math.asin(sin_phi0 * Math.cos(r) - cos_phi0 * Math.sin(r) * Math.cos(alpha));
let lambda = Math.sign(x_c) * Math.atan(Math.sin(alpha) *
Math.sin(r) * cos_phi0 / (Math.cos(r) - sin_phi0 * Math.sin(phi)));
lambda *= Math.abs(phi) != Math.PI / 2 // Handle edge case
return [lambda, phi];
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment