31:T414,md`# Adams World in a Square I + II The two projections of the *world in a square* presented here were developed by Oscar Sherman Adams (1874-1962) in 1929. The projections are *conformal*, except at the poles and at two points on each of the 180th meridians / at the intersections of the 180th meridians with the Equator. Together with the [Peirce Quincuncial](https://observablehq.com/@d3/peirce-quincuncial) and the [Guyou projection](https://observablehq.com/@d3/guyou) they belong to a family of projections, that require to compute an [elliptic integral](https://en.wikipedia.org/wiki/Elliptic_integral) of the first kind. For this task I've implemented an efficiently converging algorithm from Roland Bulirsch (as suggested by John P. Snyder) besides the forward and inverse of the two projections. The *Adams World in a Square II* projection (without the 45° rotation in the plane) serves as the basis for the [Spilhaus World Ocean Map in a Square](https://observablehq.com/@toja/spilhaus-world-ocean-map-in-a-square) projection. `32:T720,{ const w = width / 4, h = height / 4; const map = DOM.context2d(w, h); const proj = adamsSquareII().fitSize([w, h], sphere); const path = d3.geoPath(proj, map); map.save(); map.beginPath(); map.moveTo(w / 2, 0), map.lineTo(w, h / 2), map.lineTo(w / 2, h), map.lineTo(0, h / 2), map.closePath(); map.clip() map.fillStyle = "#faf8f6"; map.fillRect(0, 0, w, h); map.restore(); map.lineWidth = 0.25; map.strokeStyle = "#234"; map.fillStyle = "#234"; map.beginPath(); path(graticule); map.stroke(); map.beginPath(); path(sphere); map.stroke(); map.beginPath(); path(land); map.fill(); const image = map.canvas; // Don't look at the code! const context = DOM.context2d(width, height / 2); context.translate(-w / 2, -h / 2); context.save(); context.scale(1, -1); for (let i = 0; i < 5; ++i) { context.drawImage(image, w * i, -h * (i & 1 ? 2 : 1), w, h); context.drawImage(image, w * i, -h * (i & 1 ? 2 : 1) * 3, w, h); } context.restore(); context.save(); context.translate(0, h); for (let i = 0; i < 5; ++i) { context.drawImage(image, w * i, h * (i & 1 ? 0 : 1), w, -h); context.drawImage(image, w * i, h + (i & 1 ? 0 : h), w, h); } context.restore(); context.save(); context.translate(width / 2, height / 4); context.rotate(Math.PI / 2); context.translate(w / 2, h / 2); for (let i = 0; i < 5; ++i) { context.drawImage(image, -w * (i & 1 ? 1 : 0), -h * i, w, h); } context.restore(); context.save(); context.translate(width / 2, height / 4); context.rotate(-Math.PI / 2); context.translate(-w * 1.5, -h * 2.5); for (let i = 0; i < 5; ++i) { context.drawImage(image, w * (i & 1 ? 1 : 0), h * i, w, h); } context.restore(); return context.canvas; }33:T532,/* * Computes the elliptic integral of the first kind. * Algorithm from Bulirsch(1965), the implementation follows Snyder(1989), p. 239. * A faster alternative for m = 0.5 is presented in: * Gerald I. Evenden (2008), libproj4: A Comprehensive Library of * Cartographic Projection Functions (Preliminary Draft), p. 123. */ ellipticF = function(phi, m) { const { abs, atan, ln, PI: pi, sin, sqrt } = Math; const C1 = 10e-4, C2 = 10e-10, TOL = 10e-6; const sp = sin(phi); let k = sqrt(1 - m), h = sp * sp; // "complete" elliptic integral if (h >= 1 || abs(phi) === pi / 2) { if (k <= TOL) return sp < 0 ? -Infinity : Infinity; m = 1, h = m, m += k; while(abs(h - k) > C1 * m) { k = sqrt(h * k); m /= 2, h = m, m += k; } return sp < 0 ? -pi / m : pi / m; } // "incomplete" elliptic integral else { if (k <= TOL) return ln((1 + sp) / (1 - sp)) / 2; let g, n, p, r, y; m = 1, n = 0, g = m, p = m * k, m += k; y = sqrt((1 - h) / h); if (abs(y -= p / y) <= 0) y = C2 * sqrt(p); while (abs(g - k) > C1 * g) { k = 2 * sqrt(p), n += n; if (y < 0) n += 1; p = m * k, g = m, m += k; if (abs(y -= p / y) <= 0) y = C2 * sqrt(p); } if (y < 0) n += 1; r = (atan(m / y) + pi * n) / m; return sp < 0 ? -r : r; } }34:T6ce,/* * Newton-Raphson inversion, based on code from PROJ written by Gerald Evenden. * https://github.com/OSGeo/PROJ/blob/master/src/projections/adams.cpp */ function inverse(x, y, lam, phi, proj) { const { abs, max, min, PI } = Math; const HALFPI = Math.PI / 2; let dLamX = 0, dLamY = 0, dPhiX = 0, dPhiY = 0; let x2, y2, lam2, phi2; for(let i = 0; i < 15; ++i) { const [xAppr, yAppr] = proj(lam, phi); const dX = xAppr - x; const dY = yAppr - y; if(abs(dX) < 1e-10 && abs(dY) < 1e-10) { return [lam, phi]; } if(abs(dX) > 1e-6 || abs(dY) > 1e-6) { // Compute Jacobian matrix const dLam = lam > 0 ? -1e-6 : 1e-6; lam2 = lam + dLam; phi2 = phi; [x2, y2] = proj(lam2, phi2); const dXLam = (x2 - xAppr) / dLam; const dYLam = (y2 - yAppr) / dLam; const dPhi = phi > 0 ? -1e-6 : 1e-6; lam2 = lam; phi2 = phi + dPhi; [x2, y2] = proj(lam2, phi2); const dXPhi = (x2 - xAppr) / dPhi; const dYPhi = (y2 - yAppr) / dPhi; // Inverse of Jacobian matrix const det = dXLam * dYPhi - dXPhi * dYLam; if(det != 0) { dLamX = dYPhi / det; dLamY = -dXPhi / det; dPhiX = -dYLam / det; dPhiY = dXLam / det; } } // Limit the amplitude of correction if(x != 0) { const dLam = max(min(dX * dLamX + dY * dLamY, 0.3), -0.3); lam -= dLam; if(lam < -PI) lam = -PI; else if(lam > PI) lam = PI; } if(y != 0) { const dPhi = max(min(dX * dPhiX + dY * dPhiY, 0.3), -0.3); phi -= dPhi; if(phi < -HALFPI) phi = -HALFPI; else if(phi > HALFPI) phi = HALFPI; } } return [lam, phi]; }35:[] 2f:{"title":"Adams World in a Square I + II","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":997,"value":"{\n const context = DOM.context2d(width, height);\n const proj = adamsSquareI().fitExtent([[50 , 50], [width - 50, height - 50]], sphere).rotate([-10,0])\n const path = d3.geoPath(proj, context);\n \n context.fillStyle = \"#faf8f6\";\n context.fillRect(0, 0, width, height);\n\n context.strokeStyle = \"#234\";\n context.fillStyle = \"#234\";\n\n context.beginPath();\n path(graticule);\n context.stroke();\n\n context.beginPath();\n path(land);\n context.fill();\n\n context.beginPath();\n path(sphere);\n context.stroke();\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"{\n const context = DOM.context2d(width, height);\n const proj = adamsSquareII().fitExtent([[50 , 50], [width - 50, height - 50]], sphere)\n const path = d3.geoPath(proj, context);\n\n context.fillStyle = \"#faf8f6\";\n context.fillRect(0, 0, width, height);\n\n context.strokeStyle = \"#234\";\n context.fillStyle = \"#234\";\n\n context.beginPath();\n path(graticule);\n context.stroke();\n\n context.beginPath();\n path(land);\n context.fill();\n\n context.beginPath();\n path(sphere);\n context.stroke();\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":158,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1013,"value":"md`---\n\nAdams, O. S. (1929), *Conformal Projection of the Sphere Within a Square*. Washington: U.S. Coast and Geodetic Survey Special Publication 153.\n\n\n\nBulirsch, R., (1965), *Numerical calculation of elliptic integrals and elliptic functions*. Numerische Mathematik 7, p. 78-90.\n\n\n\nSnyder, J. P. and Voxland, P. M. (1989). *An Album of Map Projections*. U.S. Geological Survey Professional Paper 1453. Washington, DC: United States Government Printing Office.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1008,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":538,"value":"adamsSquareI = {\n \n const { abs, max, min, sin, cos, asin, acos, tan } = Math;\n const pi = Math.PI, halfPi = pi / 2;\n\n const adamsSquareIRaw = function(lambda, phi) {\n let a, b, c, sm, sn;\n const sp = tan(0.5 * phi);\n c = cos(asin(sp)) * sin(0.5 * lambda);\n a = acos((c - sp) * Math.SQRT1_2);\n b = acos((c + sp) * Math.SQRT1_2);\n sm = lambda < 0;\n sn = phi < 0;\n return ellipticFactory(a, b, sm, sn);\n }\n \n adamsSquareIRaw.invert = function(x, y) {\n // Very rough initial guess (especially at high northings)\n // The magic values are got with:\n // forward(0, Math.PI / 2)\n // forward(Math.PI, 0)\n let phi = max(min(y / 1.8540746957596883, 1), -1) * halfPi;\n let lam = abs(phi) < pi ? max(min(x / 1.854074716833181, 1), -1) * pi : 0;\n return inverse(x, y, lam, phi, adamsSquareIRaw);\n }\n \n return () => d3.geoProjection(adamsSquareIRaw).scale(134.838125);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"adamsSquareII = {\n \n const { abs, max, min, sin, cos, asin, acos, tan } = Math;\n const pi = Math.PI, halfPi = pi / 2;\n \n const adamsSquareIIRaw = function(lambda, phi) {\n let a, b, sm, sn, xy;\n const sp = tan(0.5 * phi);\n a = cos(asin(sp)) * sin(0.5 * lambda);\n sm = (sp + a) < 0;\n sn = (sp - a) < 0;\n b = acos(sp);\n a = acos(a);\n\n xy = ellipticFactory(a, b, sm, sn);\n\n /* rotate by 45deg. */\n return [\n Math.SQRT1_2 * (xy[0] - xy[1]),\n Math.SQRT1_2 * (xy[0] + xy[1])\n ];\n }\n \n adamsSquareIIRaw.invert = function(x, y) {\n // Very rough initial guess (especially at high northings)\n // The magic values are got with:\n // forward(0, Math.PI / 2)\n // forward(Math.PI, 0)\n let phi = max(min(y / 2.62181347, 1), -1) * halfPi;\n let lam = abs(phi) < pi ? max(min(x / 2.62205760 / cos(phi), 1), -1) * pi : 0;\n return inverse(x, y, lam, phi, adamsSquareIIRaw);\n }\n \n return () => d3.geoProjection(adamsSquareIIRaw).scale(95.35384885);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":531,"value":"ellipticFactory = function(a, b, sm, sn) {\n let m = Math.asin(Math.sqrt(1 + Math.min(0, Math.cos(a + b))));\n if (sm) m = -m;\n\n let n = Math.asin(Math.sqrt(Math.abs(1 - Math.max(0, Math.cos(a - b)))));\n if (sn) n = -n;\n \n return [\n ellipticF(m, 0.5),\n ellipticF(n, 0.5)\n ];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":443,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":606,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":882,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"height = width","pinned":false,"hidden":false,"mode":"ojs"},{"id":38,"value":"sphere = ({ type: \"Sphere\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":30,"value":"graticule = d3.geoGraticule10()","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"land = topojson.feature(world, world.objects.land)","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"world = fetch(\"https://cdn.jsdelivr.net/npm/world-atlas@2/land-110m.json\").then(response => response.json())","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"d3 = require(\"d3-geo@1\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}