2b:null 2c:["$","title",null,{"children":"Robust Raster Reprojection | Jeshurun Hembd | Observable"}] 31:T537,viewof mapContext = { let div = html`

`; // The map canvas will have its drawingbuffer resized by tile-frame. // For display here, we just want to make sure the CSS size fills the viewport const displaySize = Math.max(width, globeDisplayHeight); const context = DOM.context2d(displaySize, displaySize); // Make sure the display is roughly centered within the viewport // Any overflow will be clipped since the parent div has overflow: hidden context.canvas.style.position = "relative"; let yShift = -0.35 * (displaySize - globeDisplayHeight); let xShift = -0.5 * (displaySize - width); context.canvas.style.top = yShift + "px"; context.canvas.style.left = xShift + "px"; // Add the canvas to the display div, and make sure it is underneath the attribution context.canvas.style.zIndex = 1; div.appendChild(context.canvas); let attribution = html`

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL

`; div.appendChild(attribution); div.value = context; return div; }32:T478,viewof camPosition = { const form = html``; form.value = [0, 0, 0]; form.oninput = () => { form.value[0] = form.a.valueAsNumber / degrees; //form.oa.value = `${form.a.valueAsNumber.toFixed(2)}°`; form.value[1] = form.b.valueAsNumber / degrees; //form.ob.value = `${form.b.valueAsNumber.toFixed(2)}°`; form.value[2] = Math.pow(2, form.c.valueAsNumber); form.oc.value = `${form.value[2].toString().substring(0, 5)} km`; let hNorm = form.value[2] / earthRadius; let latScale = 3 / (Math.cos(form.value[1]) + 0.1); form.a.step = (latScale * hNorm).toFixed(5); form.b.step = (2 * hNorm).toFixed(5); }; form.oninput(); return form; }33:T572,viewof glslFunctions = { const form = html``; form.value = {}; form.oninput = () => { form.value.diffSqrt = form.a.checked ? diffSqrt.optimal.src : diffSqrt.direct.src; form.oa.innerHTML = ""; form.oa.append(form.a.checked ? diffSqrt.optimal.tex : diffSqrt.direct.tex); form.value.latChange = form.b.checked ? latChange.optimal.src : latChange.direct.src; form.ob.innerHTML = ""; form.ob.append( form.b.checked ? latChange.optimal.tex : latChange.direct.tex ); form.value.log1plusX = form.c.checked ? log1plusX.optimal.src : log1plusX.direct.src; form.oc.innerHTML = ""; form.oc.append( form.c.checked ? log1plusX.optimal.tex : log1plusX.direct.tex ); form.value.smallTan = form.d.checked ? smallTan.optimal.src : smallTan.direct.src; form.od.innerHTML = ""; form.od.append(form.d.checked ? smallTan.optimal.tex : smallTan.direct.tex); }; form.oninput(); return form; }34:T445,md`### How does it work? The basic workflow is fairly simple. As [Jake Low](https://observablehq.com/@jake-low) explained in [Raster reprojection](https://observablehq.com/@jake-low/raster-reprojection), we loop over the pixels in the output image, and do the following for each pixel: 1. Inverse project the screen x/y coordinates to longitude and latitude 2. Project longitude/latitude to x/y in the input map coordinates (Web Mercator) 3. Get the color of the input map at that location, and copy to the output pixel In this notebook, step 1 uses the [General Perspective Projection](https://en.wikipedia.org/wiki/General_Perspective_projection), since we want a realistic view of the Earth [as it appears from space](https://observablehq.com/@jjhembd/the-view-from-space). The flow is also rewritten in WebGL, for faster performance. A direct implementation of all 3 steps works well, *as long as we don't zoom in too far*. And WebGL makes it very fast. Even though it is looping over thousands of pixels, it uses your device's GPU to process hundreds of pixels at a time, in parallel.`35:T5ab,md`### Why is this a problem? In WebGL, numbers are stored in 32-bit blocks of memory, in [single-precision floating point format](https://en.wikipedia.org/wiki/Single-precision_floating-point_format). This means each number has 7-8 [significant decimal digits](https://en.wikipedia.org/wiki/Significant_figures), plus an exponent. For example, in the first formula above: ${tex`1 - \sqrt{1 - \epsilon}`}, the values ${tex`1`} and ${tex`\epsilon`} might be stored something like this: ${tex.block`1 \mapsto +1.0000000 \times 10^{+000}`} ${tex.block`\epsilon \mapsto +2.2537194 \times 10^{-005}`} Note how our example ${tex`\epsilon`} has 8 digits of information. To subtract these two numbers, we have to start by making the exponents the same, so the decimal points will line up: ${tex.block`\epsilon_{equalized} \mapsto 0.0000225 \times 10^{+000}`} ${tex.block`1 - \epsilon_{equalized} \mapsto \quad \begin{aligned} 1.0000000 \times 10^{+000} \\ \underline{- 0.0000225 \times 10^{+000}}\\ = 0.9999774 \times 10^{+000} \end{aligned}`} Equalizing the exponents already truncates most of the digits of ${tex`\epsilon`}. The result of the subtraction is written to a new 32-bit block of memory before computing the square root and subtracting from 1. The information in the truncated digits is therefore permanently lost. This is a classic case of [subtractive cancellation](https://en.wikipedia.org/wiki/Round-off_error#Subtractive_cancellation). `36:T4a3,md`### What can we do about it? We could do the calculations at higher precision. Since WebGL doesn't support double precision math, we would need to do some complicated tricks, like using two floating point numbers for each variable. Henry Thasler has a [nice demo of this approach](https://www.thasler.com/blog/blog/glsl-part2-emu). But besides making the code more complicated, this "emulated double precision" approach would also be [significantly slower](http://blog.hvidtfeldts.net/index.php/2012/07/double-precision-in-opengl-and-webgl/). Let's take a step back. The problems arise when we have two very differently-sized numbers in the same formula. What if the formula could be re-written using only small terms? This is where a [Taylor series](https://en.wikipedia.org/wiki/Taylor_series) could help. Consider this expansion around ${tex`\epsilon = 0`}: ${tex.block`1 - \sqrt{1-\epsilon} \space \approx \space \frac{\epsilon}{2} + \frac{\epsilon^2}{8} + \frac{\epsilon^3}{16} + \mathcal{O}(\epsilon^4).`} The large constants (the 1's) disappeared, leaving only a sum of small terms. Even in single-precision math, all the information in ${tex`\epsilon`} will be preserved.`37:T49e,// Default math font is too big html``38:T612,drawGlobe = { const progInfo = yawgl.initShaderProgram( gl, vertexShaderSrc, fragmentShaderSrc ); // Load data into GPU for shaders: attribute buffers, indices, textures const buffers = yawgl.initQuadBuffers(gl); const texture = yawgl.initTexture(gl, mapSize, mapSize); // Store links to uniform arrays const uniforms = { uMaxRay: new Float64Array(2), uTextureSampler: texture.sampler, uCamMapPos: new Float64Array(2), uMapScale: new Float64Array(2) }; return function (camPos) { // Update uniforms related to camera position uniforms.uHnorm = camPos[2] / earthRadius; uniforms.uLat0 = camPos[1]; uniforms.uCosLat0 = Math.cos(camPos[1]); uniforms.uSinLat0 = Math.sin(camPos[1]); uniforms.uTanLat0 = Math.tan(camPos[1]); [uniforms.uExpY0, uniforms.uLatErr] = getProjFactors(camPos[1]); let aspect = width / globeDisplayHeight; uniforms.uMaxRay.set([tanHalfFOV * aspect, tanHalfFOV]); // Set uniforms and update texture for the map uniforms.uCamMapPos.set(camMapPosition); uniforms.uMapScale.set([map.getScale(0), map.getScale(1)]); texture.update(map.canvas); // Draw the globe let pixRatio = devicePixelRatio > 1 ? 2 : 1; yawgl.resizeCanvasToDisplaySize(gl.canvas, pixRatio); yawgl.drawScene(gl, progInfo, buffers, uniforms); return; }; function getProjFactors(camLatitude) { var clipLat = clipLatitude(camLatitude); var expY = Math.tan(0.25 * Math.PI + 0.5 * clipLat); return [expY, camLatitude - clipLat]; } }39:[] 2f:{"id":"623ba5541ee7e1b3","slug":"robust-raster-reprojection","title":"Robust Raster Reprojection","thumbnail":"0f0952da45dfb2233ad196afbbc5e806b93cbfe43c3a3a17a69104082eef6663","access_level":"public","create_time":"$D2019-10-10T23:30:17.175Z","update_time":"$D2022-07-26T13:11:48.602Z","view_time":null,"first_public_time":"$D2019-10-18T19:48:34.212Z","fork":false,"stars":43,"starred":false,"owner":{"type":"team","id":"54b0aa2cda468b5a","login":"jjhembd","name":"Jeshurun Hembd","avatar_url":"https://avatars-next.observableusercontent.com/d42857c40a13c8179e0491536477452e54d4f6808ec32ed4d8fbdad81c03758d","tier":"starter_2024"},"version":2284,"publish_version":2284,"role":null,"authored":false,"authors":[{"type":"individual","id":"4e65ca212a169cf2","login":"jjhembd","name":"Jeshurun Hembd","avatar_url":"https://avatars-next.observableusercontent.com/d42857c40a13c8179e0491536477452e54d4f6808ec32ed4d8fbdad81c03758d","tier":null}],"members":[{"type":"individual","id":"4e65ca212a169cf2","login":"jjhembd","name":"Jeshurun Hembd","avatar_url":"https://avatars-next.observableusercontent.com/d42857c40a13c8179e0491536477452e54d4f6808ec32ed4d8fbdad81c03758d","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":18,"collection":null,"type":"notebook","body":{"title":"Robust Raster Reprojection","cells":[{"id":0,"value":"md`# Robust Raster Reprojection`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":695,"value":"md`We want to take a tiled web map, like this:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":807,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":698,"value":"md`and wrap it around a globe, like this:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":818,"value":"viewof gl = {\n const canvas = html`

`;\n const ctx = canvas.getContext(\"webgl\");\n ctx.getExtension('OES_standard_derivatives');\n canvas.value = ctx;\n return canvas;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":944,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1340,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":701,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":708,"value":"md`### Zooming into trouble\nBy the time the camera gets close enough to see streets and buildings, the image starts to fall apart. Try zooming in to your neighborhood on the above globe, with all the \"optimizations\" turned off. Here is an example I found, zooming in to the coordinates of the [Tower of London](https://en.wikipedia.org/wiki/Tower_of_London):`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1616,"value":"{\n let img1 = html`

Direct implementation
(click to toggle)

`;\n let img2 = html`

With optimizations
(click to toggle)

`;\n let toggleBox = html`

`;\n toggleBox.appendChild(img1);\n toggleBox.onclick = () => {\n let oldImg = toggleBox.removeChild(toggleBox.firstChild);\n let newImg = oldImg == img2 ? img1 : img2;\n toggleBox.appendChild(newImg);\n };\n return toggleBox;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":714,"value":"md`The direct implementation gives an image in the wrong place—somewhere in the River Thames. And the image itself is strangely aliased.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1811,"value":"md`### What is going on?\n\nThe math for our \"simple\" workflow requires us to calculate quantities like these:\n* ${tex`1 - \\sqrt{1 - \\epsilon}`}\n* ${tex`\\sin^{-1}\\left(\\phi_0 + \\epsilon\\right) - \\phi_0`}\n* ${tex`\\ln\\left(1 + \\epsilon/e^{Y_0}\\right) - \\ln\\left(1 - \\epsilon \\space e^{Y_0}\\right)`}\n\nIn each case, ${tex`\\epsilon`} can become very small as the camera approaches the Earth surface. Since the ${tex`\\epsilon`} factors are inside another expression (like ${tex`\\sqrt{1-\\epsilon}`}), they will be first combined with a large number. Then, the last step of the calculation is a subtraction of two similarly-sized numbers.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1755,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":754,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2015,"value":"md`### Approximations and accuracy\nThis truncated Taylor series is an approximation. Will it be accurate enough? \n\nFor full single-precision accuracy, the first neglected term needs to be \"small enough\" that including it would not affect the result. The single-precision [significand](https://en.wikipedia.org/wiki/Significand) has 24 bits, so \"small enough\" is defined as follows:\n${tex.block`\\frac{5\\epsilon^4}{128} < \n2^{-24}\\left(\\frac{\\epsilon}{2} + \\frac{\\epsilon^2}{8} + \\frac{\\epsilon^3}{16}\\right)\n\\implies \\epsilon < 0.01`}\n\nBut how accurate does it really need to be? It just needs to find the correct pixel from the input map. The input map is less than 10,000 pixels across. Three terms of the Taylor series will give 1/10,000 accuracy, even for inputs as large as ${tex`\\epsilon < 0.1`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2035,"value":"md`### Putting it all together\nThe above example outlines a strategy for robust reprojection in WebGL. The rest of this notebook goes through all the equations in the flow and translates them into [GLSL](https://webglfundamentals.org/webgl/lessons/webgl-shaders-and-glsl.html) code snippets. \n\nWhere the equations are numerically unstable, we rework them with Taylor expansions to isolate the small terms. We then create \"optimal\" versions of the GLSL snippets, that switch to the Taylor expansion when the small term is below a threshold.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1188,"value":"md`---\n## Math and Code\nHere is the basic workflow in [GLSL](https://webglfundamentals.org/webgl/lessons/webgl-shaders-and-glsl.html):`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":861,"value":"fragmentMainSrc = `\nvoid main(void) {\n // 0. Pre-compute some values\n float sinC = calcSinC(length(vRayParm));\n float cosC = sqrt(1.0 - sinC * sinC);\n\n // 1. Invert for longitude and latitude perturbations relative to camera\n vec2 dLonLat = xyToLonLat(vRayParm, sinC, cosC);\n\n // 2. Project dLon, dLat to perturbations in world Mercator coordinates\n vec2 dMercator = projMercator(dLonLat);\n // Transform to current map zoom & pan\n vec2 texCoord = uCamMapPos + uMapScale * dMercator;\n \n // 3. Lookup color from the texture, add cosine shading and horizon tapering\n vec4 texelColor = texture2D(uTextureSampler, texCoord);\n float taper = horizonTaper(length(vRayParm));\n gl_FragColor = vec4(texelColor.rgb * cosC * taper, texelColor.a);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1306,"value":"md`Step 1 uses the [General Perspective Projection](https://en.wikipedia.org/wiki/General_Perspective_projection). The inverse equations for this projection come from John Snyder's [Map Projections, A Working Manual](https://pubs.usgs.gov/pp/1395/report.pdf), page 175. Following Snyder, we begin by pre-computing a ${tex`\\sin c`} factor.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1311,"value":"md`### 0. Pre-computing sinC\nSnyder defines the angular distance ${tex`c`} along the sphere from the map center with this formula:\n${tex.block`c = \\sin^{-1}\n \\left[\n \\frac{ P - \\sqrt{1 - \\frac{\\rho^2(P+1)}{R^2(P-1)}} }\n { \\frac{R(P-1)}{\\rho} + \\frac{\\rho}{R(P-1)} }\n \\right],`}\nwhere ${tex`\\rho = \\sqrt{x^2 + y^2}`} is the radial distance in map coordinates from the map center, and ${tex`P = 1 + h/R`} is the camera's distance from the Earth center, normalized by the Earth radius ${tex`R`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":145,"value":"md`The formula for ${tex`c`} is complicated, and undefined at ${tex`\\rho = 0`}. It can be simplified by defining the tangents of the ray angle relative to the vertical,\n${tex.block`\\vec{p} = \\frac{1}{h}\\left(x,y\\right), \\quad p = |\\vec{p}| = \\sqrt{p_x^2 + p_y^2} = \\frac{\\rho}{h},`} \nand the normalized altitude ${tex`h_n = h / R`}. A little algebra now yields`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":166,"value":"tex.block`\\sin c = \n \\frac{p}{1 + p^2} \n \\left(h_n + 1 - \\sqrt{1 - p^2 \\thinspace h_n \\left(2 + h_n\\right)}\\right).`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1331,"value":"md`Translating into GLSL:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1094,"value":"calcSinCSrc = `\nuniform float uHnorm;\nfloat calcSinC(float p) {\n float p2 = p * p;\n float gamma = p2 * uHnorm * (2.0 + uHnorm);\n return (uHnorm + diffSqrt(gamma)) * p / (1.0 + p2);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1115,"value":"md`where diffSqrt(x) is a custom function that computes ${tex`1 - \\sqrt{1 - x}`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1335,"value":"diffSqrt_direct = `\nfloat diffSqrt(float x) {\n return 1.0 - sqrt(1.0 - x);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1375,"value":"md`The argument of diffSqrt(x) is the term ${tex`p^2 h_n (2 + h_n)`}. The square root of this term is the ray tangent of the current pixel divided by the tangent of the angle to the horizon. When the camera is close to the surface (and pointing straight down), this term can be very small. In this case, it is more robust to use a few terms of the Taylor series:\n${tex.block`1 - \\sqrt{1-x} \\space \\approx \\space \\frac{x}{2} + \\frac{x^2}{8} + \\frac{x^3}{16} + \\mathcal{O}(x^4).`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":870,"value":"diffSqrt_optimal = `\nfloat diffSqrt(float x) {\n float halfx = 0.5 * x;\n return (x < 0.1)\n ? halfx * (1.0 + 0.5 * halfx * (1.0 + halfx))\n : 1.0 - sqrt(1.0 - x);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1388,"value":"diffSqrt = ({\n direct: {\n src: diffSqrt_direct,\n tex: tex`1-\\sqrt{1-x}`\n },\n optimal: {\n src: diffSqrt_optimal,\n tex: tex`\\frac{x}{2} + \\frac{x^2}{8} + \\frac{x^3}{16}`\n }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":18,"value":"md`### 1. Inverse projection to longitude and latitude\nWe assume that the globe view is being constructed for a \"point of perspective\" (camera position) at longitude, latitude, and altitude ${tex`(\\lambda_0, \\phi_0, h)`}. Now, for a point on the map with coordinates ${tex`(x, y)`}, we want to *invert* the projection to find the corresponding longitude ${tex`\\lambda`} and latitude ${tex`\\phi`} on the globe.\n\nHere are Snyder's inverse equations, restated as *perturbations* relative to the camera coordinates:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":288,"value":"tex.block`\n\\Delta \\lambda = \\lambda - \\lambda_0 = \n \\tan^{-1}\\left(\\frac{\\hat{p}_x \\sin c}{\\quad \\cos \\phi_0 \\cos c - \\hat{p}_y \\sin \\phi_0 \\sin c}\\right)`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":302,"value":"tex.block`\\Delta \\phi = \\phi - \\phi_0 = \n \\sin^{-1}\\left(\\sin \\phi_0 \\cos c + \\hat{p}_y \\sin c \\cos \\phi_0\\right) - \\phi_0`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":276,"value":"md`where ${tex`\\hat{p}`} is the normalized map coordinate\n${tex.block`\\hat{p} \\stackrel{\\tiny{\\text{def}}}{=} \\frac{\\vec{p}}{|p|} = \\frac{1}{\\rho} \\begin{pmatrix}x \\\\ y\\end{pmatrix}.`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":209,"value":"md`Translating to GLSL, note that ${tex`\\phi_0`}, ${tex`\\cos \\phi_0`}, and ${tex`\\sin \\phi_0`} are set as [uniforms](https://webglfundamentals.org/webgl/lessons/webgl-shaders-and-glsl.html#uniforms), which are passed along to every fragment as constant values.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":873,"value":"xyToLonLatSrc = `\nvec2 xyToLonLat(vec2 xy, float sinC, float cosC) {\n vec2 pHat = normalize(xy);\n float dLon = atan(pHat.x * sinC, uCosLat0 * cosC - pHat.y * uSinLat0 * sinC);\n float dLat = latChange(pHat.x, pHat.y, sinC, cosC);\n return vec2(dLon, dLat);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1402,"value":"latChange_direct = `\nfloat latChange(float x, float y, float sinC, float cosC) {\n return asin(uSinLat0 * cosC + y * uCosLat0 * sinC) - uLat0;\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1409,"value":"md`The latitude formula has been delegated to a separate function, since it has potential for numerical issues. When ${tex`\\sin c`} is small, the formula is essentially ${tex`\\sin^{-1} \\left( \\sin \\phi_0 + \\epsilon \\right) - \\phi_0`}. In this case, a Taylor expansion better preserves the precision:\n${tex.block`\\Delta \\phi \\approx \\hat{p}_y\\sin c - \\frac{1}{2}\\hat{p}_x^2\\tan\\phi_0\\sin^2 c \n- \\left( \\frac{1}{2}\\tan^2\\phi_0\\hat{p}_y\\hat{p}_x^2 - \\frac{1}{6}\\hat{p}_y^3 \\right) \\sin^3 c\n+ \\mathcal{O}(\\sin^4 c)`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1400,"value":"latChange_optimal = `\nfloat latChange(float x, float y, float sinC, float cosC) {\n float xtan = x * uTanLat0;\n float curveTerm = 0.5 * y * (xtan * xtan - y * y / 3.0);\n return (max( sinC, abs(sinC * uTanLat0) ) < 0.1)\n ? sinC * (y - sinC * (0.5 * xtan * x + curveTerm * sinC))\n : asin(uSinLat0 * cosC + y * uCosLat0 * sinC) - uLat0;\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1422,"value":"latChange = ({\n direct: {\n src: latChange_direct,\n tex: tex`\\sin^{-1}\\left(\\sin \\phi_0 \\cos c + \\hat{p}_y \\sin c \\cos \\phi_0\\right) - \\phi_0`\n },\n optimal: {\n src: latChange_optimal,\n tex: tex`\\hat{p}_y\\sin c - \\frac{1}{2}\\hat{p}_x^2\\tan\\phi_0\\sin^2 c - \\left( \\frac{1}{2}\\tan^2\\phi_0\\hat{p}_y\\hat{p}_x^2 - \\frac{1}{6}\\hat{p}_y^3 \\right) \\sin^3 c`\n }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":181,"value":"md`### 2. Projection to map x, y\nWe define our global X-coordinate as ranging from ${tex`X_M\\nobreak=\\nobreak0`} at \n${tex`\\lambda\\nobreak=\\nobreak-180\\degree`} to ${tex`X_M\\nobreak=\\nobreak1`} \nat ${tex`\\lambda\\nobreak=\\nobreak180\\degree`}. The perturbation in X is therefore\n${tex.block`\\Delta X_M = \\frac{\\Delta \\lambda}{2 \\pi}.`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":342,"value":"md`The Mercator Y coordinate is given by\n${tex`Y_M = \\frac{1}{2\\pi} \\ln \\left[ \\tan\\left(\\frac{\\pi}{4} + \\frac{\\phi}{2} \\right) \\right],`}\nso a perturbation relative to the camera position ${tex`\\phi_0`} will be\n${tex.block`2 \\pi \\space \\Delta Y_M = \n \\ln \\left[ \\tan\\left(\\frac{\\pi}{4} + \\frac{\\phi_0 + \\Delta\\phi}{2}\\right) \\right] \n - \\ln \\left[ \\tan\\left(\\frac{\\pi}{4}+\\frac{\\phi_0}{2}\\right) \\right].`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":355,"value":"md`After defining \n${tex`e^{Y_0} \\stackrel{\\tiny{\\text{def}}}{=} \\tan\\left(\\frac{\\pi}{4} + \\frac{\\phi_0}{2}\\right)`}\nand \n${tex`\\delta \\stackrel{\\tiny{\\text{def}}}{=} \\tan\\left(\\frac{\\Delta\\phi}{2}\\right)`} and using the identities\n${tex`\\tan\\left(\\alpha + \\beta\\right) = \\frac{\\tan\\alpha + \\tan\\beta}{1 - \\tan\\alpha \\tan\\beta}`} \nand\n${tex`\\ln a - \\ln b = \\ln \\frac{a}{b}`}, \nwe can rearrange as follows:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":411,"value":"tex.block`2 \\pi \\space \\Delta Y_M = \\ln\\left(1 + \\delta / e^{Y_0}\\right) - \\ln\\left(1 - \\delta e^{Y_0}\\right).`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":425,"value":"md`In this way, the ${tex`e^{Y_0}`} factor can be set as a uniform, and only ${tex`\\delta`} needs to be computed in the fragment shader. Also, the ${tex`\\ln(1 + x)`} form naturally lends itself to a robust expansion when ${tex`\\delta`} is small.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":877,"value":"projMercatorSrc = `\nconst float ONEOVERTWOPI = 0.15915493667125702;\nuniform float uExpY0;\nuniform float uLatErr; // Difference of clipping to map limit\n\nvec2 projMercator(vec2 dLonLat) {\n float tandlat = smallTan( 0.5 * (dLonLat[1] + uLatErr) );\n float p = tandlat * uExpY0;\n float q = tandlat / uExpY0;\n return vec2(dLonLat[0], log1plusX(-p) - log1plusX(q)) * ONEOVERTWOPI;\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1451,"value":"log1plusX_direct = `\nfloat log1plusX(float x) {\n return log( 1.0 + max(x, -0.999) ); // Avoid infinity at x = -1\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2048,"value":"md`When ${tex`x`} is small, information will be lost to roundoff error in the sum ${tex`(1 + x)`} inside the log. In this case, the Taylor expansion will better preserve the precision:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1456,"value":"tex.block`\\ln \\left( 1 + x \\right) \\approx x - \\frac{x^2}{2} + \\frac{x^3}{3} - \\frac{x^4}{4} + \\mathcal{O}(x^5)`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":881,"value":"log1plusX_optimal = `\nfloat log1plusX(float x) {\n return (abs(x) < 0.15)\n ? x * (1.0 - x * (0.5 - x / 3.0 + x * x / 4.0))\n : log( 1.0 + max(x, -0.999) );\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1454,"value":"log1plusX = ({\n direct: {\n src: log1plusX_direct,\n tex: tex`\\ln \\left( 1 + x \\right)`\n },\n optimal: {\n src: log1plusX_optimal,\n tex: tex`x - \\frac{x^2}{2} + \\frac{x^3}{3} - \\frac{x^4}{4}`\n }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2057,"value":"md`A somewhat surprising source of numerical problems is the GLSL tangent function, which has poor precision with small arguments. The Taylor series proves useful again:\n${tex.block`\\tan x \\approx x + \\frac{x^3}{3} + \\mathcal{O}(x^5)`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1482,"value":"smallTan_direct = `\nfloat smallTan(float x) {\n return tan(x);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":880,"value":"smallTan_optimal = `\nfloat smallTan(float x) {\n return (abs(x) < 0.1)\n ? x * (1.0 + x * x / 3.0)\n : tan(x);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1485,"value":"smallTan = ({\n direct: {\n src: smallTan_direct,\n tex: tex`\\tan x`\n },\n optimal: {\n src: smallTan_optimal,\n tex: tex`x + \\frac{x^3}{3}`\n }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2118,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":769,"value":"md`---\n## Appendices`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1885,"value":"md`### Miscellaneous constants`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":915,"value":"earthRadius = 6371","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":940,"value":"degrees = 180.0 / Math.PI","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":968,"value":"tanHalfFOV = Math.tan(verticalFieldOfView / 2 / degrees)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":970,"value":"verticalFieldOfView = 25","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":965,"value":"globeDisplayHeight = Math.min(width * 1.5, 600)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":60,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":900,"value":"md`### Rendering\nWe need to redraw when the map changes:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":794,"value":"cacheSize = {\n while (true) {\n let mapChanged = map.drawTiles(); // drawTiles does nothing if map is complete\n let numCachedTiles = tileCache.prune(map.tileDistance, 1.0);\n if (mapChanged) drawGlobe(camPosition);\n yield numCachedTiles;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1293,"value":"md`...and also when the camera position changes:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1287,"value":"drawGlobe(camPosition)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1296,"value":"md`The drawGlobe function uses [yawgl](https://www.npmjs.com/package/yawgl) to manage the tedious parts of the WebGL API.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":903,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1207,"value":"camMapPosition = {\n let center = lonLatToXY(camPosition);\n\n let hNorm = camPosition[2] / earthRadius;\n let dTanDPixel = (2 * tanHalfFOV) / globeDisplayHeight;\n let dMap = hNorm * dTanDPixel * getMercatorScale(camPosition[1]);\n let zoom = Math.round(-Math.log2(tileSize * dMap));\n\n map.setCenterZoom(center, zoom);\n\n let mapXY = [];\n map.toLocal(mapXY, center);\n\n return mapXY;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":993,"value":"function lonLatToXY(lonLat) {\n let x = 0.5 + (0.5 * lonLat[0]) / Math.PI;\n\n let expY = Math.tan(Math.PI / 4.0 + clipLatitude(lonLat[1]) / 2.0);\n let y = 0.5 - (0.5 / Math.PI) * Math.log(expY); // Note sign flip\n y = Math.min(Math.max(0.0, y), 1.0); // y doesn't wrap\n\n return [x, y];\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":998,"value":"getMercatorScale = latitude =>\n 1 / (2 * Math.PI * Math.cos(clipLatitude(latitude)))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1014,"value":"clipLatitude = {\n const maxMercatorLat = 2.0 * Math.atan(Math.exp(Math.PI)) - Math.PI / 2.0;\n return latitude =>\n Math.min(Math.max(-maxMercatorLat, latitude), maxMercatorLat);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":839,"value":"md`### WebGL stuff\nHere is where we put all the pieces of GLSL code together:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":886,"value":"fragmentShaderSrc = fragmentShaderHeader +\n glslFunctions.diffSqrt +\n calcSinCSrc +\n glslFunctions.latChange +\n xyToLonLatSrc +\n glslFunctions.smallTan +\n glslFunctions.log1plusX +\n projMercatorSrc +\n horizonTaperSrc +\n fragmentMainSrc","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":843,"value":"fragmentShaderHeader = `\n#extension GL_OES_standard_derivatives : enable\nprecision highp float;\nprecision highp sampler2D;\nvarying vec2 vRayParm;\n\nuniform sampler2D uTextureSampler;\nuniform vec2 uCamMapPos;\nuniform vec2 uMapScale;\n\nuniform float uLat0;\nuniform float uCosLat0;\nuniform float uSinLat0;\nuniform float uTanLat0;`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2136,"value":"horizonTaperSrc = `\nfloat horizonTaper(float p) {\n float p2 = p * p;\n float gamma = p2 * uHnorm * (2.0 + uHnorm);\n\n // sqrt(gamma) = tan(ray_angle) / tan(horizon)\n float horizonRatio = sqrt(gamma);\n float delta = 2.0 * fwidth(horizonRatio);\n return 1.0 - smoothstep(1.0 - delta, 1.0, horizonRatio);\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":824,"value":"vertexShaderSrc = `\nattribute vec4 aVertexPosition;\nuniform vec2 uMaxRay;\nvarying highp vec2 vRayParm;\n\nvoid main(void) {\n vRayParm = uMaxRay * aVertexPosition.xy;\n gl_Position = aVertexPosition;\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":658,"value":"yawgl = import(\"yawgl@0.1.2\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":832,"value":"md`### Map setup`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":786,"value":"map = tileFrame.init({\n getTile: tileCache.retrieve,\n context: mapContext,\n width: mapSize,\n height: mapSize,\n tileSize: tileSize,\n maxZoom: 22\n})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":784,"value":"tileCache = tileRack.initRasterCache(tileSize, tileURL)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2121,"value":"mapSize = {\n let screenSize = Math.max(screen.width, screen.height);\n let power = Math.ceil(Math.log2(screenSize));\n return 2 ** power;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":788,"value":"tileSize = 256","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":780,"value":"tileURL = (z, x, y) =>\n `https://stamen-tiles-${\n \"abc\"[Math.abs(x + y) % 3]\n }.a.ssl.fastly.net/toner-lite/${z}/${x}/${y}.png`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":767,"value":"tileFrame = import(\"tile-frame@1.0.0\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":778,"value":"tileRack = import(\"tile-rack@0.2.1\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W39","stdlib":"1","titling":"auto"}} 30:"$2f:body"

Optimizations

Direct implementation(click to toggle)

With optimizations(click to toggle)

Direct implementation
(click to toggle)

With optimizations
(click to toggle)