30:T9da,rhumb_map = { const mapwidth = (width/2 | 0) - 1, mapheight = mapwidth; const azimuthal = (() => { const forward_projection = d3.geoAzimuthalEquidistant().rotate([52]).precision(0.1) .fitExtent([[mapwidth*.02,mapheight*.02], [mapwidth*.98, mapheight*.98]], {type: "Sphere"}) const context = DOM.context2d(mapwidth, mapheight); const path = d3.geoPath(forward_projection, context); context.lineCap = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#eef7ff", context.fillRect(0, 0, mapwidth, mapheight); context.beginPath(), path(land50), context.fillStyle = "#9b8", context.fill(); context.beginPath(), path(graticule), context.strokeStyle = `rgba(70, 70, 70, 0.13)`, context.stroke(); context.restore(); context.save(); context.lineWidth = 3; context.beginPath(), path(great_circle_arc), context.strokeStyle = "#e44", context.stroke(); context.beginPath(), path(rhumb_line), context.strokeStyle = "#26b", context.stroke(); context.restore(); context.beginPath(), path({type: "Sphere"}), context.strokeStyle = "#000", context.stroke(); return context.canvas; })() const mercator = (() => { const forward_projection = d3.geoMercator().rotate([52, 1e-30]) .fitExtent([[mapwidth*.06,mapheight*.06], [mapwidth*.94, mapheight*.94]], {type: "Sphere"}) const context = DOM.context2d(mapwidth, mapheight); const path = d3.geoPath(forward_projection, context); context.lineCap = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#eef7ff", context.fillRect(0, 0, mapwidth, mapheight); context.beginPath(), path(land50), context.fillStyle = "#9b8", context.fill(); context.beginPath(), path(graticule), context.strokeStyle = `rgba(70, 70, 70, 0.13)`, context.stroke(); context.restore(); context.save(); context.lineWidth = 3; context.beginPath(), path(great_circle_arc), context.strokeStyle = "#e44", context.stroke(); context.beginPath(), path({type: "Sphere"}), context.clip(); context.beginPath(), path(rhumb_line), context.strokeStyle = "#26b", context.stroke(); context.restore(); context.beginPath(), path({type: "Sphere"}), context.strokeStyle = "#000", context.stroke(); return context.canvas; })(); return html` ${azimuthal}${mercator} ` }31:T52e,{ const mapwidth = Math.min(width, 680), mapheight = mapwidth; const forward_projection = d3.geoOrthographic().rotate([-150, -60, 0]).precision(.1) .fitExtent([[mapwidth*.02,mapheight*.02], [mapwidth*.98, mapheight*.98]], {type: "Sphere"}) const context = DOM.context2d(mapwidth, mapheight); const path = d3.geoPath(forward_projection, context); var render = (land) => { context.clearRect(0, 0, mapwidth, mapheight); context.lineCap = context.lineJoin = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#eef7ff", context.fillRect(0, 0, mapwidth, mapheight); context.beginPath(), path(land), context.fillStyle = "#bdb", context.fill(); context.lineWidth = 1.2; context.beginPath(), path(diagonals), context.strokeStyle = "#c44", context.stroke(); context.beginPath(), path(precise_graticule), context.strokeStyle = "#15a", context.stroke(); context.restore(); context.beginPath(), path({type: "Sphere"}), context.strokeStyle = "#000", context.stroke(); } d3.select(context.canvas) .call(zoom(forward_projection) .on("zoom.render", () => render(land110)) .on("end.render", () => render(land50))) .call(() => render(land50)) return context.canvas; }32:T5e6,{ const mapwidth = Math.min(width, 680), mapheight = mapwidth; const forward_projection = d3.geoOrthographic().rotate([-70, -35, 0]).precision(.1) .fitExtent([[mapwidth*.02,mapheight*.02], [mapwidth*.98, mapheight*.98]], {type: "Sphere"}) const context = DOM.context2d(mapwidth, mapheight); const path = d3.geoPath(forward_projection, context); var render = (countries) => { context.clearRect(0, 0, mapwidth, mapheight); context.lineCap = context.lineJoin = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#e0f0ff", context.fillRect(0, 0, mapwidth, mapheight); context.beginPath(), path(countries), context.fillStyle = "#f4fff7", context.fill(), context.strokeStyle = 'rgba(0,0,0,.3)', context.stroke(); context.beginPath(), path(graticule), context.strokeStyle = `rgba(30, 30, 30, 0.08)`, context.stroke(); context.fillStyle = 'rgba(255,40,60,.18)', context.strokeStyle = 'rgba(100,0,10,.3)', boxes.geometries.forEach(box => { context.beginPath(), path(box), context.fill(), context.stroke(); }); context.restore(); context.beginPath(), path({type: "Sphere"}), context.strokeStyle = "#000", context.stroke(); } d3.select(context.canvas) .call(zoom(forward_projection) .on("zoom.render", () => render(countries110)) .on("end.render", () => render(countries50))) .call(() => render(countries50)) return context.canvas; }33:Tb44,{ const mapwidth = (width/2 | 0) - 1, mapheight = mapwidth; const forward_projection = d3.geoOrthographic().rotate([45, 35.26, 0]).precision(.1) .fitExtent([[mapwidth*.02,mapheight*.02], [mapwidth*.98, mapheight*.98]], {type: "Sphere"}) let grid_line_width = .5; const sac = DOM.context2d(mapwidth, mapheight); var sacpath = d3.geoPath(forward_projection, sac); const render_sac = function(land) { const context = sac, path = sacpath; context.clearRect(0, 0, mapwidth, mapheight); context.lineCap = context.lineJoin = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#eef7ff", context.fillRect(0, 0, mapwidth, mapheight); context.beginPath(), path(land), context.fillStyle = "#bdb", context.fill(); context.lineWidth = grid_line_width; context.beginPath(), show_triangles ? path(lonlat_trigrid) : null, context.strokeStyle = "#c44", context.stroke(); context.beginPath(), show_hexagons ? path(lonlat_hexgrid) : null, context.strokeStyle = "#15a", context.stroke(); context.restore(); context.beginPath(), path({type: "Sphere"}), context.strokeStyle = "#000", context.stroke(); } const gnomonic = DOM.context2d(mapwidth, mapheight); const gnomonicpath = d3.geoPath(forward_projection, gnomonic); const render_gnomonic = function(land) { const context = gnomonic, path = gnomonicpath; context.clearRect(0, 0, mapwidth, mapheight); context.lineCap = context.lineJoin = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#eef7ff", context.fillRect(0, 0, mapwidth, mapheight); context.beginPath(), path(land), context.fillStyle = "#bdb", context.fill(); context.lineWidth = grid_line_width; context.beginPath(), show_triangles ? path(gnomonic_lonlat_trigrid) : null, context.strokeStyle = "#c44", context.stroke(); context.beginPath(), show_hexagons ? path(gnomonic_lonlat_hexgrid) : null, context.strokeStyle = "#15a", context.stroke(); context.restore(); context.beginPath(), path({type: "Sphere"}), context.strokeStyle = "#000", context.stroke(); } const render = function(land) { grid_line_width = .2 + (.4 + 3/grid_resolution) * forward_projection.scale() / 250; render_sac(land); render_gnomonic(land); } d3.select(sac.canvas) .call(zoom(forward_projection) .on("zoom.render", () => render(land110)) .on("end.render", () => render(land50))) .call(() => render(land50)) d3.select(gnomonic.canvas) .call(zoom(forward_projection) .on("zoom.render", () => render(land110)) .on("end.render", () => render(land50))) .call(() => render(land50)) return html` ${sac.canvas}${gnomonic.canvas} ` }34:T5d4,{ const mapwidth = 600, mapheight = 600; const forward_projection0 = d3.geoOrthographic().rotate([43.25, 35.55, 24.8]).precision(.1) .scale(2800).translate([mapwidth/2 + 200,mapheight/2 - 200]) const forward_projection = d3.geoOrthographic().rotate([0, 0, 0]).precision(.1) .scale(2800).translate([mapwidth/2 + 200,mapheight/2 - 200]) let grid_line_width = .5; const sac = DOM.context2d(mapwidth, mapheight); var sacpath = d3.geoPath(forward_projection, sac); var sacpath0 = d3.geoPath(forward_projection0, sac); const lonlat_hexgrid = inverseResampleJSON(sac_quincuncial_degrees, .01)({ type: 'MultiLineString', coordinates: hexgrid(16) }); const render_sac = function() { const context = sac, path = sacpath; context.clearRect(0, 0, mapwidth, mapheight); context.lineCap = context.lineJoin = "round"; context.save(); context.beginPath(), path({type: "Sphere"}), context.clip(), context.fillStyle = "#eef7ff", context.fillRect(0, 0, mapwidth, mapheight); context.lineWidth = grid_line_width; context.beginPath(), sacpath0(lonlat_hexgrid), context.strokeStyle = "#c44", context.stroke(); context.beginPath(), sacpath(lonlat_hexgrid), context.strokeStyle = "#15a", context.stroke(); context.restore(); } const render = function() { grid_line_width = .2 + (.4 + 3/grid_resolution) * forward_projection.scale() / 250; render_sac(); } render(); return sac.canvas; }35:T4eb,md` ## Introduction [D3-geo](https://github.com/d3/d3-geo/) approximates a projected great-circle arc (on a planar map, this is an arbitrary curve) by a “[polygonal chain](https://en.wikipedia.org/wiki/Polygonal_chain)” of straight line segments, using a technique it calls [“adaptive resampling”](https://bost.ocks.org/mike/example/): It decides how many segments to use by repeatedly bisecting the spherical arc, and testing how far the new projected midpoint is from the segment between previously projected endpoints. When this distance heuristic is satisfied, the curve is considered to be close enough to straight to visibly render as a line. (D3-geo also includes 2 additional secondary heuristic criteria to work around two classes of curves which foil this primary criterion.) If we instead want to approximate an arbitrary curve on the sphere (e.g. the inverse projection of a line segment on a map) by a chain of great-circle arcs, we can use the same technique with similar criteria. If we use a triple of Cartesian coordinates or a pair of stereographic coordinates to represent points on a sphere, we can avoid classical spherical trigonometry in implementing heuristic criteria similar to the ones D3-geo uses for planar resampling. `36:T86f,md`## Planar heuristic D3-geo’s adaptive resampling routine can be found in [\`projection/resample.js\`](https://github.com/d3/d3-geo/blob/v2.0.1/src/projection/resample.js#L27-L41). The heuristic for when to stop subdividing can be written as a compound of three predicate criteria ${$`\mathbf{a} \parallel \left( \mathbf{b} \mathrel{\&} \mathbf{c} \right)`} where ${$`\parallel`} means [logical *or*](https://en.wikipedia.org/wiki/Logical_disjunction) (inclusive disjunction) and ${$`\mathrel{\&}`} means [logical *and*](https://en.wikipedia.org/wiki/Logical_conjunction) (conjunction). If we call the projected start point ${$`p_0`}, the projected end point ${$`p_1`}, and the projected midpoint ${$`p_2`}, with displacement vectors ${$`v_1 = p_1 - p_0`} and ${$`v_2 = p_2 - p_0`}, then the predicates can be written: ${$$` \tfrac{1}{2}|v_1| < \delta \tag{\textbf{a}} `} ${$$` \frac{|v_2 \wedge v_1|}{|v_1|} < \delta \tag{\textbf{b}} `} ${$$` \tfrac{1}{5} < \frac{v_2 \cdot v_1}{|v_1|^2} < \tfrac{4}{5} \tag{\textbf{c}} `} Where ${$`|u \wedge v|`} is the area of a parallelogram formed by the two vectors ${$`u`} and ${$`v`}, ${$`|v|`} is the length of vector ${$`v`}, and ${$`u \cdot v`} is the [dot product](https://en.wikipedia.org/wiki/Dot_product) of vectors ${$`u`} and ${$`v`}. These conditions check whether ${$`(\mathbf{a})`} the length of the segment ${$`\overline{p_0p_1}`} is short enough that two segments of half its length would each be shorter than the visibility threshold ${$`\delta`}, ${$`(\mathbf{b})`} the perpendicular distance from segment ${$`\overline{p_0p_1}`} to the point ${$`p_2`} is sufficiently small that two segments would not be visually distinguishable from one, and ${$`(\mathbf{c})`} the projection of ${$`v_2`} onto the segment ${$`\overline{p_0p_1}`} is not too close to either of the endpoints. In combination, predicates ${$`(\mathbf{b} \mathrel{\&} \mathbf{c})`} specify that recursion-ending midpoints should fall into a segment-aligned rectangle. (D3-geo’s adaptive resampling also includes a fourth criterion that the spherical arc being projected must be smaller than 30° before subdivision stops.) `37:T623,md` ### Combining predicates in an ellipse Tests ${$`(\mathbf{b})`} and ${$`(\mathbf{c})`} can equivalently be expressed in a more symmetrical form. If we define ${$`v_3 = p_1 - p_2`}, we can write: ${$$` \frac{|v_2 \wedge v_3|}{|v_1|} < \delta \tag{\textbf{b}} `} ${$$` -\tfrac{3}{5} < \frac{ |v_3|^2 - |v_2|^2 }{|v_1|^2} < \tfrac{3}{5} \tag{\textbf{c}} `} The expression in the middle of the new test ${$`(\mathbf{c})`} above transforms us to a coordinate system where ${$`p_0`} is at ${$`-1`}, ${$`p_1`} is at ${$`+1`}, and we look at the projection of the point ${$`p_2`} onto that axis. This value is still linear in terms of the projected coordinate. Instead of absolute value we can use squaring for our tests (and indeed, D3-geo does, because finding lengths otherwise requires square roots): ${$$` \frac{|v_2 \wedge v_3|^2}{\delta^2 |v_1|^2} < 1 \tag{\textbf{b}} `} ${$$` \left(\tfrac{5}{3}\right)^{2} \left( \frac{ |v_3|^2 - |v_2|^2 }{|v_1|^2}\right)^{\!\! 2} < 1 \tag{\textbf{c}} `} Once we have these quadratic tests, we can combine them into a single new test ${$`(\mathbf{d})`}, turning a rectangle into an ellipse. I am replacing the width constant ${$`\bigl(\tfrac{5}{3}\bigr)^2`} with the slightly more lenient width constant ${$`2`} because we no longer need to worry about the corners of the rectangle: ${$$` \frac{|v_2 \wedge v_3|^2}{\delta^2 |v_1|^2} + 2 \left( \frac{ |v_3|^2 - |v_2|^2 }{|v_1|^2}\right)^{\!\! 2} < 1 \tag{\textbf{d}} `} To see this graphically, we can plot it using Desmos: https://www.desmos.com/calculator/vw04nfnidx `38:T960,Sphere = md` ## Adaptation for the sphere We can adapt this heuristic to work for approximating a curve on the sphere by “straight” great-circle arcs. ### Altitude predicate The primary predicate in the planar case is ${$`(\mathbf{b})`} which checks the altitude of the curve midpoint from the chord between start and end points. In the plane we calculated this by dividing the area of a parallelogram by the length of its side. On a sphere, if we call our points ${$`{P_0, P_1, P_2}`} we can calculate a good-enough approximation for the half-angle tangent of altitude ${$`h`} by dividing the half-angle tangent of spherical excess ${$`A`} of the spherical triangle by the half-angle tangent of great-circle arclength (${$`\ell`}): ${$$`\begin{aligned} \tan \tfrac{1}{2} \operatorname{Altitude}(P_2, \overline{P_0P_1}) &\approx \frac{\tan \tfrac{1}{2} \operatorname{Excess}(\triangle P_0P_1P_2)} {\tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0P_1})} \\[20mu] h &\approx \frac{A}{\ell} \end{aligned}`} This is not an exact formula unless triangle ${$`\triangle P_0P_1P_2`} has a right angle at ${$`P_0`} or ${$`P_1`}, and is a slight *underestimate* for the expected inputs with the projected midpoint roughly between the projected endpoints. If we want an approximation which is more accurate but still pretty computationally cheap, we can adjust our formula by multiplying it by a quadratic Taylor expansion of the error under the assumption that ${$`P_2`} is on the perpendicular bisector of segment ${$`\overline{P_0P_1}`}. Our spherical analog of planar criterion ${$`(\mathbf{b})`} becomes: ${$$` h \approx \frac{A\bigl(1 + \tfrac{1}{4}\ell^2 \bigr)}{\ell} \lt \tan \tfrac12 \delta \tag{$ \mathbf{b'} $} `} This is now an *overestimate*, which becomes very accurate when segment ${$`\overline{P_0P_1}`} is of moderate size and ${$`P_2`} is roughly halfway between them. ` // A useful identity from spherical trigonometry tells us that for a right-angled spherical triangle ${$`\triangle ABC`} with its right angle at point ${$`C`}, if we define ${$$`\begin{aligned} // a &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{BC}) \\[5mu] // b &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{AC}) // \end{aligned}`} // then we can compute the stereographic area as the product: // ${$$`ab = \tan \tfrac{1}{2} \operatorname{Excess}(\triangle ABC)`} 39:T1136,hidden_details(md` Let’s investigate more carefully how close the approximation is. First, we define a new operator for composing stereographic distances and areas: ${$$` x \oplus y = \frac {x + y} {1 - xy} `} If we call the foot of the altitude ${$`Q`} (somewhere along the great-circle between ${$`P_0`} and ${$`P_1`}), then we can break our spherical triangle ${$`\triangle P_0P_1P_2`} into two (right-triangle) parts ${$`\triangle P_0QP_2`} and ${$`\triangle QP_1P_2`}. If we call the stereographic distances: ${$$`\begin{aligned} h &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{QP_2}) \\[5mu] \ell_0 &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0Q}) \\[5mu] \ell_1 &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{QP_1}) \\[5mu] \ell = \ell_0 \oplus \ell_1 &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0P_1}) \end{aligned}`} Then the resulting stereographic areas are: ${$$`\begin{aligned} \ell_0h &= \tan \tfrac{1}{2} \operatorname{Excess}(\triangle P_0QP_2) \\[5mu] \ell_1h &= \tan \tfrac{1}{2} \operatorname{Excess}(\triangle QP_1P_2) \\[5mu] (\ell_0h) \oplus (\ell_1h) &= \tan \tfrac{1}{2} \operatorname{Excess}(\triangle P_0P_1P_2) \end{aligned}`} (Notice that the stereographic area of a right-angled spherical triangle is the product of the stereographic leg lengths. Also notice that both stereographic lengths and areas sum with the ${$`\oplus`} operation.) We want to find an approximation for ${$`h`}. If we start with stereographic coordinates for our points we can cheaply compute the area ${$`(\ell_0h) \oplus (\ell_1h)`} and the length ${$`\ell = \ell_0 \oplus \ell_1`}; we would like to avoid computing the precise position of ${$`Q`} or any other metrical quantities. The quotient of these is: ${$$`\begin{aligned} \tilde{h} = \frac{(\ell_0h) \oplus (\ell_1h)}{\ell_0 \oplus \ell_1} &= \frac{h(\ell_0 + \ell_1)}{1 - h^2\ell_0\ell_1} \Bigg/ \frac{\ell_0 + \ell_1}{1 - \ell_0\ell_1} \\[20mu] &= \frac{1 - \ell_0\ell_1}{1 - h^2\ell_0\ell_1} h \approx h \end{aligned}`} The approximation is exact whenever ${$`\ell_0 = 0`} or ${$`\ell_1 = 0`} (i.e. we started with a right-angled triangle), but is an underestimate if the projected midpoint is somewhere roughly between the endpoints. However, it gets better and better as the interval shrinks: spherical geometry is approximately Euclidean as we zoom closer. The worst case (greatest inaccuracy in this approximation) occurs where ${$`Q`} is on the perpendicular bisector of segment ${$`\overline{P_0P_1}`}. In that case, ${$`\ell_0 = \ell_1`}. We can assume that ${$`h \ll 1`} so the denominator ${$`1 - h^2\ell_0\ell_1 \approx 1`}. Then at worst our estimate ${$`\tilde{h}`} is an underestimate by approximately: ${$$` \frac{\tilde{h}}{h} \approx \frac 2 {1 + \sqrt{1 + \ell^2}} `} When ${$`\ell = 1`}, i.e. a quarter of the circle (90° angle, the distance from the pole to the equator), this worst-case ratio is ${$`\tilde{h} / h \approx 2 / (1 + \sqrt{2}) \approx 0.83.`} If we bisect once so that ${$`\ell = 1 / (1 + \sqrt{2})`} represents an eighth of a circle (45° angle), then ${$`\tilde{h} / h \approx 0.96.`} Bisect one more time (22.5° angle), and ${$`\tilde{h} / h \approx 0.99.`} The above is probably good enough, but if we want to guarantee that our true altitude is smaller than the specified tolerance, we can modify our approximation to always overestimate it instead. If we make a Taylor expansion of the reciprocal of our ${$`\ell_0 = \ell_1`} relative error estimate, we get: ${$$` \frac{h}{\tilde{h}} \approx \frac{1 + \sqrt{1 + \ell^2}}{2} = 1 + \frac{\ell^2}{4} - \frac{\ell^4}{16} + \mathcal{O}\bigl(\ell^6\bigr) `} We can multiply our approximation ${$`\tilde{h}`} by the first two terms of this series to refine it, turning it from an underestimate to an overestimate: ${$$`h \approx \frac{A\bigl(1 + \tfrac{1}{4}\ell^2 \bigr)}{\ell} `} At ${$`\ell = 1`} this formula overestimates ${$`h`} by less than 4% when ${$`P_2`} is on the perpendicular bisector of segment ${$`\overline{P_0P_1}`}, and the approximation improves very fast as ${$`\ell`} shrinks. `) // Relative to this composition operation, we can define a kind of halving operation: // ${$$`\begin{aligned} // x &= \sqrt[{\scriptsize \oplus}]{x} \oplus \!\sqrt[{\scriptsize \oplus}]{x} \\[5mu] // \sqrt[{\scriptsize \oplus}]{x} &= \frac x {1 + \sqrt{1 + x^2}} // \end{aligned}`} 3a:T57e,md` ### Width predicate For an analog of test ${$`(\mathbf{c})`}, we can also use stereographic lengths: ${$$`\begin{aligned} \ell_1 &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0P_1}) \\[5mu] \ell_2 &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0P_2}) \\[5mu] \ell_3 &= \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_2P_1}) \end{aligned}`} Then our spherical predicate can be (for some constant ${$`k_1`}): ${$$` \left|\ell_3^2 - \ell_2^2\right| < k_1\ell_1^2 \tag{$\mathbf{c'}$} `} ### Distance predicate One edge case that our predicates don’t properly handle is when the projected midpoint on the sphere is close to the midpoint of the *longer* great-circle arc between the projected endpoints; we want to guarantee it is close to the midpoint of the *shorter* great-circle arc instead. We also want to always bisect if the distance between points is large relative to the sphere. We can handle both of these by adding one more criterion ${$`\ell_3^2 + \ell_2^2 < k_2`} for some constant ${$`k_2`}. ### Combined predicate As in the planar case, we can combine the first two tests into one, defining an approximately elliptical shape (setting the two arbitrary width constant to ${$`2`}): ${$$` \frac{A^2\bigl(1 + \tfrac{1}{4}\ell_1^2 \bigr)^2}{\delta^2\ell_1^2} + 2 \biggl(\frac{\ell_3^2 - \ell_2^2}{\ell_1^2}\biggr)^{\!\!2} \lt 1 \tag{$ \mathbf{d'} $} `} `3b:T829,md` ## Computing stereographic length and area Given points on a unit-radius sphere represented by displacement vectors from the sphere center, e.g. the point ${$`P_0`} represented by vector ${$`u_0 = \langle x_0, y_0, z_0 \rangle`} with ${$`x_0^2 + y_0^2 + z_0^2 = 1`}, we can find the stereographic length between two points using the formula: ${$$`\begin{aligned} \tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0P_2}) &= \frac{|u_0 \wedge u_2|}{1 + u_0 \cdot u_2}\\[20mu] &= \frac{|(u_0 - u_2) \wedge u_2|}{(u_0 + u_2) \cdot u_2} \end{aligned}`} We can likewise find the stereographic area between three points using the formula: ${$$`\begin{aligned} \tan \tfrac{1}{2} \operatorname{Excess}(\triangle P_0P_1P_2) &= \frac{\left| u_0 \wedge u_1 \wedge u_2 \right|}{1 + u_0\cdot u_1 + u_1 \cdot u_2 + u_2 \cdot u_0}\\[20mu] &= \frac{\left| (u_0 - u_2) \wedge (u_1 - u_2) \wedge u_2 \right|}{(u_0 + u_2) \cdot (u_1 + u_2)} \end{aligned}`} Where ${$`\left| a \wedge b \wedge c \right|`} is the volume of the parallelepiped spanned by three vectors. If we are instead representing our points using [stereographic coordinates](/@jrus/planisphere), we need to define a new difference operator ${$`\boxminus`} for the position of a point after we have rotated the sphere to place another point at the north pole (where concatenation represents the [geometric product](https://en.wikipedia.org/wiki/Geometric_algebra)): ${$$`\begin{aligned} s_1 \boxminus s_0 &= \frac{s_1 - s_0}{1 + s_0 s_1} \\[20mu] &= \frac{(s_1 - s_0)(1 + s_0 \cdot s_1 - s_0 \wedge s_1)}{(1 + s_0 \cdot s_1)^2 + (s_0 \wedge s_1)^2} \end{aligned}`} Then we can easily find the stereographic length between two points: ${$$`\tan \tfrac{1}{2} \operatorname{Arclength}(\overline{P_0P_2}) = |s_0 \boxminus s_2|`} And we can find the stereographic area between three points using the formula: ${$$`\tan \tfrac{1}{2} \operatorname{Excess}(\triangle P_0P_1P_2) = \frac{\left| (s_0 \boxminus s_2) \wedge (s_1 \boxminus s_2) \right|}{1 + (s_0 \boxminus s_2) \cdot (s_1 \boxminus s_2)}`} `3c:Ta5d,inverseResampleJSON = function (projection, delta) { // delta = angle-measure tolerance, in degrees const maxDepth = 16, radians = Math.PI/180, dd = Math.tan(radians*delta/2) ** 2; var resampleLineTo = function( w0, u0, // w = input coordinates w1, u1, // u = cartesian coordinates on the sphere ll01, // ll = stereographic distance squared depth, array) { if (depth--) { var w2 = planar_midpoint(w0, w1), λφ2 = projection.invert(...w2), u2 = cartesian(λφ2), ll02 = stereo_length2(u2, u0), ll12 = stereo_length2(u2, u1), AA = stereo_area2(u2, u0, u1), hh = AA * (1 + 0.25*ll01)*(1 + 0.25*ll01) / (dd * ll01), // perpendicular projected distance ww = 2 * ((ll02 - ll12) / ll01) * ((ll02 - ll12) / ll01); // midpoint roughly in the middle // (ll02 + ll12 > 0.25) forces bisection for over-long segments, and handles ll01 == Infinity edge case // (ll02 + ll12 > dd)) stops bisection when the segment gets tiny if (((hh + ww > 1) & (ll02 + ll12 > dd)) | (ll02 + ll12 > 0.25)) { resampleLineTo(w0, u0, w2, u2, ll02, depth, array); array.push(λφ2); resampleLineTo(w2, u2, w1, u1, ll12, depth, array); } } } var resampleChain = function(pointarray) { let outarray = []; let w0 = pointarray[0], λφ0 = projection.invert(...w0), u0 = cartesian(λφ0); outarray.push(λφ0); for (var i = 1, n = pointarray.length; i < n; i++) { let w1 = pointarray[i], λφ1 = projection.invert(...w1), u1 = cartesian(λφ1); resampleLineTo(w0, u0, w1, u1, stereo_length2(u0, u1), maxDepth, outarray) outarray.push(λφ1); w0 = w1, u0 = u1; } return outarray; } var project = w => projection.invert(...w); var mapInPlace = fn => array => array.forEach((e, i) => array[i] = fn(e)); var convert, convertType = ({ Point: o => o.coordinates = project(o.coordinates), MultiPoint: o => mapInPlace(project)(o.coordinates), LineString: o => o.coordinates = resampleChain(o.coordinates), Polygon: o => mapInPlace(resampleChain)(o.coordinates), MultiLineString: o => mapInPlace(resampleChain)(o.coordinates), MultiPolygon: o => o.coordinates.forEach(mapInPlace(resampleChain)), Feature: o => convert(o.geometry), GeometryCollection: o => o.geometries.forEach(convert), FeatureCollection: o => o.features.forEach(convert), }); convert = o => (convertType?.[o?.type]?.(o), o); return function(json) { json = JSON.parse(JSON.stringify(json)); // make deep copy return convert(json); } }3d:T404,// apply the provided projection to every coordinate pair of a GeoJSON object projectJSON = function (projection) { var project = w => projection(...w); var projectChain = arr => arr.map(project); var mapInPlace = fn => array => array.forEach((e, i) => array[i] = fn(e)); var convert, convertType = ({ Point: o => o.coordinates = project(o.coordinates), MultiPoint: o => mapInPlace(project)(o.coordinates), LineString: o => o.coordinates = projectChain(o.coordinates), Polygon: o => mapInPlace(projectChain)(o.coordinates), MultiLineString: o => mapInPlace(projectChain)(o.coordinates), MultiPolygon: o => o.coordinates.forEach(mapInPlace(projectChain)), Feature: o => convert(o.geometry), GeometryCollection: o => o.geometries.forEach(convert), FeatureCollection: o => o.features.forEach(convert), }); convert = o => (convertType?.[o?.type]?.(o), o); return function(json) { json = JSON.parse(JSON.stringify(json)); // make deep copy return convert(json); } }3e:[] 2e:{"title":"Adaptive resampling on the sphere","cells":[{"id":0,"value":"md`# Adaptive resampling on the sphere`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1499,"value":"md`\nIf we want to plot a close approximation to a curve on the sphere corresponding to a straight line in some map projection, we can use this notebook’s \\`inverseResampleJSON\\` function to fit a chain of great-circle-arc segments approximating the inverse projection on the sphere of each original linear segment (in a [GeoJSON](https://geojson.org) object) to within a specified tolerance (here, expressed in degrees).\n\nFor example, it is straightforward to plot a [rhumb line](https://en.wikipedia.org/wiki/Rhumb_line) (line of constant compass bearing, in blue below) between [Santiago de Chile](https://en.wikipedia.org/wiki/Santiago) and [Reykjavík](https://en.wikipedia.org/wiki/Reykjav%C3%ADk) which circles the globe nearly 4 times westward along the way by inverse-projecting a line from the Mercator projection. Compare to the great-circle arc between the two cities (in red below).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":769,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":744,"value":"mercator_segment = projectJSON(mercator)({\n type: \"LineString\",\n coordinates: [\n [-70.67, -33.45], // Santiago\n [-21.88 - 360*4, 64.13]] }) // Reykjavík - 4 full turns","pinned":true,"hidden":false,"mode":"ojs"},{"id":1014,"value":"great_circle_arc = projectJSON(mercator.invert)(mercator_segment) // red line above","pinned":true,"hidden":false,"mode":"ojs"},{"id":886,"value":"rhumb_line = inverseResampleJSON(mercator, .05)(mercator_segment) // blue line above","pinned":true,"hidden":false,"mode":"ojs"},{"id":1670,"value":"md`\n* * *\nWe can easily produce a graticule matching a target tolerance or plot a diagonal grid by inverse-projecting straight segments drawn on an equirectangular projection:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1590,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1608,"value":"precise_graticule = inverseResampleJSON(d3.geoEquirectangularRaw, .01)({\n type: 'MultiLineString',\n coordinates: [].concat(\n d3.range(-180, 180, 10).map(λ => [[λ, -80], [λ, 80]]),\n d3.range(-80, 81, 10).map(φ => [[-180, φ], [180, φ]])) })","pinned":true,"hidden":false,"mode":"ojs"},{"id":1698,"value":"diagonals = inverseResampleJSON(d3.geoEquirectangularRaw, .01)({\n type: 'MultiLineString',\n coordinates: d3.range(-180, 180, 20).map(λ => [[λ, -90], [λ+180, 90], [λ+360, -90]]) })","pinned":true,"hidden":false,"mode":"ojs"},{"id":2716,"value":"md`\n* * *\n\nStill using the equirectangular projection, we can recreate something similar to Jason Davies’s [“Geographic Bounding Boxes”](https://www.jasondavies.com/maps/bounds/) plot:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2742,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":2718,"value":"bounds = countries110.features.map(feature => d3.geoBounds(feature))","pinned":true,"hidden":false,"mode":"ojs"},{"id":2720,"value":"boxes = inverseResampleJSON(d3.geoEquirectangularRaw, .02)({\n type: 'GeometryCollection',\n geometries: bounds.map(([[λ0, φ0], [λ1, φ1]]) => ({\n type: 'Polygon',\n coordinates: (φ0 === -90)\n ? [[[λ0, φ1], [λ1, φ1]]] // Antarctica\n : [[[λ0, φ0], [λ0, φ1], [λ1 += (λ1 < λ0) * 360, φ1], [λ1, φ0], [λ0, φ0]]] })) })","pinned":true,"hidden":false,"mode":"ojs"},{"id":1982,"value":"md`\n* * *\n\nResampling can be very useful for projecting a [global grid](https://en.wikipedia.org/wiki/Discrete_global_grid) back onto the globe. Compare spherical area coordinates ([Praun & Hoppe 2003](http://hhoppe.com/sphereparam.pdf); [Lei, Qi & Tian 2020](https://www.mdpi.com/2076-3417/10/2/655/htm)) for each octant on the left vs. a [gnomonic projection](https://en.wikipedia.org/wiki/Gnomonic_projection) for each octant on the right, below.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1919,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2032,"value":"viewof grid_resolution = Inputs.range([1, 32], {value: 16, step: 1, transform: Math.sqrt, label: \"Grid resolution\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2542,"value":"viewof show_hexagons = Inputs.toggle({value: true, label: \"Show hexagon grid\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2114,"value":"viewof show_triangles = Inputs.toggle({value: false, label: \"Show triangle grid\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3225,"value":"md`If we look at the spherical area coordinate grid near the center of an ocahedron face (red below) vs. near an octahedron vertex (blue below) we can see that the projection is not quite area-preserving, but does approximately preserve the length of the longest diagonal in each hexagon:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3125,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1840,"value":"hexgrid = {\n var map2 = (arr, cb) => arr.map(pts => pts.map(cb));\n // draw hexagon edges within one prototypical square:\n var square = [[[0, 0], [1/3, 1/3]], [[1, 0], [1/3, 1/3]], [[0, 1], [1/3, 1/3]]];\n return (n) => map2([].concat(\n ...d3.range(n).map(i => [].concat(\n ...d3.range(n).map(j => [].concat(\n map2(square, ([x, y]) => [i+x, j+y]), map2(square, ([x, y]) => [-i-x, -j-y]),\n map2(square, ([x, y]) => [i+1-x, -j-1+y]), map2(square, ([x, y]) => [-i-1+x, j+1-y]) )) ))\n ), ([x, y]) => [x/n, y/n]) }","pinned":true,"hidden":false,"mode":"ojs"},{"id":1922,"value":"lonlat_hexgrid = inverseResampleJSON(sac_quincuncial_degrees, .01)({\n type: 'MultiLineString',\n coordinates: hexgrid(grid_resolution) })","pinned":true,"hidden":false,"mode":"ojs"},{"id":8,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":23,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":447,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":144,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":328,"value":"md`\n#### Explanation \n\nJustifying the above formula gets a bit into the weeds. Click the button below to see a derivation:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":320,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":541,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":560,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":2677,"value":"spherical = {\n const degrees = 180/Math.PI;\n return ([x, y, z]) => [\n degrees*Math.atan2(y, x),\n degrees*Math.asin(z)]}","pinned":true,"hidden":false,"mode":"ojs"},{"id":646,"value":"cartesian = {\n const radians = Math.PI/180;\n return ([λ, φ]) => [\n Math.cos(radians*φ) * Math.cos(radians*λ),\n Math.cos(radians*φ) * Math.sin(radians*λ),\n Math.sin(radians*φ)]}","pinned":true,"hidden":false,"mode":"ojs"},{"id":719,"value":"md`\n${$$`\\begin{aligned}\n\\tan \\tfrac{1}{2} \\operatorname{Arclength}(\\overline{P_0P_1})\n&= \\frac{|u_0 \\wedge (u_1 - u_0)|}{u_0 \\cdot (u_1 + u_0)}\n\\end{aligned}`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":657,"value":"stereo_length2 = ([x0, y0, z0], [x1, y1, z1]) => {\n var \n pxy = x0 * (y1 - y0) - (x1 - x0) * y0,\n pyz = y0 * (z1 - z0) - (y1 - y0) * z0,\n pzx = z0 * (x1 - x0) - (z1 - z0) * x0,\n q = x0*(x1 + x0) + y0*(y1 + y0) + z0*(z1 + z0);\n return (pxy*pxy + pyz*pyz + pzx*pzx + !(q*q)) / (q*q); // adding !(q*q) means q==0 => return Infinity\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":689,"value":"md`\n${$$`\\begin{aligned}\n\\tan \\tfrac{1}{2} \\operatorname{Excess}(\\triangle P_0P_1P_2) \n&= \\frac{\\left| u_0 \\wedge (u_1 - u_0) \\wedge (u_2 - u_0) \\right|}{(u_0 + u_1) \\cdot (u_0 + u_2)}\n\\end{aligned}`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":682,"value":"stereo_area2 = ([x0, y0, z0], [x1, y1, z1], [x2, y2, z2]) => {\n var\n p = x0*((y1 - y0)*(z2 - z0) - (y2 - y0)*(z1 - z0)) +\n y0*((z1 - z0)*(x2 - x0) - (z2 - z0)*(x1 - x0)) +\n z0*((x1 - x0)*(y2 - y0) - (x2 - x0)*(y1 - y0)),\n q = (x0 + x2)*(x0 + x1) + (y0 + y2)*(y0 + y1) + (z0 + z2)*(z0 + z1);\n return (p*p + !(q*q)) / (q*q); // adding !(q*q) means q==0 => return Infinity\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":662,"value":"planar_midpoint = ([x0, y0], [x1, y1]) => [\n 0.5 * (x0 + x1),\n 0.5 * (y0 + y1)]","pinned":true,"hidden":false,"mode":"ojs"},{"id":697,"value":"md`\n${$$`\n\\frac{A^2\\bigl(1 + \\tfrac{1}{4}\\ell_1^2 \\bigr)^2}{\\delta^2\\ell_1^2}\n+ 2 \\biggl(\\frac{\\ell_3^2 - \\ell_2^2}{\\ell_1^2}\\biggr)^{\\!\\!2}\n\\lt 1\n`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1016,"value":"$3c","pinned":true,"hidden":false,"mode":"ojs"},{"id":875,"value":"$3d","pinned":true,"hidden":false,"mode":"ojs"},{"id":726,"value":"md`\n* * *\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1051,"value":"// mercator that takes lon/lat input in degrees\nmercator = {\n var radians = Math.PI/180;\n var degrees = 180/Math.PI;\n var proj = (λ, φ) => d3.geoMercatorRaw(λ*radians, φ*radians);\n proj.invert = (x, y) => {x = d3.geoMercatorRaw.invert(x, y); return [x[0]*degrees, x[1]*degrees]};\n return proj;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2089,"value":"import {zoom} from '@d3/versor-zooming'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1940,"value":"trigrid = {\n var map2 = (arr, cb) => arr.map(pts => pts.map(cb));\n var square = [[[1/3, 0], [0, 1/3]], [[1/3, 1], [1, 1/3]], [[2/3, 0], [2/3, 1]], [[0, 2/3], [1, 2/3]]];\n return (n) => map2([].concat(\n ...d3.range(n).map(i => [].concat(\n ...d3.range(n).map(j => [].concat(\n map2(square, ([x, y]) => [+i+x, +j+y]),\n map2(square, ([x, y]) => [-i-x, -j-y]),\n map2(square, ([x, y]) => [+i+1-x, -j-1+y]),\n map2(square, ([x, y]) => [-i-1+x, +j+1-y])\n ))\n ))\n ), ([x, y]) => [x/n, y/n])\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1907,"value":"sac_quincuncial_degrees = {\n var radians = Math.PI/180;\n var degrees = 180/Math.PI;\n var projection = (λ, φ) => sac_quincuncial(λ*radians, φ*radians);\n projection.invert = (x, y) => {x = sac_quincuncial.invert(x, y); return [x[0]*degrees, x[1]*degrees]}\n return projection;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2012,"value":"gnomonic_quincuncial_degrees = {\n var radians = Math.PI/180;\n var degrees = 180/Math.PI;\n var projection = (λ, φ) => gnomonic_quincuncial(λ*radians, φ*radians);\n projection.invert = (x, y) => {x = gnomonic_quincuncial.invert(x, y); return [x[0]*degrees, x[1]*degrees]}\n return projection;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2010,"value":"gnomonic_lonlat_trigrid = inverseResampleJSON(gnomonic_quincuncial_degrees, .01)({\n type: 'MultiLineString',\n coordinates: trigrid(grid_resolution)\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2015,"value":"gnomonic_lonlat_hexgrid = inverseResampleJSON(gnomonic_quincuncial_degrees, .01)({\n type: 'MultiLineString',\n coordinates: hexgrid(grid_resolution)\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1953,"value":"lonlat_trigrid = inverseResampleJSON(sac_quincuncial_degrees, .01)({\n type: 'MultiLineString',\n coordinates: trigrid(grid_resolution)\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1899,"value":"import {quincuncial as sac_quincuncial} from '@jrus/sac-quincuncial'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2006,"value":"import {quincuncial as gnomonic_quincuncial} from '@jrus/gnomonic-quincuncial'","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"import {math_css, $, $$} with {$css as $css} from \"@jrus/misc\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"$$css = html`${math_css}`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":316,"value":"hidden_details = (content) => {\n var button = Object.assign(\n html`