31:T60b,rain = ({ DZX, DZY, invalidation }) => { const random = d3.randomNormal(0, 200); const nodes = Array.from({ length: Math.min(w * 10, npoints * 2) }, () => ({ x: w / 2 + random(), y: h / 2 + random() })); const context = DOM.context2d(w, h, 1); const path = d3.geoPath(d3.geoIdentity()).context(context); function gradient(x, y) { x = Math.floor(x * w); y = Math.floor(y * h); const k = x + w * y; return [DZX[k] || 0, DZY[k] || 0]; } const f = -0.06; const simulation = d3 .forceSimulation(nodes) .alphaDecay(0) .force("wind", () => { nodes.forEach((d) => { if ( !d.age || Math.random() > d.age || d.x < 0 || d.x >= w - 1 || d.y < 0 || d.y >= h - 1 ) { d.x = w * Math.random(); d.y = h * Math.random(); d.vx = 0; d.vy = 0; d.age = 1; } else { const [vx, vy] = gradient(d.x / w, d.y / h); d.vx += f * vx; d.vy += f * vy; d.age *= 0.9999; } }); }) .on("tick", () => { fade(context, 0.1); context.beginPath(); context.strokeStyle = "steelblue"; context.lineWidth = 3; for (const { x, y, vx, vy } of nodes) { context.moveTo(x + vx / 2, y + vy / 2); context.lineTo(x - vx, y - vy); } context.stroke(); }); invalidation.then(() => simulation.stop()); context.canvas.style = `max-width:${w / 2}px`; return context.canvas; }32:T64d,streaklines = (grid, DZX, DZY, { avoidance = 0.01, seed = 42 } = {}) => { function altitude(x, y) { const k = Math.floor(x) + w * Math.floor(y); return grid[k]; } function gradient(x, y) { const k = Math.floor(x) + w * Math.floor(y); return [DZX[k], DZY[k]]; } const rng = d3.randomLcg(seed); const cW = w >> 3; const cH = h >> 3; const visited = new Float32Array(cW * cH).fill(-1); let maxstep = 200000; const streaklines = []; const order = d3.shuffle(d3.range(cW * cH)); for (const i of order) { // already visited? drop if (visited[i] > -1) continue; const streakline = []; visited[i] = i; let x; let y; let j = i; let vx, vy, dx, dy; for (const direction of [-1, 1]) { x = (i % cW) * (w / cW); y = Math.floor(i / cW) * (h / cH); vx = 0; vy = 0; do { dx = vx; dy = vy; if (altitude(x, y) <= 0) break; [vx, vy] = gradient(x, y).map((d) => d * direction); const delta = Math.hypot(vx, vy); if (delta === 0) break; x += vx / delta; y += vy / delta; j = Math.floor((x / w) * cW) + cW * Math.floor((y / w) * cW); streakline.push([x, y, i, j]); if (visited[j] < 0) visited[j] = i; if (visited[j] !== i && rng() > 1 - avoidance) break; } while ( vx * dx + vy * dy >= 0 && x > 0 && y > 0 && x < w && y < h && maxstep-- > 0 ); streakline.reverse(); } if (streakline.length > 5) streaklines.push(streakline); } return streaklines; }33:[{"id":"424841d90f27b2816ada638ca903aac1110e1719a7f0e610a20fd6e3f67b0ff7a012bcfed2d14b2793abc89851356527914e1b4bbf30c8f39bcd2701cfde8ed5","href":"https://static.observableusercontent.com/files/424841d90f27b2816ada638ca903aac1110e1719a7f0e610a20fd6e3f67b0ff7a012bcfed2d14b2793abc89851356527914e1b4bbf30c8f39bcd2701cfde8ed5","name":"pantelleria-dem (1).png","type":"image/png","lastModified":1637309243261.249,"size":797494}] 2f:{"title":"Deriving Pantelleria","cells":[{"id":0,"value":"# Deriving Pantelleria","pinned":false,"hidden":false,"mode":"md"},{"id":1138,"value":"*#30DayMapChallenge — Day 21. Elevation.*","pinned":false,"hidden":false,"mode":"md"},{"id":1003,"value":"html`${smaller(\n await FileAttachment(\"pantelleria-dem (1).png\").image()\n)}${smaller(\n heatmap(grid, {\n width: w,\n color: (t) => (t === 0 ? \"lightblue\" : d3.interpolatePiYG(1 - t)),\n dither: false\n })\n)}${smaller(contourMap(grid))}${smaller(\n hillShadingMap(grid, { elevation: 45, azimuth, color: d3.interpolateBrBG })\n)}${smaller(\n heatmap(DZY, {\n width: w,\n color: d3.interpolateRdBu,\n extent: [-18, 18],\n clamp: true\n })\n)}${smaller(streak({ strokeWidth: 0.25, stroke: \"#777\" }))}${smaller(\n heatmap(laplacian, {\n width: w,\n color: (t) => d3.interpolateRdYlBu(1 - t)\n })\n)}${smaller(rain({ DZX, DZY, invalidation }))}${smaller(\n heatmap(\n laplacian.map((d) => Math.min(0, d)),\n {\n width: w,\n extent: [0, -2]\n }\n )\n)}\n\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3,"value":"---\n\n## DEM\nThis notebook starts with a Digital elevation model (DEM) of the Italian island of [Pantelleria](https://it.wikipedia.org/wiki/Pantelleria), halfway between Sicily and the coast of Tunisia in the Mediterranean Sea. I’ve never been there, but it looks like the perfect shape of an island… anyway… let’s begin. I captured this DEM on [nextzen’s terrarium](https://observablehq.com/@fil/reproject-elevation-tiles-detail).","pinned":false,"hidden":false,"mode":"md"},{"id":8,"value":"dem0 = FileAttachment(\"pantelleria-dem (1).png\").image()","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"w = dem0.naturalWidth","pinned":true,"hidden":false,"mode":"ojs"},{"id":48,"value":"h = dem0.naturalHeight","pinned":true,"hidden":false,"mode":"ojs"},{"id":15,"value":"The DEM is ${w} pixels wide by ${h} pixels high. For reference, it corresponds to the following projection (we won’t use it in the notebook).","pinned":false,"hidden":false,"mode":"md"},{"id":23,"value":"projection = d3.geoMercator()\n .scale(297582.21102247556)\n .translate([-61857.49875166787, 206233.6124812287])","pinned":true,"hidden":false,"mode":"ojs"},{"id":28,"value":"Each pixel of the image is an indication of the elevation (altitude) of that point. We can read their red, green and blue values and apply the following formula to obtain an altitude (above the sea level), expressed in meters:","pinned":false,"hidden":false,"mode":"md"},{"id":32,"value":"pixelAltitude = (r, g, b) => 256 * (r - 128) + g + b / 256","pinned":true,"hidden":false,"mode":"ojs"},{"id":1060,"value":"{\n const { r, g, b } = d3.rgb(rgb);\n return html`→ pixelAltitude( ${r}, ${g}, ${b} ) = ${pixelAltitude(\n r,\n g,\n b\n ).toFixed(1)} m.`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1101,"value":"viewof rgb = Inputs.color({ value: d3.rgb(128, 10, 30).formatHex() })","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"If we convert the whole grid, we get an array of a million numbers. (A closer inspection of the DEM reveals that the numbers on the shoreline are varying between 0 and -55, we’ll normalize this to just below 0.)","pinned":false,"hidden":false,"mode":"md"},{"id":39,"value":"grid0 = {\n const w = dem0.naturalWidth;\n const h = dem0.naturalHeight;\n const context = DOM.context2d(w, h, 1);\n context.drawImage(dem0, 0, 0);\n const { data } = context.getImageData(0, 0, w, h);\n return Float32Array.from({ length: w * h }, (_, i) => {\n const pixel = data.subarray(4 * i, 4 * i + 3);\n const a = pixelAltitude(...pixel);\n return a < 0 ? -0.5 : a;\n });\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":131,"value":"The highest point is the Montagna Grande, which in our database appears at an altitude of ${peak.toFixed(2)}m (836m in actuality, I don't know where the difference comes from).","pinned":false,"hidden":false,"mode":"md"},{"id":41,"value":"peak = d3.max(grid0)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1156,"value":"Its color can be computed with the inverse formula:","pinned":false,"hidden":false,"mode":"md"},{"id":1158,"value":"color = (alt) => {\n const r = Math.floor(128 + alt / 256);\n alt -= 256 * (r - 128);\n const g = Math.floor(alt);\n const b = Math.round(256 * (alt - g));\n return [r, g, b];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1163,"value":"html`→ color( ${peak.toFixed(\n 2\n)}m ) = ${color(\n peak\n).join(\", \")} `","pinned":false,"hidden":false,"mode":"ojs"},{"id":57,"value":"## Heightmap\nUsing a simple heatmap of the altitude, we create a recognizable image of the island:","pinned":false,"hidden":false,"mode":"md"},{"id":54,"value":"heatmap(grid, {\n width: w,\n color: (t) => (t === 0 ? \"lightblue\" : d3.interpolatePiYG(1 - t)),\n dither: false\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":53,"value":"import {heatmap} from \"@fil/heatmap\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":1268,"value":"(*Note:* grid *is a slightly blurred copy of *grid0,* making the images smoother. See appendix.*)","pinned":false,"hidden":false,"mode":"md"},{"id":161,"value":"## Line\n\nNow, let’s take a straight line through the island, from west to east, going through that peak, and plot it:","pinned":false,"hidden":false,"mode":"md"},{"id":207,"value":"peakJ = d3.greatestIndex(grid) % w","pinned":true,"hidden":false,"mode":"ojs"},{"id":168,"value":"line = grid.subarray(w * peakJ, w * (peakJ + 1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":166,"value":"Plot.lineY(line).plot({\n height: 180,\n x: { axis: null },\n y: { label: \"elevation (m)\", ticks: 5, grid: true }\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":180,"value":"## Joy Plot\nOk that was a fun use of [Observable Plot](https://observablehq.com/@observablehq/plot)! Let’s push it a bit and make a gratuitous Joy Plot. (This will be my entry for [Day 19. Island](https://observablehq.com/@visionscarto/pantelleria-joy?collection=@visionscarto/30daymapchallenge).)","pinned":false,"hidden":false,"mode":"md"},{"id":205,"value":"joy = Plot.plot({\n marks: d3.range(peakJ - 400, peakJ + 450, 25).map((j) => {\n const line = grid.subarray(w * j, w * (j + 1));\n const dx = j / 8;\n const dy = -110 + j / 4;\n return [\n Plot.areaY(line, { dx, dy, fill: j }),\n Plot.lineY(line, { dx, dy, strokeWidth: 0.75 })\n ];\n }),\n height: 240,\n x: { axis: null },\n y: { axis: null },\n color: { scheme: \"magma\" },\n marginBottom: 150,\n marginRight: 135,\n width\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":374,"value":"## Contours, gradient, hill-shading\nFor an image of this topographic information, we can create [*contours*](https://github.com/d3/d3-contour), or [*hill-shading*](@visionscarto/crater-shading).","pinned":false,"hidden":false,"mode":"md"},{"id":435,"value":"html`\n${small(contourMap(grid))}${small(\n hillShadingMap(grid, { elevation: 45, azimuth, color: d3.interpolateBrBG })\n)}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":457,"value":"viewof azimuth = Inputs.range([0, 360], {\n label: \"azimuth\",\n value: 315,\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":379,"value":"contourMap = grid => heatmap(grid, {\n width: w,\n color: d3.interpolateViridis,\n contours: true\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":479,"value":"For hill-shading, we need to *derive* the altitude matrix into a gradient that encodes the slope. A gradient has two components: DZX, the change in altitude when we move one pixel to the right, and DZY, change of altitude when we move one pixel down.\n\nTwo small technical notes:\n* the units in the gradient are in meters per pixel—we *should* compute everything in a more natural way by dividing by the length of one pixel in meters, but I’m taking this shortcut for now, as it has no incidence on the visuals.\n* since we’re on an island, we can “wrap” the grid—it’s sea level all around the grid.","pinned":false,"hidden":false,"mode":"md"},{"id":475,"value":"DZX = Float32Array.from(grid, (d, i) => Math.max(grid[i + 1] || 0, 0) - d)","pinned":true,"hidden":false,"mode":"ojs"},{"id":474,"value":"DZY = Float32Array.from(grid, (d, i) => Math.max(grid[i + w] || 0, 0) - d)","pinned":true,"hidden":false,"mode":"ojs"},{"id":578,"value":"html`\n${small(\n heatmap(DZX, {\n width: w,\n color: d3.interpolateRdBu,\n extent: [-18, 18],\n clamp: true\n })\n)}${small(\n heatmap(DZY, {\n width: w,\n color: d3.interpolateRdBu,\n extent: [-18, 18],\n clamp: true\n })\n)}`","pinned":true,"hidden":false,"mode":"ojs"},{"id":408,"value":"hillShadingMap = (grid, { elevation = 45, azimuth = 315, color }) => {\n const strength = 0.2;\n\n const r = elevation * (Math.PI / 180);\n const a = -azimuth * (Math.PI / 180);\n const cr = Math.cos(r);\n const sr = Math.sin(r);\n const ca = Math.cos(a);\n const sa = Math.sin(a);\n\n const lz = sr;\n const lx = cr * ca;\n const ly = cr * sa;\n\n const DZX = Float32Array.from(\n grid,\n (d, i) => Math.max(grid[i + 1] || 0, 0) - d\n );\n const DZY = Float32Array.from(\n grid,\n (d, i) => Math.max(grid[i + w] || 0, 0) - d\n );\n\n return heatmap(\n Float32Array.from(grid, (_, i) => {\n const dzx = DZX[i];\n const dzy = DZY[i];\n const norm = Math.hypot(1, dzx, dzy);\n const dot = (lz + dzx * lx + dzy * ly) / norm;\n return 1 + strength * (dot - cr);\n }),\n { width: w, extent: [0.7, 1.15], color }\n );\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":647,"value":"## Rain\n\nNow that we have the gradient, isn’t is very natural to send a bit of rain on the island, and see where the water flows?","pinned":false,"hidden":false,"mode":"md"},{"id":666,"value":"(await visibility(), small(rain({ DZX, DZY, invalidation })))","pinned":true,"hidden":false,"mode":"ojs"},{"id":654,"value":"npoints = w * h","pinned":true,"hidden":false,"mode":"ojs"},{"id":659,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":657,"value":"function fade(context, alpha) {\n const { globalAlpha, globalCompositeOperation } = context;\n context.globalCompositeOperation = \"destination-out\";\n context.fillStyle = \"white\";\n context.globalAlpha = alpha || 1;\n context.fillRect(0, 0, context.canvas.width, context.canvas.height);\n context.globalAlpha = globalAlpha;\n context.globalCompositeOperation = globalCompositeOperation;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":752,"value":"## Streaklines\n\nWe can follow the individual particles of rain floating downhill, and draw their paths. This is called a *streakline* or *streamline* (the two are equivalent on a terrain that doesn’t change with time). We can compute these lines with a simple algorithm loosely based on Whittaker, 1977 ([PDF](https://www.semanticscholar.org/paper/Automated-Streamline-Analysis-Whittaker/768839999d6156f555cf1978a8ca1b4963fdcc86)).","pinned":false,"hidden":false,"mode":"md"},{"id":973,"value":"(await visibility(),\nhtml`

\n${small(streak({ strokeWidth: 0.5, avoidance: 0.5 }))}${small(\n rain({ DZX, DZY, invalidation }),\n \"position: absolute; top: 0;\"\n)}${small(streak({ strokeWidth: 0.25, stroke: \"#777\", avoidance: 0.01 }))}\n`)","pinned":true,"hidden":false,"mode":"ojs"},{"id":849,"value":"function streak({ stroke = \"black\", strokeWidth = 1, ...options } = {}) {\n const context = DOM.context2d(w, h, 1);\n const path = d3.geoPath(d3.geoIdentity()).context(context);\n\n for (const line of streaklines(grid, DZX, DZY, options)) {\n context.moveTo(...line[0]);\n for (const p of line) context.lineTo(...p);\n }\n context.lineWidth = strokeWidth;\n context.strokeStyle = stroke;\n context.stroke();\n\n return context.canvas;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":758,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":470,"value":"## Laplacian\n\nAnother useful derivation of the DEM data is to find peaks and depressions. In this case, the tool to use is the **laplacian**—which has a fancy mathematical definition of second derivatives—or curvature—, that should be explained in the simplest of ways by asking the following question: **is this point higher, or lower, than the mean of its neighbors?** To answer the question, we could use a fancy convolution with a matrix and what not, but I prefer to stick to my simple explanation and compute the difference between the altitude of the point and its value of the *blurred* grid—a way as good as any to measure the mean of neighbors. For this, we’ll use [array-blur](https://github.com/Fil/array-blur).","pinned":false,"hidden":false,"mode":"md"},{"id":495,"value":"B = require(\"array-blur@1\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":499,"value":"laplacian = B.blur()\n .width(w)(grid)\n .map((d, i) => grid[i] - d)\n .map((d) => Math.atan(d / 10)) // dampen \"infinite\" values near the shoreline","pinned":true,"hidden":false,"mode":"ojs"},{"id":620,"value":"Here’s a map of the laplacian, in signed (blue = depression, red = peak), and in absolute value. If you look carefully, you can see that some of the depressions are on the top of the peaks: yes, that is a volcanic island.","pinned":false,"hidden":false,"mode":"md"},{"id":511,"value":"(await visibility(),\nhtml`

\n${small(\n heatmap(laplacian, {\n width: w,\n color: (t) => d3.interpolateRdBu(1 - t)\n })\n)}${small(\n rain({ DZX, DZY, invalidation }),\n \"position: absolute; top: 0;\"\n)}${small(\n heatmap(\n laplacian.map((d) => Math.min(0, d)),\n {\n width: w,\n extent: [0, -2]\n }\n )\n)}\n`)","pinned":true,"hidden":false,"mode":"ojs"},{"id":628,"value":"As a very simple model of hydrography, the laplacian map shows the regions where water will tend to evacuate quickly (in red), or flow into (in blue). The white parts have steady flows (as much water comes in as goes out). In other words, rivers will tend to form in the blue parts of the map. To go beyond that, and measure the volume of flows, we'd have to compute the *accumulation* along the streaklines. (Maybe that will be for an another extension of this notebook.)","pinned":false,"hidden":false,"mode":"md"},{"id":503,"value":"grid = B.blur().radius(0.5).width(w)(grid0) // I blur the grid a little to make the drawings nicer, try setting radius = 5 to make it smoooooth","pinned":true,"hidden":false,"mode":"ojs"},{"id":1250,"value":"---\n*display functions*","pinned":false,"hidden":false,"mode":"md"},{"id":1007,"value":"function small(canvas, style = \"\") {\n return Object.assign(canvas, {\n style: `vertical-align: bottom;max-width: ${\n ((width / 2) | 0) - 1\n }px;${style}`\n });\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1048,"value":"function smaller(canvas, style = \"\") {\n return Object.assign(canvas, {\n style: `vertical-align: bottom;max-width: ${((width / 3) | 0) - 1}px;${style}`\n });\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}