diff --git a/src/data/firmsHeatRaster.js b/src/data/firmsHeatRaster.js new file mode 100644 index 000000000..52a020671 --- /dev/null +++ b/src/data/firmsHeatRaster.js @@ -0,0 +1,176 @@ +/** + * Snap Map–style heat raster for the Active Fires far-zoom bands: detections + * are binned on an equirectangular grid (row 0 = north, column 0 = 180°W), + * blurred, normalized on a log scale and mapped through HEAT_PALETTE. Pure: + * no Cesium, no DOM, so the whole pipeline is unit-tested in node. + */ + +/** Below this normalized heat nothing is painted: empty land stays empty. */ +export const HEAT_FLOOR = 0.03; + +export const HEAT_PALETTE = Object.freeze([ + { at: 0.0, rgb: [43, 108, 255], alpha: 0 }, + { at: 0.08, rgb: [43, 108, 255], alpha: 0.35 }, + { at: 0.25, rgb: [25, 195, 255], alpha: 0.55 }, + { at: 0.45, rgb: [47, 224, 122], alpha: 0.65 }, + { at: 0.65, rgb: [245, 230, 66], alpha: 0.75 }, + { at: 0.82, rgb: [255, 138, 31], alpha: 0.82 }, + { at: 1.0, rgb: [255, 42, 26], alpha: 0.88 }, +]); + +export function heatColorAt(t) { + if (!(t > HEAT_FLOOR)) return [0, 0, 0, 0]; + const v = Math.min(1, t); + let i = 1; + while (i < HEAT_PALETTE.length - 1 && HEAT_PALETTE[i].at < v) i += 1; + const a = HEAT_PALETTE[i - 1]; + const b = HEAT_PALETTE[i]; + const k = (v - a.at) / (b.at - a.at || 1); + const mix = (p, q) => Math.round(p + (q - p) * k); + return [ + mix(a.rgb[0], b.rgb[0]), + mix(a.rgb[1], b.rgb[1]), + mix(a.rgb[2], b.rgb[2]), + Math.round((a.alpha + (b.alpha - a.alpha) * k) * 255), + ]; +} + +const defaultWeight = (fire) => + 1 + Math.log10(1 + Math.max(0, Number(fire.frp) || 0)); + +export function binFires(fires, { width, height, weight = defaultWeight }) { + const grid = new Float32Array(width * height); + for (const fire of fires) { + if (!Number.isFinite(fire.lat) || !Number.isFinite(fire.lon)) continue; + const x = Math.min( + width - 1, + Math.max(0, Math.floor(((fire.lon + 180) / 360) * width)), + ); + const y = Math.min( + height - 1, + Math.max(0, Math.floor(((90 - fire.lat) / 180) * height)), + ); + grid[y * width + x] += weight(fire); + } + return grid; +} + +function kernel(sigma) { + const radius = Math.max(1, Math.ceil(sigma * 3)); + const k = new Float32Array(radius * 2 + 1); + let sum = 0; + for (let i = -radius; i <= radius; i++) { + k[i + radius] = Math.exp(-(i * i) / (2 * sigma * sigma)); + sum += k[i + radius]; + } + for (let i = 0; i < k.length; i++) k[i] /= sum; + return { k, radius }; +} + +export function blurGaussian(grid, width, height, sigma) { + const { k, radius } = kernel(sigma); + const tmp = new Float32Array(grid.length); + const out = new Float32Array(grid.length); + for (let y = 0; y < height; y++) { + const row = y * width; + for (let x = 0; x < width; x++) { + const v = grid[row + x]; + if (v === 0) continue; + for (let i = -radius; i <= radius; i++) { + const xx = (((x + i) % width) + width) % width; // longitude wraps + tmp[row + xx] += v * k[i + radius]; + } + } + } + for (let y = 0; y < height; y++) { + for (let x = 0; x < width; x++) { + const v = tmp[y * width + x]; + if (v === 0) continue; + for (let i = -radius; i <= radius; i++) { + const yy = Math.min(height - 1, Math.max(0, y + i)); // poles clamp + out[yy * width + x] += v * k[i + radius]; + } + } + } + return out; +} + +/** Bins of the log-density histogram used to find the reference quantile. */ +const QUANTILE_BINS = 4096; + +/** + * The p-quantile of the non-zero values, from a histogram of log1p(value): + * one linear pass instead of sorting ~350k burning pixels (a sort cost + * 200 ms+ and dominated the build). Bin width on the log axis is + * max/4096, far finer than the 256-step palette it feeds. + */ +function percentileOfNonZero(values, p) { + let count = 0; + let max = 0; + for (let i = 0; i < values.length; i++) { + const v = values[i]; + if (v > 0) { + count += 1; + if (v > max) max = v; + } + } + if (!count) return 0; + const logMax = Math.log1p(max); + const bins = new Uint32Array(QUANTILE_BINS); + const toBin = (QUANTILE_BINS - 1) / logMax; + for (let i = 0; i < values.length; i++) { + const v = values[i]; + if (v > 0) bins[Math.floor(Math.log1p(v) * toBin)] += 1; + } + const target = Math.floor(p * (count - 1)); + let seen = 0; + for (let b = 0; b < QUANTILE_BINS; b++) { + seen += bins[b]; + if (seen > target) return Math.expm1((b + 1) / toBin); + } + return max; +} + +/** HEAT_PALETTE sampled at 256 steps, so the per-pixel walk is one lookup. */ +const LUT_STEPS = 256; +let paletteLut = null; +function heatLut() { + if (!paletteLut) { + paletteLut = new Uint8ClampedArray(LUT_STEPS * 4); + for (let i = 0; i < LUT_STEPS; i++) { + paletteLut.set(heatColorAt(i / (LUT_STEPS - 1)), i * 4); + } + } + return paletteLut; +} + +export function buildHeatRaster( + fires, + { width, height, sigmaPx, percentile = 0.995 }, +) { + const density = blurGaussian( + binFires(fires, { width, height }), + width, + height, + sigmaPx, + ); + const reference = percentileOfNonZero(density, percentile); + const rgba = new Uint8ClampedArray(width * height * 4); + if (reference > 0) { + const lut = heatLut(); + const scale = (LUT_STEPS - 1) / Math.log1p(reference); + for (let i = 0; i < density.length; i++) { + const d = density[i]; + if (d <= 0) continue; + const step = Math.min(LUT_STEPS - 1, Math.round(Math.log1p(d) * scale)); + const l = step * 4; + if (lut[l + 3] === 0) continue; + const o = i * 4; + rgba[o] = lut[l]; + rgba[o + 1] = lut[l + 1]; + rgba[o + 2] = lut[l + 2]; + rgba[o + 3] = lut[l + 3]; + } + } + return { width, height, rgba, reference }; +} diff --git a/src/data/firmsHeatRaster.test.mjs b/src/data/firmsHeatRaster.test.mjs new file mode 100644 index 000000000..704c83427 --- /dev/null +++ b/src/data/firmsHeatRaster.test.mjs @@ -0,0 +1,73 @@ +import test from 'node:test'; +import assert from 'node:assert/strict'; +import { + HEAT_PALETTE, + heatColorAt, + binFires, + blurGaussian, + buildHeatRaster, +} from './firmsHeatRaster.js'; + +const fire = (lat, lon, frp = 10) => ({ lat, lon, frp }); + +test('palette runs cold to hot and starts fully transparent', () => { + assert.equal(HEAT_PALETTE[0].at, 0); + assert.equal(HEAT_PALETTE.at(-1).at, 1); + assert.deepEqual(heatColorAt(0), [0, 0, 0, 0]); + assert.equal(heatColorAt(0.01)[3], 0, 'below the floor nothing is painted'); + const low = heatColorAt(0.15); + const high = heatColorAt(1); + assert.ok(low[2] > low[0], 'low heat is blue'); + assert.ok(high[0] > high[2], 'the core is red'); + assert.ok(high[3] > low[3], 'hotter is more opaque'); +}); + +test('binFires drops each detection into its equirectangular pixel', () => { + const grid = binFires([fire(89.9, -179.9, 0), fire(-89.9, 179.9, 0)], { + width: 4, + height: 2, + }); + assert.equal(grid[0], 1, 'north-west corner'); + assert.equal(grid[7], 1, 'south-east corner'); + assert.equal( + grid.reduce((a, b) => a + b, 0), + 2, + ); +}); + +test('binFires weights strong fires slightly more', () => { + const grid = binFires([fire(0.1, 0.1, 999)], { width: 4, height: 2 }); + assert.ok(grid[2] > 3 && grid[2] < 5, String(grid[2])); +}); + +test('blurGaussian keeps mass and wraps across the antimeridian', () => { + const w = 16; + const h = 8; + const grid = new Float32Array(w * h); + grid[3 * w + 0] = 1; + const out = blurGaussian(grid, w, h, 1.5); + const sum = out.reduce((a, b) => a + b, 0); + assert.ok(Math.abs(sum - 1) < 1e-3, `mass ${sum}`); + assert.ok(out[3 * w + (w - 1)] > 0, 'spills to the east edge'); +}); + +test('buildHeatRaster is transparent where there are no fires', () => { + const r = buildHeatRaster([fire(40, -120), fire(40.1, -120.1)], { + width: 360, + height: 180, + sigmaPx: 2, + }); + assert.equal(r.rgba.length, 360 * 180 * 4); + const alphaAt = (lat, lon) => { + const x = Math.floor(((lon + 180) / 360) * 360); + const y = Math.floor(((90 - lat) / 180) * 180); + return r.rgba[(y * 360 + x) * 4 + 3]; + }; + assert.ok(alphaAt(40, -120) > 0, 'blob over the fires'); + assert.equal(alphaAt(-30, 60), 0, 'empty far away'); +}); + +test('buildHeatRaster with no fires is fully transparent', () => { + const r = buildHeatRaster([], { width: 8, height: 4, sigmaPx: 1 }); + assert.ok(r.rgba.every((v) => v === 0)); +}); diff --git a/src/layers/firms/heatImagery.js b/src/layers/firms/heatImagery.js new file mode 100644 index 000000000..f2f99fbe6 --- /dev/null +++ b/src/layers/firms/heatImagery.js @@ -0,0 +1,71 @@ +/** + * Drapes the Active Fires heat raster over the globe as one full-extent + * imagery tile on the active imagery host (globe or 3D-tiles), like the wind + * layer's scalar field. Owns the canvas and the imagery layer; nothing else + * in the fires layer touches imagery. + */ + +const defaultCanvas = (width, height) => { + const canvas = document.createElement('canvas'); + canvas.width = width; + canvas.height = height; + return canvas; +}; + +export function createHeatImagery({ + cesium, + getHost, + createCanvas = defaultCanvas, +}) { + let raster = null; + let layer = null; + let collection = null; + let show = true; + + function detach() { + if (layer && collection) collection.remove(layer, true); + layer = null; + collection = null; + } + + function attach() { + detach(); + if (!raster) return; + const host = getHost?.() ?? { collection: null, kind: 'none' }; + if (!host.collection) return; + const canvas = createCanvas(raster.width, raster.height); + const ctx = canvas.getContext('2d'); + const image = ctx.createImageData(raster.width, raster.height); + image.data.set(raster.rgba); + ctx.putImageData(image, 0, 0); + const provider = new cesium.SingleTileImageryProvider({ + url: canvas.toDataURL('image/png'), + tileWidth: raster.width, + tileHeight: raster.height, + rectangle: cesium.Rectangle.MAX_VALUE, + }); + collection = host.collection; + layer = collection.addImageryProvider(provider); + layer.show = show; + collection.raiseToTop?.(layer); + } + + return { + setRaster(next) { + raster = next; + attach(); + }, + setShow(value) { + show = Boolean(value); + if (layer) layer.show = show; + }, + rehome() { + const host = getHost?.(); + if ((host?.collection ?? null) !== collection) attach(); + }, + destroy() { + detach(); + raster = null; + }, + }; +} diff --git a/src/layers/firms/heatImagery.test.mjs b/src/layers/firms/heatImagery.test.mjs new file mode 100644 index 000000000..0fca7b97f --- /dev/null +++ b/src/layers/firms/heatImagery.test.mjs @@ -0,0 +1,87 @@ +import test from 'node:test'; +import assert from 'node:assert/strict'; +import { createHeatImagery } from './heatImagery.js'; + +function harness() { + const layers = []; + const removed = []; + const collection = { + addImageryProvider(provider) { + const layer = { provider, show: true, alpha: 1 }; + layers.push(layer); + return layer; + }, + remove(layer, destroy) { + layers.splice(layers.indexOf(layer), 1); + removed.push({ layer, destroy }); + }, + raiseToTop() {}, + }; + let host = { collection, kind: 'globe' }; + const cesium = { + Rectangle: { MAX_VALUE: { max: true } }, + SingleTileImageryProvider: class { + constructor(options) { + this.options = options; + } + }, + }; + const createCanvas = (w, h) => ({ + width: w, + height: h, + getContext: () => ({ + createImageData: (cw, ch) => ({ + data: new Uint8ClampedArray(cw * ch * 4), + }), + putImageData() {}, + }), + toDataURL: () => 'data:image/png;base64,AAAA', + }); + const imagery = createHeatImagery({ + cesium, + getHost: () => host, + createCanvas, + }); + return { imagery, layers, removed, collection, setHost: (h) => (host = h) }; +} + +const raster = { width: 4, height: 2, rgba: new Uint8ClampedArray(32) }; + +test('setRaster drapes one full-globe tile and replaces it on the next build', () => { + const h = harness(); + h.imagery.setRaster(raster); + assert.equal(h.layers.length, 1); + assert.deepEqual(h.layers[0].provider.options.rectangle, { max: true }); + assert.equal(h.layers[0].provider.options.tileWidth, 4); + h.imagery.setRaster(raster); + assert.equal(h.layers.length, 1, 'never two heat layers'); + assert.equal(h.removed[0].destroy, true); +}); + +test('setShow toggles without rebuilding', () => { + const h = harness(); + h.imagery.setRaster(raster); + const layer = h.layers[0]; + h.imagery.setShow(false); + assert.equal(layer.show, false); + h.imagery.setShow(true); + assert.equal(h.layers[0], layer); + assert.equal(layer.show, true); +}); + +test('no imagery host keeps nothing draped; rehome brings it back', () => { + const h = harness(); + h.setHost({ collection: null, kind: 'none' }); + h.imagery.setRaster(raster); + assert.equal(h.layers.length, 0); + h.setHost({ collection: h.collection, kind: 'globe' }); + h.imagery.rehome(); + assert.equal(h.layers.length, 1); +}); + +test('destroy removes the layer', () => { + const h = harness(); + h.imagery.setRaster(raster); + h.imagery.destroy(); + assert.equal(h.layers.length, 0); +});