2a:null 2b:["$","title",null,{"children":"Mapbox hillshade and satellite map blending | Armand Emamdjomeh | Observable"}] 31:T5bc,md`I've been inspired lately by Sahil Chinoy's [hillshading](https://observablehq.com/@sahilchinoy/a-faster-hillshader) [explanations](https://observablehq.com/@sahilchinoy/hillshader), which are in turn based off of Andy Woodruff's [examples](https://github.com/awoodruff/canvas-shaded-relief). Sahil's workflow made me think, if we could do this with a static PNG, maybe we could do this with map tiles? Andy already has an explainer [using Leaflet and translating RGB-terrain tiles to a hillshade layer](https://observablehq.com/@awoodruff/diy-hillshade), but I've been experimenting with the idea of blending satellite and hillshade layers in the browser, using Mapbox [satellite](https://www.mapbox.com/maps/satellite/) and [terrain-rgb](https://docs.mapbox.com/help/troubleshooting/access-elevation-data/#mapbox-terrain-rgb) data. This builds on a workflow for [making static maps using vector tile data](https://observablehq.com/@aboutaaron/mapbox-map-maker) that my colleague [Aaron Williams](https://observablehq.com/@aboutaaron) and I have been working on. This is similar to a more traditional QGIS/Photoshop workflow we use at the Post to make [hillshaded maps](https://www.washingtonpost.com/graphics/2019/world/india-pakistan-conflict/) with a [satellite base layer](https://www.washingtonpost.com/graphics/2018/national/california-wildfires-maps/?utm_term=.fa2697d65ebd), and give you more of an impression of the terrain that is being hillshaded. `32:T7af,map = { const w = width; const h = width * 2/3; const canvas = DOM.canvas(width, width * 2/3) const context = canvas.getContext('2d') const path = d3.geoPath(projection, context) // This could be used to build a matrix if we need it. const xMax = d3.max(satellite_tiles.map(d => d.x)); const xMin = d3.min(satellite_tiles.map(d => d.x)); const yMax = d3.max(satellite_tiles.map(d => d.y)); const yMin = d3.min(satellite_tiles.map(d => d.y)); const xRows = xMax - xMin + 1; const yRows = yMax - yMin + 1; // Drawing the satellite raster layers satellite_tiles.map(d => { const dx = (d.x + tiles.translate[0]) * tiles.scale; const dy = (d.y + tiles.translate[1]) * tiles.scale; context.drawImage(d.data, dx, dy, tileDim, tileDim) }); const rasterImageData = context.getImageData(0, 0, w, h).data; const hillshadeData = hillshader(elevationData, azimuth, elevation, z, w, h) const compositeImageData = blendImage(rasterImageData, hillshadeData); const compositeImage = new ImageData(compositeImageData, w, h); context.putImageData(compositeImage, 0, 0, 0, 0, w, h); // Drawing road vectors context.save(); context.lineWidth = 1.0; context.strokeStyle = "#666666"; context.globalAlpha = 0.5; vector_tiles.map(d => { context.beginPath(); path(filter(d.data.road, d => isLayer(d, 'motorway'))); context.stroke(); }); // Drawing water context.globalAlpha = 1.0; context.fillStyle = "#acc8cc"; // context.fillStyle = "#ffffff"; vector_tiles.map(d => { context.beginPath(); path(filter(d.data.water, d => isLayer(d, 'water'))); context.filter = "blur(2px)"; context.fill("evenodd"); }); context.filter = "none"; const mapboxLogo = new Image(); mapboxLogo.onload = function() { context.drawImage(mapboxLogo, 0, height - 30) } mapboxLogo.src = `data:image/svg+xml;base64,${mbxbase64}` return context.canvas; }33:T79e,blendImage = (layer1, layer2) => { const compositeImageData = new Uint8ClampedArray(layer1.length); for (let n = 0, l = layer1.length; n < l; n += 4) { let srcR = layer1[n]; let srcG = layer1[n + 1]; let srcB = layer1[n + 2]; // Alpha is treated as a decimal 0 to 1 instead of 0 to 255 // let srcA = layer1[n + 3] / 255; let srcA = satelliteAlpha; let deltR = layer2[n]; let deltG = layer2[n + 1]; let deltB = layer2[n + 2]; // Alpha is treated as a decimal 0 to 1 instead of 0 to 255 // let deltA = layer2[n + 3] / 255; let deltA = hillshadeAlpha; let newR, newG, newB; switch(blendMethod) { case 'multiply': newR = srcR * deltR / 255; newG = srcG * deltG / 255; newB = srcB * deltB / 255; break; case 'screen': newR = 255 - ( ( (255 - srcR) * (255 - deltR) ) / 255); newG = 255 - ( ( (255 - srcG) * (255 - deltG) ) / 255); newB = 255 - ( ( (255 - srcB) * (255 - deltB) ) / 255); break; case 'color dodge': newR = srcR + deltR; newG = srcG + deltG; newB = srcB + deltB; break; case 'overlay' : newR = deltR < 128 ? (2 * srcR * deltR / 255) : (255 - ( ( 2 * (255 - srcR) * (255 - deltR) ) / 255)); newG = deltG < 128 ? (2 * srcG * deltG / 255) : (255 - ( ( 2 * (255 - srcG) * (255 - deltG) ) / 255)); newB = deltB < 128 ? (2 * srcB * deltB / 255) : (255 - ( ( 2 * (255 - srcB) * (255 - deltB) ) / 255)); break; } // Composite layers based on alpha newR = newR * srcA + deltR * deltA; newG = newG * srcA + deltG * deltA; newB = newB * srcA + deltB * deltA; let newA = 255; compositeImageData[n] = newR; compositeImageData[n + 1] = newG; compositeImageData[n + 2] = newB; compositeImageData[n + 3] = newA; } return compositeImageData; }34:[] 2e:{"id":"5eccb672c04640f4","slug":"mapbox-hillshade-and-satellite-map-blending","title":"Mapbox hillshade and satellite map blending","thumbnail":"883b66e37be7f4281cf8f8c0e6982ac51367f4639705558aa5ce543c8a108d85","access_level":"public","create_time":"$D2019-03-04T19:20:02.062Z","update_time":"$D2024-06-13T14:15:10.172Z","view_time":null,"first_public_time":"$D2019-03-21T18:20:23.869Z","fork":false,"stars":59,"starred":false,"owner":{"type":"team","id":"775f5069ff441bf9","login":"emamd","name":"Armand Emamdjomeh","avatar_url":"https://avatars-next.observableusercontent.com/7a932cbd387d1f561d8d1c405bf9e61ef1b6a679e1fa54aaa262a54380ab2f9b","tier":"starter_2024"},"version":1776,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"aac33ad2144eec5b","login":"emamd","name":"Armand Emamdjomeh","avatar_url":"https://avatars-next.observableusercontent.com/7a932cbd387d1f561d8d1c405bf9e61ef1b6a679e1fa54aaa262a54380ab2f9b","tier":null}],"members":[{"type":"individual","id":"aac33ad2144eec5b","login":"emamd","name":"Armand Emamdjomeh","avatar_url":"https://avatars-next.observableusercontent.com/7a932cbd387d1f561d8d1c405bf9e61ef1b6a679e1fa54aaa262a54380ab2f9b","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":28,"collection":null,"type":"notebook","body":{"title":"Mapbox hillshade and satellite map blending","cells":[{"id":0,"value":"md`# Mapbox hillshade and satellite map blending`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1096,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1276,"value":"md`## Blend settings`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1136,"value":"md`The map works by requesting a satellite raster, vector, and terrain-rgb tile layers and combining them into one map. You can see the results below. Changing the different blend modes changes how the satellite and hillshade layers affect one another, so play around with them. I typically like \"multiply\" best since it tends to just punch up the terrain a bit.\n\nBy changing the opacity of each layer, you can affect how much of each one comes through in the final image.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1327,"value":"md`Change the opacity of the satellite and hillshade layers, and see how the combined image reacts.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1248,"value":"md`## Location`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1417,"value":"md`You can use the map below to center your view. Zoom is set further down but I'd like to eventually incorporate the map being set by the locator bounds as well.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1646,"value":"viewof centerpoint = {\n const container = html`

`;\n yield container;\n\n let map = new mapboxgl.Map({\n container,\n center: [-118.2437, 34.0522],\n zoom: 8.5,\n style: \"mapbox://styles/mapbox/light-v9\",\n trackResize: false\n });\n\n map.on(\"moveend\", (e) => {\n container.value = {\n center: map.getCenter(),\n zoom: map.getZoom()\n };\n container.dispatchEvent(new CustomEvent(\"input\"));\n });\n\n container.value = {\n center: map.getCenter(),\n zoom: map.getZoom()\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"viewof zoomscale = slider({\n min: 50000,\n max: 1000000,\n value: 500000,\n step: 1000,\n title: \"Zoom Scale\", \n description: \"Use slider to adjust values\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":82,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":998,"value":"viewof blendMethod = select({\n title: \"Blend method\",\n description: \"Please choose how you would like to blend layers. Multiply or overlay are recommended.\",\n options: [\"screen\", \"multiply\", \"overlay\", \"color dodge\"],\n value: \"multiply\"\n})\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1007,"value":"viewof satelliteAlpha = slider({\n min: 0,\n max: 1,\n value: 0.7,\n step: 0.1,\n description: \"Satellite layer opacity\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1008,"value":"viewof hillshadeAlpha = slider({\n min: 0,\n max: 1,\n value: 0.6, \n step: 0.1,\n description: \"Hillshade layer alpha\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1281,"value":"md`## Light settings`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1133,"value":"md`The below settings will change the direction the \"light\" is shining on the scene, as well as the angle and strength of shadows in the scene. Generally, the defaults are fine but it's fun to play around to see their effect.` ","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"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":426,"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":428,"value":"viewof z = slider({\n min: 0,\n max: 2,\n value: 1.0, \n step: 0.1,\n description: \"Z-factor\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1242,"value":"md`There are still a couple issues here that I'll call out as we go along, mainly some striations in the hillshading that look like contour lines. Shout out and thanks to [@noodlecube](https://observablehq.com/@noodlecube) for fixing a bug with one row of black pixels the canvas.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":387,"value":"md`## How this works\n\n### Tile layers`","pinned":false,"hidden":false,"mode":"ojs"},{"id":145,"value":"md`This uses three different types of tiles - satellite raster, terrain-rgb heightmaps, and vector data tiles. \n\nFor a satellite base layer raster, all we need is the plain image.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":77,"value":"getSatelliteData = function(url) {\n return new Promise((resolve, reject) => { \n const image = new Image();\n image.crossOrigin = \"anonymous\";\n image.onerror = reject;\n\n image.onload = () => {\n resolve(image);\n }\n image.src = url\n })\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":76,"value":"satellite_tiles = Promise.all(tile().map(async t => {\n t.data = await getSatelliteData(`https://api.mapbox.com/v4/mapbox.satellite/${t.z}/${t.x}/${t.y}.pngraw?access_token=pk.eyJ1IjoiZW1hbWQtai1hcm1hbmQiLCJhIjoiY2p0aTBwY3drMGt0ajQ0cGYyMWozZTUyMCJ9.VcYxgVc6cJMaYNC3w-KUyA`)\n return t\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1215,"value":"satpreview = {\n const canvas = DOM.canvas(width, height / 2);\n const context = canvas.getContext('2d');\n \n satellite_tiles.map(d => {\n const dx = (d.x + tiles.translate[0]) * tiles.scale;\n const dy = (d.y + tiles.translate[1]) * tiles.scale;\n context.drawImage(d.data, dx, dy, tileDim, tileDim) \n });\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":149,"value":"md`Hillshade layers need the digital elevation model parsed out as an array. Currently, because I don't know how to apply these algorithms onto an RGB-scaled terrain, I convert to a grayscale DEM. This winds up in a large loss of detail, since grayscale is limited to 255 values and RGB has... many times more. Something to figure out in the future.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"getDemData = function(url, transform=true) {\n return new Promise((resolve, reject) => { \n const image = new Image();\n image.crossOrigin = \"anonymous\";\n image.onerror = reject;\n image.src = url\n\n image.onload = () => {\n const canvas = DOM.canvas(tileSize, tileSize);\n const ctx = canvas.getContext('2d');\n\n ctx.drawImage(image, 0, 0, image.width, image.height, 0, 0, tileDim, tileDim);\n const rawData = ctx.getImageData(0, 0, image.width, image.height).data;\n \n let data\n if (transform === true) {\n data = convertToGrayScale(rawData);\n } else {\n data = rawData;\n }\n\n resolve(data);\n }\n })\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1209,"value":"md`This is what the terrain-RGB tiles look like without any transformations.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1202,"value":"terrainrgb_tiles = Promise.all(tile().map(async t => {\n t.data = await getDemData(`https://api.mapbox.com/v4/mapbox.terrain-rgb/${t.z}/${t.x}/${t.y}.pngraw?access_token=pk.eyJ1IjoiZW1hbWQtai1hcm1hbmQiLCJhIjoiY2p0aTBwY3drMGt0ajQ0cGYyMWozZTUyMCJ9.VcYxgVc6cJMaYNC3w-KUyA`, false)\n \n return t\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1204,"value":"rgbpreview = {\n const canvas = DOM.canvas(width, height / 2);\n const context = canvas.getContext('2d');\n \n terrainrgb_tiles.map(d => {\n const dx = (d.x + tiles.translate[0]) * tiles.scale;\n const dy = (d.y + tiles.translate[1]) * tiles.scale;\n const imageData = new ImageData(d.data, tileSize, tileSize);\n context.putImageData(imageData, dx, dy, 0, 0, tileDim, tileDim) \n });\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1450,"value":"getElevation = (index, demData) => {\n if (!demData) return 0;\n if (index < 0 || demData[index] === undefined) return undefined;\n let r = demData[index],\n g = demData[index + 1],\n b = demData[index + 2]\n // Converting Mapbox terrain-rgb tiles\n // -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)\n return -10000 + ((r * 256 * 256 + g * 256 + b) * 0.1);\n // Looking at the mapbox hillshade shader, it looks like they're using this to determine elevation?\n // https://github.com/mapbox/mapbox-gl-js/blob/a3deabf2cc1df16809112fb7b9e37fdd0615a429/src/shaders/hillshade_prepare.fragment.glsl#L12-L13\n // vec4 data = texture2D(u_image, coord) * 255.0;\n // return (r + g * 256.0 + b * 256.0 * 256.0) / 4.0;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1440,"value":"elevationData = {\n const canvas = DOM.canvas(width, height);\n const context = canvas.getContext('2d');\n\n terrainrgb_tiles.map(d => {\n const dx = (d.x + tiles.translate[0]) * tiles.scale;\n const dy = (d.y + tiles.translate[1]) * tiles.scale;\n const imageData = new ImageData(d.data, tileSize, tileSize);\n context.putImageData(imageData, dx, dy, 0, 0, tileSize, tileSize) \n });\n \n const demImageData = context.getImageData(0, 0, width, height);\n const demData = demImageData.data;\n const elevationValues = [];\n \n for (let y = 0; y < height; y++) {\n for (let x = 0; x < width; x++) {\n // x + y*width is the position for point (x, y) in a flat array of pixel data\n // multiply by 4 for imageData position as each pixel is 4 (rgba) values\n elevationValues[x + y * width] = getElevation(4 * (x + y * width), demData);\n }\n }\n return elevationValues;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":31,"value":"hillshade_tiles = Promise.all(tile().map(async t => {\n t.data = await getDemData(`https://api.mapbox.com/v4/mapbox.terrain-rgb/${t.z}/${t.x}/${t.y}.pngraw?access_token=pk.eyJ1IjoiZW1hbWQtai1hcm1hbmQiLCJhIjoiY2p0aTBwY3drMGt0ajQ0cGYyMWozZTUyMCJ9.VcYxgVc6cJMaYNC3w-KUyA`)\n \n return t\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1152,"value":"md`This is what the terrain-rgb data looks like parsed into a grayscale image.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1021,"value":"dempreview = {\n const canvas = DOM.canvas(width, height / 2);\n const context = canvas.getContext('2d');\n \n hillshade_tiles.map(d => {\n const dx = (d.x + tiles.translate[0]) * tiles.scale;\n const dy = (d.y + tiles.translate[1]) * tiles.scale;\n const imageData = new ImageData(d.data, tileSize, tileSize);\n context.putImageData(imageData, dx, dy, 0, 0, tileDim, tileDim) \n });\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":151,"value":"md`And here is the vector data that goes on top. Features are parsed out of their protobuf format and converted into geojson. What's useful about this approach is all the data that we might need is in the tile data, and we can filter only what we want to draw onto the map.` ","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"vector_tiles = Promise.all(tile().map(async t => {\n let data = new vt.VectorTile(new Protobuf(await d3.buffer(`https://a.tiles.mapbox.com/v4/mapbox.mapbox-streets-v7/${t.z}/${t.x}/${t.y}.mvt?access_token=pk.eyJ1IjoiZW1hbWQtai1hcm1hbmQiLCJhIjoiY2p0aTBwY3drMGt0ajQ0cGYyMWozZTUyMCJ9.VcYxgVc6cJMaYNC3w-KUyA`))).layers\n \n t.data = Object.keys(data).reduce((initialData, key) => {\n initialData[key] = {\n type: \"FeatureCollection\",\n features: [...Array(data[key].length).keys()].map((item, i) => {\n let feature = data[key].feature(i).toGeoJSON(t.x, t.y, t.z)\n \n if (key === 'water') {\n feature.properties['class'] = 'water'\n }\n \n return feature\n })\n }\n\n return initialData\n }, {})\n \n return t;\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":397,"value":"md`### Hillshader functions`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1219,"value":"md`The hillshading functions are largely adapted from Sahil's notebook, with some small tweaks. By merging the tiles into one image before calculating the hillshade, we avoid any issues with having to detect tile edges.\n\nNote the striations that emerge in flatter areas like the San Fernando valley in the top-left, I assume this is due to subtle elevation changes and looking for ways to smooth those out. Maybe also something to do with the RGB-grayscale conversion.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":842,"value":"hillshade = {\n const canvas = DOM.canvas(width, height / 2);\n const context = canvas.getContext('2d');\n\n const hillshadeData = hillshader(elevationData, azimuth, elevation, z, width, height)\n\n const imageData = new ImageData(hillshadeData, width, height);\n context.putImageData(imageData, 0, 0, 0, 0, width, height)\n\n return context.canvas\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":910,"value":"md`Here's what a hillshaded individual tile of the terrain-rgb data would look like.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":687,"value":"tinymap = {\n const canvas = DOM.canvas(tileDim, tileDim);\n const context = canvas.getContext('2d')\n context.fillColor = '#ffffff';\n context.fillRect(0, 0, tileDim, tileDim);\n \n const tinyData = terrainrgb_tiles[8].data\n\n const tinyElevationValues = []\n for (let y = 0; y < tileSize; y++) {\n for (let x = 0; x < tileSize; x++) {\n // x + y*width is the position for point (x, y) in a flat array of pixel data\n // multiply by 4 for imageData position as each pixel is 4 (rgba) values\n tinyElevationValues[x + y * tileSize] = getElevation(4 * (x + y * tileSize), tinyData);\n }\n }\n \n const hillshadeData = hillshader(tinyElevationValues, azimuth, elevation, z, tileSize, tileSize);\n const imageData = new ImageData(hillshadeData, tileSize, tileSize);\n\n context.putImageData(imageData, 0, 0, 0, 0, tileSize, tileSize) \n return context.canvas\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1460,"value":"hillshader = (demData, azimuth, elevation, z=.2, canvwidth=width, canvheight=height) => {\n\n const {a1, a2, a3} = getConstants(azimuth, elevation);\n let data = new Uint8ClampedArray(demData.length * 4);\n \n for (var n = 0, l = demData.length; n < l; n++) {\n const idx = n * 4; \n \n // X, Y coordinates from imageData index\n const x = n % canvwidth;\n const y = Math.floor(n / canvwidth);\n \n const sm = getSubmatrix(x, y, demData, canvwidth, canvheight);\n \n const { dzdx, dzdy } = getDerivativesHorn(sm, z);\n \n let luminance = (a1 - a2 * dzdx - a3 * dzdy) / Math.sqrt(1 + dzdx ** 2 + dzdy ** 2);\n \n // Luminance cannot be less than zero.\n if (luminance < 0) {\n luminance = 0;\n }\n \n luminance = brighten(luminance);\n \n data[idx] = parseInt(luminance * 255);\n data[idx+1] = parseInt(luminance * 255);\n data[idx+2] = parseInt(luminance * 255);\n data[idx+3] = 255;\n }\n \n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":409,"value":"// Create a submatrix of the pixel values surrounding a given X,Y point in a DEM data set.\ngetSubmatrix = (x, y, demData, w, h) => {\n const left = Math.max(0, x - 1);\n const right = Math.min(w - 1, x + 1);\n const top = Math.max(0, y - 1);\n const bottom = Math.min(h - 1, y + 1);\n\n const ix1 = getIndexForCoordinates(left, top, w);\n const ix2 = getIndexForCoordinates(x, top, w);\n const ix3 = getIndexForCoordinates(right, top, w);\n const ix4 = getIndexForCoordinates(left, y, w);\n const ix6 = getIndexForCoordinates(right, y, w);\n const ix7 = getIndexForCoordinates(left, bottom, w);\n const ix8 = getIndexForCoordinates(x, bottom, w);\n const ix9 = getIndexForCoordinates(right, bottom, w);\n\n const sm = {};\n\n sm.z1 = demData[ix1];\n sm.z2 = demData[ix2];\n sm.z3 = demData[ix3];\n sm.z4 = demData[ix4];\n sm.z6 = demData[ix6];\n sm.z7 = demData[ix7];\n sm.z8 = demData[ix8];\n sm.z9 = demData[ix9];\n \n return sm;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":399,"value":"getConstants = (azimuth, elevation) => {\n // convert to radians \n const alpha = Math.PI / 180 * (azimuth - 90);\n const beta = Math.PI / 180 * elevation; \n \n const a1 = Math.sin(beta);\n const a2 = Math.cos(beta) * Math.sin(alpha);\n const a3 = Math.cos(beta) * Math.cos(alpha);\n \n return {\n a1, a2, a3\n }\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":406,"value":"getIndexForCoordinates = (x, y, w) => w * y + x","pinned":false,"hidden":false,"mode":"ojs"},{"id":414,"value":"brighten = n => Math.sqrt(n * .8 + .2);","pinned":false,"hidden":false,"mode":"ojs"},{"id":1163,"value":"md`Hillshading algorithms adapted from [Sahil](https://observablehq.com/@sahilchinoy/a-faster-hillshader) and [Andy's](https://observablehq.com/@awoodruff/diy-hillshade) work.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":421,"value":"function getDerivativesHorn(sm, z) {\n const dzdx = z * ((sm.z3 + 2 * sm.z6 + sm.z9) - (sm.z1 + 2 * sm.z4 + sm.z7)) / px_size; \n const dzdy = z * ((sm.z7 + 2 * sm.z8 + sm.z9) - (sm.z1 + 2 * sm.z2 + sm.z3)) / px_size; \n \n return {\n dzdx,\n dzdy\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1030,"value":"md`### Image compositing`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1225,"value":"md`The image compositing works by taking two layers and blending each pixel according to the selected blend method. The idea and the blending formulas are from [this demo](http://bseth99.github.io/projects/canvas/3-html5-canvas-layers-blending.html).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1035,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":392,"value":"md`## Utilities`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1194,"value":"convertToGrayScale = (rgbdata) => {\n const grayscaledata = new Uint8ClampedArray(rgbdata.length);\n\n // Convert RGB-terrain values into a proportional greyscale value\n for (var n = 0, l = rgbdata.length; n < l; n += 4) {\n const height = -10000 + ((rgbdata[n] * 256 * 256 + rgbdata[n + 1] * 256 + rgbdata[n + 2]) * 0.1);\n const pct = height / maxHeight;\n const rgbVal = 255 * pct \n\n grayscaledata[n] = rgbVal;\n grayscaledata[n + 1] = rgbVal;\n grayscaledata[n + 2] = rgbVal;\n grayscaledata[n + 3] = rgbdata[n + 3]\n } \n return grayscaledata;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"getLayers = key => {\n const flatten = (arr, depth = 1) => arr.reduce((a, v) => a.concat(depth > 1 && Array.isArray(v) ? flatten(v, depth - 1) : v), [])\n \n return [ ...new Set(\n flatten(\n tiles.map(d => {\n if (d.data[key]) {\n return d.data[key].features.map(d => d.properties.class)\n }\n })\n )\n )\n ]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"filter = function(d, test) {\n if (!d || !d.features) { return }\n return { type: \"FeatureCollection\", features: d.features.filter(test) };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":165,"value":"isLayer = (feature, layer) => feature.properties.class === layer","pinned":false,"hidden":false,"mode":"ojs"},{"id":522,"value":"tileDim = Math.ceil(256 * scale);","pinned":false,"hidden":false,"mode":"ojs"},{"id":651,"value":"tileSize = 256","pinned":false,"hidden":false,"mode":"ojs"},{"id":1525,"value":"zoomLevel = tiles[0].z","pinned":false,"hidden":false,"mode":"ojs"},{"id":1638,"value":"exaggeration = zoomLevel < 2.0 ? 0.4 : zoomLevel < 4.5 ? 0.35 : 0.3;","pinned":false,"hidden":false,"mode":"ojs"},{"id":1527,"value":"px_size = 8 * 40075016.6855785 / (512 * Math.pow(2, zoomLevel))\n// px_size = Math.pow(2.0, (zoomLevel - 14.0) * exaggeration + 19.2562 - zoomLevel)","pinned":false,"hidden":false,"mode":"ojs"},{"id":516,"value":"scale = tiles.scale / 256;\n ","pinned":false,"hidden":false,"mode":"ojs"},{"id":582,"value":"maxHeight = 2000","pinned":false,"hidden":false,"mode":"ojs"},{"id":179,"value":"tiles = tile()","pinned":false,"hidden":false,"mode":"ojs"},{"id":23,"value":"tile = d3.tile()\n .size([width, height])\n .scale(projection.scale() * 2 * Math.PI)\n .translate(projection([0,0]))\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"projection = {\n let lng = centerpoint.center.lng;\n let lat = centerpoint.center.lat;\n \n return d3.geoMercator()\n .center([lng, lat])\n .scale(zoomscale / (2 * Math.PI))\n .translate([width / 2, height / 2])\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":812,"value":"width = 960","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"height = width","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"json = async url => {\n const { features } = await d3.json(url);\n return features;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":390,"value":"md`## Dependencies`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"d3 = require(\"d3-geo@1\", \"d3-selection@1\", \"d3-fetch@1\", \"d3-tile@0.0\", \"d3-scale@2.1.2\", \"d3-array@2\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":1649,"value":"mapboxgl = {\n const gl = await require(\"mapbox-gl@0.53\");\n if (!gl.accessToken) {\n gl.accessToken =\n \"pk.eyJ1IjoiZW1hbWQtai1hcm1hbmQiLCJhIjoiY2p0aTBwY3drMGt0ajQ0cGYyMWozZTUyMCJ9.VcYxgVc6cJMaYNC3w-KUyA\";\n const href = await require.resolve(\"mapbox-gl@0.53/dist/mapbox-gl.css\");\n document.head.appendChild(html``);\n }\n return gl;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"vt = require('https://bundle.run/@mapbox/vector-tile@1.3.1')","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"Protobuf = require('https://bundle.run/pbf@3.1.0')","pinned":false,"hidden":false,"mode":"ojs"},{"id":112,"value":"import {slider, text, select, coordinates, worldMapCoordinates, usaMapCoordinates} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1053,"value":"import {serialize} from \"@mbostock/saving-svg\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1296,"value":"import {mbxLogo as mbxLogo} from '@aboutaaron/mapbox-map-maker'\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1300,"value":"mbxbase64 = window.btoa(new XMLSerializer().serializeToString(mbxLogo))","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}} 2f:"$2e:body"