30:T707,{ const svg = d3 .select(DOM.svg(250, 120)) .style("font-style", "italic") .style("font-size", "12px"); svg .append("text") .attr("x", 50) .attr("y", 15) .text("Azimuth") .style("text-anchor", "middle"); svg .append("circle") .attr("cx", 50) .attr("cy", 70) .attr("r", 40) .style("stroke", "grey") .style("fill", "none"); svg .append("rect") .attr("x", 46) .attr("y", 66) .attr("width", 8) .attr("height", 8) .style("fill", "black"); const sunAzimuth = (azimuth * Math.PI) / 180; const dx = Math.cos(sunAzimuth - Math.PI / 2) * 40; const dy = Math.sin(sunAzimuth - Math.PI / 2) * 40; svg .append("line") .attr("x1", 50) .attr("x2", 50 + dx) .attr("y1", 70) .attr("y2", 70 + dy) .style("stroke", "black"); svg .append("circle") .attr("cx", 50 + dx) .attr("cy", 70 + dy) .attr("r", 5) .attr("fill", "gold"); const g = svg.append("g").attr("transform", "translate(150, 30)"); g.append("text") .attr("x", 40) .attr("y", -15) .text("Elevation") .style("text-anchor", "middle"); g.append("path") .attr("d", "M80 80 a 80 80 0 0 0 -80 -80") .style("fill", "none") .style("stroke", "grey"); const sunElev = (elevation * Math.PI) / 180; const dx2 = Math.cos(sunElev) * 80; const dy2 = Math.sin(sunElev) * 80; g.append("line") .attr("x1", 0) .attr("y1", 80) .attr("x2", dx2) .attr("y2", 80 - dy2) .style("stroke", "black"); g.append("circle") .attr("cx", dx2) .attr("cy", 80 - dy2) .attr("r", 5) .attr("fill", "gold"); g.append("rect") .attr("x", -4) .attr("y", 76) .attr("width", 8) .attr("height", 8) .style("fill", "black"); return svg.node(); }31:T5c5,hillshader = (demData, azimuth, elevation) => { var data = new Uint8ClampedArray(demData.length); for (var n = 0; n < demData.length; n += 4) { // value 0 is ocean, make it white if (demData[n] < 1) { data[n] = 255; data[n + 1] = 255; data[n + 2] = 255; data[n + 3] = 255; continue; } // x, y coordinates from imageData index const x = (n / 4) % width; const y = Math.floor(n / 4 / width); // image data index for surrounding pixels const top = getIndexForCoordinates(x, Math.max(0, y - 1)); const left = getIndexForCoordinates(Math.max(0, x - 1), y); const right = getIndexForCoordinates(Math.min(width - 1, x + 1), y); const bottom = getIndexForCoordinates(x, Math.min(height, y + 1)); // get slope values const topValue = demData[top]; const leftValue = demData[left]; const rightValue = demData[right]; const bottomValue = demData[bottom]; const { dzdx, dzdy } = getDerivatives( topValue, leftValue, bottomValue, rightValue ); const slope = getSlope(dzdx, dzdy); const aspect = getAspect(dzdx, dzdy); var L = getReflectance( aspect, slope, (Math.PI / 180) * azimuth, (Math.PI / 180) * elevation ); if (L < 0) L = 0; L = brighten(L); data[n] = parseInt(L * 255); data[n + 1] = parseInt(L * 255); data[n + 2] = parseInt(L * 255); data[n + 3] = 255; } return data; }32:[] 2e:{"title":"Hillshader","cells":[{"id":445,"value":"md`# Hillshader\n\nI came across Andy Woodruff’s [demonstration](https://github.com/awoodruff/canvas-shaded-relief) of hillshading the Hawaiian island of Moloka’i in the browser and thought I’d learn how it works by adapting his example to an Observable notebook. Much of the code is from Andy’s original page and his [earlier](http://andywoodruff.com/blog/shaded-relief-in-as3/) post about shaded relief in ActionScript.\n\nHere’s the final result: a textured terrain layer combined with a [hypsometric tint](https://en.wikipedia.org/wiki/Hypsometric_tints).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":506,"value":"html`\n\n

${drawMap(colorData)}

${drawHillshade(azimuth, elevation)}

Moloka’i

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":522,"value":"md`Hillshading algorithms have two important parameters: the *azimuth* of the sun, which controls the direction that light comes from, and its *elevation*, which controls how high it is in the sky. Normally, the azimuth is 315 degrees, which corresponds to light from the northwest, and the elevation is 45 degrees.\n\nTry changing those parameters.` ","pinned":false,"hidden":false,"mode":"ojs"},{"id":238,"value":"viewof azimuth = slider({\n min: 0,\n max: 360,\n value: 315,\n step: 1,\n description: \"Azimuth in degrees\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":231,"value":"viewof elevation = slider({\n min: 0,\n max: 90,\n value: 45,\n step: 1,\n description: \"Elevation in degrees\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":266,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":367,"value":"drawHillshade(azimuth, elevation)","pinned":false,"hidden":false,"mode":"ojs"},{"id":928,"value":"md`To better see the effects of the two parameters, here’s a grid of hillshades computed for each 45 degree increment of the azimuth and each of three elevation values: 30, 45 and 60 degrees.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":888,"value":"html`\n\n

\nElevation = 30°\n45°\n60°\n

\n${d3.range(0, 360, 45).map((d) => drawRow(d))}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":494,"value":"md`---\n## How it works\n\nThis simple hillshade algorithm works by calculating the local slope at each pixel using just the four pixels immediately above, below, left and right of it.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":472,"value":"html`Given the elevation ${tex`Z`} at location ${tex`(i, j)`}, the local partial derivatives in the ${tex`x`} and ${tex`y`} directions are approximated by:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":455,"value":"tex`dz/dx = Z(i+1, j) - Z(i-1, j)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":462,"value":"tex`dz/dy = Z(i, j+1) - Z(i, j-1)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":483,"value":"md`Here’s the code to do it. (I think there’s a faster way using a [convolution matrix](https://en.wikipedia.org/wiki/Kernel_(image_processing)).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1087,"value":"function getDerivatives(top, left, bottom, right) {\n const dzdx = right - left;\n const dzdy = top - bottom;\n\n return {\n dzdx,\n dzdy\n };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1119,"value":"md`Then, find the slope angle, which measures how steep the slope of the land is`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1097,"value":"tex`\\theta = \\tan^{-1} ( \\sqrt{(\\frac{dz}{dx})^2 + (\\frac{dz}{dy})^2}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1094,"value":"getSlope = (dzdx, dzdy, z = 0.2) =>\n Math.atan(z * Math.sqrt(dzdx ** 2 + dzdy ** 2)) // the z factor controls how exaggerated the peaks look","pinned":true,"hidden":false,"mode":"ojs"},{"id":1123,"value":"md`And the aspect, the direction that the land faces`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1129,"value":"tex`\\phi = \\text{atan2}( -\\frac{dz}{dy}, -\\frac{dz}{dx})`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1399,"value":"md`We want the angle that is exactly opposite the direction of the steepest slope, hence the negative signs. For example, if the steepest ascent is in the positive ${tex`x`} and ${tex`y`} directions (to the northeast) then the aspect should be to the southwest.\n\nIt’s important to use the \\`atan2\\` function instead of the normal arctangent because \\`atan2\\` returns values in ${tex`[-\\pi, \\pi]`} as opposed to ${tex`[-\\frac{\\pi}{2}, \\frac{\\pi}{2}]`}.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1126,"value":"getAspect = (dzdx, dzdy) => {\n return Math.atan2(-dzdy, -dzdx);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":646,"value":"md`Finally, calculate a luminance ${tex`L`} for each pixel using a [Lambertian reflectance model](https://en.wikipedia.org/wiki/Lambertian_reflectance), which means the light scattered from the surface is a function of the dot product between the incident light vector and the surface normal.\n\nThe formula is an application of the [spherical law of cosines](https://en.wikipedia.org/wiki/Spherical_law_of_cosines), where ${tex`\\alpha`} is the sun’s azimuth and ${tex`\\beta`} is the elevation. Note that we need to subtract the aspect from 90 degrees to convert it to an azimuthal angle (measured clockwise from due north).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1010,"value":"tex`L = \\cos(90° - \\phi - \\alpha) \\sin( \\theta ) \\sin(90° - \\beta) + \\cos(\\theta) \\cos(90° - \\beta)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1049,"value":"function getReflectance(aspect, slope, sunAzimuth, sunElevation) {\n return (\n Math.cos(Math.PI * 0.5 - aspect - sunAzimuth) *\n Math.sin(slope) *\n Math.sin(Math.PI * 0.5 - sunElevation) +\n Math.cos(slope) * Math.cos(Math.PI * 0.5 - sunElevation)\n );\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":285,"value":"md`---\n## Hypsometric tint`","pinned":false,"hidden":false,"mode":"ojs"},{"id":717,"value":"md`Create a separate layer for the elevation tint by mapping the elevation value from the DEM to a color. Here’s the gradient.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":738,"value":"{\n const canvas = DOM.canvas(255 * 2, 15);\n const ctx = canvas.getContext('2d');\n \n for (var i = 0; i < 255; i++) {\n ctx.beginPath();\n ctx.rect(i * 2, 0, 2, 15);\n ctx.fillStyle = colorScale(i);\n ctx.fill();\n }\n \n return canvas;\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":1477,"value":"md`And the image.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":370,"value":"drawMap(colorData)","pinned":false,"hidden":false,"mode":"ojs"},{"id":771,"value":"md`In the combined map at the top of the page, the hillshade layer is overlaid and multiplied with the color tint layer (CSS makes this surprisingly easy!). There's also a slight blur.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1061,"value":"md`---\n\n## Further reading\n\n* [Hillshading, useful fun with digital elevation models](https://scientiaplusconscientia.wordpress.com/2014/12/16/hillshading-digital-elevation-model/)\n* [Hill Shading and the Reflectance Map](http://people.csail.mit.edu/bkph/papers/Hill-Shading.pdf)\n* [How HillShade works](https://pro.arcgis.com/en/pro-app/tool-reference/3d-analyst/how-hillshade-works.htm)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":439,"value":"md`---\n## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":792,"value":"md`Here’s what the DEM looks like. It’s a PNG, so it can only represent 255 elevation values. Black regions are the lowest, white are the highest.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":156,"value":"drawMap(demData);","pinned":false,"hidden":false,"mode":"ojs"},{"id":998,"value":"drawHillshade = (azimuth, elevation) => {\n return drawMap(hillshader(demData, azimuth, elevation));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1483,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":990,"value":"drawRow = (azimuth) => {\n return html`

\n ${\n azimuth == 0 ? \"Azi- \\n muth \\n = \" : \"\"\n }${azimuth}°\n ${drawSmallMultiple(azimuth, 30)}\n ${drawSmallMultiple(azimuth, 45)}\n ${drawSmallMultiple(azimuth, 60)}\n

\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":993,"value":"drawSmallMultiple = (azimuth, elevation) => {\n const canvas = drawHillshade(azimuth, elevation);\n canvas.style.width = `${width * 0.3}px`;\n canvas.style.height = `${height * 0.3}px`;\n return canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":827,"value":"drawMap = (data) => {\n const canvas = DOM.canvas(width, height);\n const ctx = canvas.getContext(\"2d\");\n\n const imgData = new ImageData(data, width, height);\n\n ctx.putImageData(imgData, 0, 0, 0, 0, width, height);\n\n return canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":804,"value":"demData = new Promise((resolve, reject) => {\n const image = new Image();\n image.crossOrigin = \"anonymous\";\n image.onerror = reject;\n image.onload = () => {\n const canvas = DOM.canvas(width, height);\n const ctx = canvas.getContext(\"2d\");\n\n ctx.drawImage(image, 0, 0, image.width, image.height, 0, 0, width, height);\n const data = ctx.getImageData(0, 0, width, height).data;\n\n resolve(data);\n };\n image.src = molokai;\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":808,"value":"colorScale = d3\n .scaleLinear()\n .domain([3, 53, 103, 153, 203, 255])\n .range([\"#004000\", \"#679167\", \"#81b279\", \"#bfdfa8\", \"#d0b8aa\", \"#ffffff\"])","pinned":false,"hidden":false,"mode":"ojs"},{"id":800,"value":"colorData = {\n var data = new Uint8ClampedArray(demData.length);\n\n for (var n = 0; n < demData.length; n += 4) {\n // value 0 is ocean\n\n if (demData[n] < 1) {\n data[n] = blue.r;\n data[n + 1] = blue.g;\n data[n + 2] = blue.b;\n data[n + 3] = 255;\n continue;\n }\n\n const color = d3.color(colorScale(demData[n]));\n\n data[n] = color.r;\n data[n + 1] = color.g;\n data[n + 2] = color.b;\n data[n + 3] = 255;\n }\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":783,"value":"molokai = 'https://s3-us-west-1.amazonaws.com/sahilchinoy/molokai.png'","pinned":false,"hidden":false,"mode":"ojs"},{"id":780,"value":"blue = {\n return { r: 150, g: 177, b: 228 };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":764,"value":"brighten = (n) => Math.sqrt(n * 0.8 + 0.2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":76,"value":"getIndexForCoordinates = (x, y) => width * y * 4 + 4 * x","pinned":false,"hidden":false,"mode":"ojs"},{"id":18,"value":"width = 960","pinned":false,"hidden":false,"mode":"ojs"},{"id":21,"value":"height = 351","pinned":false,"hidden":false,"mode":"ojs"},{"id":228,"value":"import {slider} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":267,"value":"d3 = require(\"d3@7\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}