40:null 41:["$","title",null,{"children":"Quasi-Spilhaus square Projection | jie chen | Observable"}] 47:T74c,# Quasi-Spilhaus square Projection [This](https://observablehq.com/d/bb73e74c1685e498) is an implementation of a quasi-Spilhaus square map in Observable Notebook with [d3-geo](https://github.com/d3/d3-geo#projections), combined from work done by [Torben](https://observablehq.com/@toja/adams-world-in-a-square-i-ii) [Jansen](https://observablehq.com/@toja/spilhaus-world-ocean-map-in-a-square). However, this map is not the actual Spilhaus square projection but quasi (gridlines of longitude and latitude are 15'). For example, the upper right and lower left corners in the Spilhaus square projection must be at 30°N/115°E and 30°S/65°W, respectively. The reason for the slight error may be that d3 does not consider transforming the ellipsoidal reference system to the spherical reference system when performing map projections. Presently, only ArcGIS Pro 2.5+ supports the Spilhaus projection (see its [StoryMap](https://storymaps.arcgis.com/stories/756bcae18d304a1eac140f19f4d5cb3d), [template package](https://nation.maps.arcgis.com/home/item.html?id=9b2ce7c8179b4744af7bf3ddb86b7804), and [examples](https://www.esri.com/arcgis-blog/products/arcgis-pro/mapping/spilhaus-more-like-thrillhaus/) by an Esri cartographer [John Nelson](https://www.esri.com/arcgis-blog/author/j_nelson/)). See the implementation of the Spilhaus projection also in [PROJ](https://github.com/OSGeo/PROJ/blob/master/src/projections/adams.cpp), [Bertin](https://github.com/neocarto/bertin/blob/main/src/projections/spilhaus.js), and a personal website of [Kerkovits Krisztián](http://mercator.elte.hu/~kerkovits/projections/megnez.php?vetulet=adams&s1=115&s2=30&s3=-28.8). See comments and discussion on the Spilhaus projection in [GitHub issue](https://github.com/OSGeo/PROJ/issues/1851). This notebook relates to maps in [figshare](https://doi.org/10.6084/m9.figshare.21229757).48: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; } }49: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]; }4a:[] 45:{"id":"bb73e74c1685e498","slug":null,"title":"Quasi-Spilhaus square Projection","thumbnail":"b87da12b2c12b37f7fd3342eba62aa5a5402f8eeee34da871db83756a7c3c3ce","access_level":"public","create_time":"$D2022-12-12T22:53:29.617Z","update_time":"$D2023-06-12T15:12:32.285Z","view_time":null,"first_public_time":"$D2022-12-12T22:53:29.617Z","has_unpublished_edits":null,"fork":{"id":"f546aa0764b3720e","slug":"spilhaus-world-ocean-map-in-a-square","title":"Spilhaus World Ocean Map in a Square","owner":{"id":"03836a98946da550","login":"toja","name":"Torben Jansen"}},"stars":1,"starred":false,"owner":{"type":"team","id":"015bfd0008c9f4af","login":"chenjie0105mail","name":"jie chen","avatar_url":"https://avatars-next.observableusercontent.com/e2f3b72fffcedf34e8f64986ff1eb49eb32d769292546e4912bb962ee179b8de","tier":"starter_2024","bio":"","following":false},"version":535,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"fa07ff01a8d35588","login":"chenjie0105mail","name":"jie chen","avatar_url":"https://avatars-next.observableusercontent.com/e2f3b72fffcedf34e8f64986ff1eb49eb32d769292546e4912bb962ee179b8de","tier":null,"bio":"","following":false}],"members":[{"type":"individual","id":"fa07ff01a8d35588","login":"chenjie0105mail","name":"jie chen","avatar_url":"https://avatars-next.observableusercontent.com/e2f3b72fffcedf34e8f64986ff1eb49eb32d769292546e4912bb962ee179b8de","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":4,"collection":null,"type":"notebook","body":{"titling":"auto","title":"Quasi-Spilhaus square Projection","cells":[{"id":237,"value":"$47","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":56,"value":"{\n const context = DOM.context2d(width, height);\n const path = d3.geoPath(projection, context);\n\n context.fillStyle = \"#faf8f6\";\n context.fillRect(0, 0, width, height);\n\n context.strokeStyle = \"#234\";\n context.fillStyle = \"#234\";\n\n\n \n context.beginPath();\n path(land);\n context.lineWidth = 0.5;\n context.fillStyle = \"#999494\";\n context.fill(\"evenodd\");\n\n context.beginPath();\n path(sphere);\n context.stroke();\n\n context.beginPath();\n path(graticuleData);\n context.stroke();\n \n context.canvas.addEventListener('mousemove', pick);\n \n\n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":106,"value":"projection = spilhausSquare()\n .fitExtent([[1, 1], [width - 1, height - 1]], sphere);","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":198,"value":"// longitude and latitude of the mouse position\nlatlon ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":279,"value":"..........................................","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":204,"value":"spilhausSquare = {\n const center_lon = -66.94970198;\n const center_lat = 49.56371678;\n const azimuth = 40.17823482;\n \n const { abs, max, min, sin, cos, asin, acos, tan } = Math;\n const pi = Math.PI, halfPi = pi / 2;\n \n const spilhausSquareRaw = 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 return ellipticFactory(a, b, sm, sn);\n }\n \n spilhausSquareRaw.invert = function(x, y) {\n // Very rough initial guess (especially at high northings)\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, spilhausSquareRaw);\n }\n \n return () => d3.geoProjection(spilhausSquareRaw)\n .rotate([center_lon,center_lat, azimuth])\n .scale(1);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":215,"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","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":221,"value":"$48","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":217,"value":"$49","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":134,"value":"mutable latlon = [0, 0]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":139,"value":"function pick(event) {\n mutable latlon = projection.invert([event.layerX, event.layerY]);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":23,"value":"height = width","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":21,"value":"sphere = ({ type: \"Sphere\" })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":428,"value":"graticuleData = graticule();","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":18,"value":"graticule = d3.geoGraticule().step([15, 15])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"land = topojson.feature(world, world.objects.land)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":12,"value":"world = fetch(\"https://cdn.jsdelivr.net/npm/world-atlas@2/land-110m.json\").then(response => response.json())","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":2,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":7,"value":"d3 = require(\"d3-geo@1\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W4a","stdlib":"1"}} 46:"$45:body"