31:T515,## *Overall preliminary outcome of this exploration* [JPEG-XL](https://jpeg.org/jpegxl/) ([Wikipedia](https://en.wikipedia.org/wiki/JPEG_XL)), especially in lossless modular mode, seems to be the most straightforward way short of coming up with a novel compression scheme. The standard supports: - channels with up to 32 bit precision, 16bit is maximum with png (16 bits are nice, nbut not enough for sub-meter precision) - having the data in a single channel make the data easier to compress, since there are less hard jumps, than with splitting up the channel into multiple ones. - significantly better lossless compression performance than png - [GDAL](https://gdal.org/drivers/raster/jpegxl.html) - native decoding in [some browsers](https://caniuse.com/jpegxl), or through polyfills like [jxl.js](https://github.com/niutech/jxl.js) that uses WASM web workers Further down the pipeline in the browser, the data needs to be still mapped to the 3 or 4 8-bit channels of a standard png or ImageData to be used in texture for WebGl where 8bit RGBA is still the maximum as of mid 2023. But since the goal is to be able to dynamically map this high bit depth value at a local value range to the limited bit depth of an output color ramp in a shader, this is not a limitiation for this use case.32:T4d0,# Motivation This Notebook explores the feasibility of storing elevation data (rearrain raster, DEM, DTM, or effectively any high channel resolution / bit depth raster channel) at maximum possible bit depth in widely supported modern web image codecs like webp and avif. Maximum target channel bit depth is 24 bit in this scenario. Goal is to find possible file size improvements over delivery as PNG at (near) lossless quality. Encoding in single channel is not suffienct for the elevation data use case since those codecs allow at maximum 12 bits per channel. 16-bit grayscale PNG is possible, but has poor compression. The idea is to have a computationally simple mapping of single channel data to multiple channels of a given image format. The compression step can than be handled by standard image codecs that might be accelarated by existing hardware implementations. One abandoned Idea was to plot directly into JPEGs Y`CbCr color space and then transform that to RGB, but to the need to handle the out of gamut values, this approach seemed to complicate thing more than needed and ad one more transformation for all values. Distributing the Value to RGB channels directly at differen bit depth seems to work quite well.33:[] 2f:{"title":"Storing lossless earth elevation data in modern image formats","cells":[{"id":0,"value":"# Storing lossless earth elevation data in modern image formats","pinned":false,"hidden":false,"mode":"md"},{"id":605,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":1672,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":600,"value":"## Basic Parameters","pinned":false,"hidden":false,"mode":"md"},{"id":1624,"value":"Optimal discovered data pre-processing settings for optimal compression / minimal file size at 20 bit resolution with 100% JPEG quality:\n- Conversion Mode: Experimental (ternary bitwise interleaving of Gray code mapped input value)\n- Stretch: None\n- Channel Mapping: no significant influence","pinned":false,"hidden":false,"mode":"md"},{"id":1631,"value":"viewof applyOptimalButton = Inputs.button(\n \"Apply JPEG Size Optimized Settings\",\n { reduce: applyOptimal }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":171,"value":"viewof elevationMin = Inputs.range([-20000, 20000], {\n value: -485.76,\n label: \"minimum\",\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":195,"value":"viewof elevationMax = Inputs.range([-100000, 100000], {\n label: \"maximum\",\n value: 10000,\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":227,"value":"At 20 bits precision it's possible to get 1cm precision for the dry surface elevation range of earth\nthe minimum value of -485.76 is chosen to get to an increment value of 1cm exactly","pinned":false,"hidden":false,"mode":"md"},{"id":203,"value":"params = ({\n bits,\n steps,\n increment: elevationInterval / steps,\n factor: steps / elevationInterval,\n elevation: {\n interval: elevationInterval,\n min: elevationMin,\n max: elevationMax\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":842,"value":"viewof conversionMode = Inputs.radio(Object.keys(CONVERSION_MODES), {\n label: \"Conversion Mode\",\n value: \"CUSTOM\",\n disabled: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":550,"value":"viewof selectedStrategy = Inputs.radio(Object.keys(STRATEGIES), {\n label: \"Custom Modes\",\n value: \"INTERLEAVE\",\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":151,"value":"viewof bits = Inputs.range([1, 24], {\n value: 20,\n step: 1,\n label: \"Bits\",\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":859,"value":"Stretch Channel Values to full channel bit depth or only RGB/YCBR safe","pinned":false,"hidden":false,"mode":"md"},{"id":835,"value":"viewof channelBitOffsetMode = Inputs.radio(Object.keys(OFFSET_MODES), {\n label: \"Stretch\",\n value: \"FULL\",\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":838,"value":"viewof channelOrderSelection = Inputs.select(channelPermutations, {\n label: \"Channel Mapping\",\n value: \"\",\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"refCanvas = canvasFromImageData(imageData)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1605,"value":"refPNGSize = refCanvas.toDataURL(\"image/png\").length","pinned":false,"hidden":false,"mode":"ojs"},{"id":1612,"value":"jpgSize = encodedImage.data.length","pinned":false,"hidden":false,"mode":"ojs"},{"id":1615,"value":"ratio = ((jpgSize / refPNGSize) * 100).toPrecision(5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1600,"value":"viewof jpgQuality = Inputs.range([1, 100], {\n value: 50,\n step: 1,\n label: \"JPEG Quality\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1589,"value":"encodedImage = jpeg.encode(imageData, jpgQuality)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1595,"value":"jpgToImg(encodedImage)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1572,"value":"jpgToImg = (encodedImage) => {\n const imageBlob = new Blob([encodedImage.data], {\n type: \"image/jpeg\"\n });\n const imageUrl = URL.createObjectURL(imageBlob);\n const imgElement = document.createElement(\"img\");\n imgElement.src = imageUrl;\n imgElement.style.width = encodedImage.width / window.devicePixelRatio + \"px\";\n return imgElement;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":709,"value":"canvasChannelsFromImageData(imageData)","pinned":false,"hidden":false,"mode":"ojs"},{"id":834,"value":"viewof elevationOffset = Inputs.range([-2000, 10000], {\n label: \"offset Input Elevation\",\n value: 0,\n step: 10\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":833,"value":"viewof elevationExagerration = Inputs.range([-3, 3], {\n label: \"Input Elevation exagerration power\",\n value: 0,\n step: 0.1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1309,"value":"## Elevation Data View","pinned":false,"hidden":false,"mode":"md"},{"id":469,"value":"viewof dataSource = Inputs.radio([\"pattern\", \"arc13\"], {\n label: \"Sample Data\",\n value: \"arc13\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1322,"value":"viewof stretchView = Inputs.toggle({\n value: true,\n label: \"'auto-expose' elevation range\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1303,"value":"viewof bitsScale = Inputs.range([1, 12], {\n value: 10,\n step: 1,\n label: \"Color Scale Bits\",\n disabled: notCustom\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1223,"value":"viewof selectedColorScale = Inputs.select(Object.keys(colorScaleBuffers), {\n label: \"Color Scale\",\n value: \"Viridis\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1231,"value":"canvasFromImageData(\n integerBufferToColorizedImageData(\n elevationBuffer,\n colorScaleBuffers[selectedColorScale],\n { stretch: stretchView }\n )\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":683,"value":"## Functions","pinned":false,"hidden":false,"mode":"md"},{"id":1055,"value":"CONVERSION_MODES = createEnum([\"CUSTOM\", \"MB\", \"TERRARIUM\"])","pinned":false,"hidden":false,"mode":"ojs"},{"id":534,"value":"STRATEGIES = createEnum([\"SPLIT8\", \"SPLIT_EQUAL\", \"INTERLEAVE\", \"EXPERIMENTAL\"])","pinned":false,"hidden":false,"mode":"ojs"},{"id":969,"value":"OFFSET_MODES = createEnum([\"NONE\", \"FULL\", \"YCbCrSAFE\"])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1476,"value":"### Elevation/Scalar to Raster","pinned":false,"hidden":false,"mode":"md"},{"id":348,"value":"elevationBuffer = elevationRasterToIntegerBuffer(elevationData)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1168,"value":"imageData = integerBufferToImageData(elevationBuffer, {\n customChannelOrder,\n strategy: STRATEGIES[selectedStrategy]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1146,"value":"function elevationRasterToIntegerBuffer(raster) {\n let { width, height, e } = raster;\n let eUint32 = {\n buffer: new Uint32Array(e.length),\n width,\n height,\n bitDepth: params.bits\n };\n\n let min = 1 << bits;\n let max = 0;\n e.forEach((v, i) => {\n v = (v + elevationOffset) * Math.pow(10, elevationExagerration);\n let valueInt = mapToIntegers(v);\n min = Math.min(valueInt, min);\n max = Math.max(valueInt, max);\n eUint32.buffer[i] = valueInt;\n });\n if (eUint32.buffer.length > 0) {\n eUint32.max = max;\n eUint32.min = min;\n }\n return eUint32;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":419,"value":"function integerBufferToImageData(\n eUint32,\n { strategy, customChannelOrder } = {}\n) {\n let { width, height, buffer } = eUint32;\n let rgb = new ImageData(width, height);\n // make opaque\n rgb.data.fill(255);\n buffer.forEach((v, i) => {\n let vRGB = integerToRGB(v, { strategy, customChannelOrder });\n const iUint32 = i * 4;\n rgb.data.set(vRGB, iUint32);\n });\n return rgb;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":38,"value":"function integerToRGB(\n valueInt,\n { strategy = STRATEGIES.INTERLEAVE, customChannelOrder = [0, 1, 2] } = {}\n) {\n const [R, G, B] = customChannelOrder;\n let triple = undefined;\n if (strategy === STRATEGIES.SPLIT8) {\n triple = split24bitsTo3at8bits(valueInt);\n } else if (strategy === STRATEGIES.SPLIT_EQUAL) {\n triple = splitEqual(valueInt);\n } else if (strategy === STRATEGIES.INTERLEAVE) {\n triple = bitwiseInterleave(valueInt);\n } else if (strategy === STRATEGIES.EXPERIMENTAL) {\n triple = experimental(valueInt);\n }\n\n if (OFFSET_MODES[channelBitOffsetMode] !== OFFSET_MODES.NONE) {\n return new Uint8ClampedArray([\n triple[R] << channelBitOffset[R],\n triple[G] << channelBitOffset[G],\n triple[B] << channelBitOffset[B]\n ]);\n } else {\n return new Uint8ClampedArray([triple[R], triple[G], triple[B]]);\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1470,"value":"### Canvas","pinned":false,"hidden":false,"mode":"md"},{"id":1156,"value":"### Integer Value Splitting methods","pinned":false,"hidden":false,"mode":"md"},{"id":827,"value":"function split24bitsTo3at8bits(valueInt) {\n return new Uint8ClampedArray([\n (valueInt >> 16) & mask8bit,\n (valueInt >> 8) & mask8bit,\n valueInt & mask8bit\n ]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":997,"value":"function splitEqual(valueInt) {\n return new Uint8ClampedArray([\n (valueInt >> outBitsOffsets[0]) & outBitMasks[0],\n (valueInt >> outBitsOffsets[1]) & outBitMasks[1],\n valueInt & outBitMasks[2]\n ]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":611,"value":"function bitwiseInterleave(value) {\n let outChannels = new Uint8ClampedArray(channels).fill(0);\n outChannels.forEach((v, i, a) => {\n const channelBits = outBits[i];\n for (let channelBitPos = 0; channelBitPos < channelBits; channelBitPos++) {\n let inputBitPos = channelBitPos * channels + i;\n const bit = (value >> inputBitPos) & 1;\n outChannels[i] |= bit << channelBitPos;\n }\n });\n return outChannels;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":819,"value":"function experimental(value) {\n // function to evaluate nother methods\n const grayValue = value ^ (value >> 1);\n let outChannels = new Array(channels).fill(0);\n\n outChannels.forEach((v, i, a) => {\n const channelBits = outBits[i];\n for (let channelBitPos = 0; channelBitPos < channelBits; channelBitPos++) {\n let inputBitPos = channelBitPos * channels + i;\n const bit = (grayValue >> inputBitPos) & 1;\n outChannels[i] |= bit << channelBitPos;\n }\n });\n return outChannels;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":618,"value":"function distributeBitsPerChannel(bits, channels = 3) {\n return [...Array(channels)].map(\n (v, i) => Math.floor(bits / channels) + (i < bits % channels ? 1 : 0)\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":276,"value":"function mapToIntegers(value) {\n let clampedValue = Math.max(Math.min(elevationMax, value), elevationMin);\n let int = parseInt((clampedValue - elevationMin) * params.factor);\n return int;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":680,"value":"### Elevation Data Providers","pinned":false,"hidden":false,"mode":"md"},{"id":7,"value":"geotiff = require(\"geotiff@2.0.7/dist-browser/geotiff.js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":476,"value":"arc13 = (\n await geotiff.fromUrl(\n \"https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n34w117/USGS_13_n34w117.tif\"\n )\n).getImage(4)","pinned":false,"hidden":false,"mode":"ojs"},{"id":31,"value":"readElevationData = async (image, options = {}) => {\n const rasters = await image.readRasters(options);\n let [e] = rasters;\n let { width, height } = rasters;\n return { e, width, height };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":493,"value":"function generateElevationSample(width = 676, height = 676) {\n const range = elevationMax - elevationMin;\n const increment = range / width;\n const e = new Float32Array(width * height);\n for (let x = 0; x < width; x++) {\n // let rowFactor = Math.sin(x / 60) * 5000 + 4800;\n for (let y = 0; y < height; y++) {\n e[y * width + x] =\n noise2DOctaves(x / width, y / height, 8, 0.5) * 2000 + 1800;\n }\n }\n return { e, width, height };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":147,"value":"elevationData = dataSource === \"arc13\"\n ? readElevationData(arc13)\n : generateElevationSample()","pinned":false,"hidden":false,"mode":"ojs"},{"id":805,"value":"### Reference RGB to Meter conversions","pinned":false,"hidden":false,"mode":"md"},{"id":804,"value":"mbdem = ({\n meterToRGB: (meter) => {\n const value = (meter + 10000) / 0.1;\n const R = value >>> 16;\n const G = (value >>> 8) & 0xff;\n const B = value & 0xff;\n return new Uint8ClampedArray([R, G, B]);\n },\n RGBToMeter: (R, G, B) => {\n return -10000 + ((R << 16) | (G << 8) | B) * 0.1;\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":803,"value":"mapzenTerrarium = ({\n meterToRGB: (meter) => {\n const value = meter + 32768;\n const R = value >>> 8;\n const G = value & 0xff;\n const B = Math.floor(((value << 8) & 0xff) * 256);\n return new Uint8ClampedArray([R, G, B]);\n },\n RGBToMeter: (R, G, B) => {\n return ((R << 8) | G) + B / 256 - 32768;\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":670,"value":"### Color Space Encoding","pinned":false,"hidden":false,"mode":"md"},{"id":922,"value":"channelNames = \"Rd,Gr,Bl\".split(\",\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":802,"value":"channels = channelNames.length","pinned":false,"hidden":false,"mode":"ojs"},{"id":1354,"value":"indexedChannels = channelNames.reduce(\n (a, v, i) => Object.assign({ [v]: i }, a),\n {}\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1422,"value":"customChannelOrder = channelNames.map((v) => channelOrderSelection.indexOf(v))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1420,"value":"channelOrder = channelOrderSelection","pinned":false,"hidden":false,"mode":"ojs"},{"id":938,"value":"channelPermutations = getPermutationsOfArrayItems(channelNames)","pinned":false,"hidden":false,"mode":"ojs"},{"id":959,"value":"channelBitDepth = 8","pinned":false,"hidden":false,"mode":"ojs"},{"id":961,"value":"channelBitDepths = Array(channels).fill(channelBitDepth)","pinned":false,"hidden":false,"mode":"ojs"},{"id":801,"value":"outBits = distributeBitsPerChannel(bits, 3)","pinned":false,"hidden":false,"mode":"ojs"},{"id":825,"value":"leastSignificantInterleavedChannelIndex = getLeastSignificantIndex(bits)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1004,"value":"outBitMasks = outBits.map(getBitmask)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1007,"value":"outBitsOffsets = outBits.map((v, i, a) => sum(a.slice(i)) - v)","pinned":false,"hidden":false,"mode":"ojs"},{"id":979,"value":"channelBitOffsets = ({\n FULL: outBits.map(getFullChannelStretchOffset),\n YCbCrSAFE: channelBitDepths.map((v, i) =>\n outBits[i] <= safeRGBchannelDepthsForYCbCr[i]\n ? v - safeRGBchannelDepthsForYCbCr[i]\n : 0\n )\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":956,"value":"channelBitOffset = channelBitOffsets[channelBitOffsetMode]","pinned":false,"hidden":false,"mode":"ojs"},{"id":799,"value":"function getFullChannelStretchOffset(v, i) {\n return channelBitDepths[i] - v;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":942,"value":"#### TBD","pinned":false,"hidden":false,"mode":"md"},{"id":905,"value":"safeRGBchannelDepthsForYCbCr = [7, 7, 7] // TBD","pinned":false,"hidden":false,"mode":"ojs"},{"id":662,"value":"### Helpers","pinned":false,"hidden":false,"mode":"md"},{"id":1207,"value":"#### Colors","pinned":false,"hidden":false,"mode":"md"},{"id":1521,"value":"colorScaleSteps = 1024","pinned":false,"hidden":false,"mode":"ojs"},{"id":1519,"value":"colorScaleBuffers = createColorScaleLUTBuffers(colorScaleSteps)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1082,"value":"#### Contextual","pinned":false,"hidden":false,"mode":"md"},{"id":796,"value":"notCustom = CONVERSION_MODES[conversionMode] !== CONVERSION_MODES.CUSTOM","pinned":false,"hidden":false,"mode":"ojs"},{"id":795,"value":"elevationInterval = elevationMax - elevationMin","pinned":false,"hidden":false,"mode":"ojs"},{"id":794,"value":"steps = 1 << bits","pinned":false,"hidden":false,"mode":"ojs"},{"id":1085,"value":"#### Generic","pinned":false,"hidden":false,"mode":"md"},{"id":1009,"value":"sum = (a) =>\n a.reduce(function (a, b) {\n return a + b;\n })","pinned":false,"hidden":false,"mode":"ojs"},{"id":538,"value":"mask8bit = getBitmask(8)","pinned":false,"hidden":false,"mode":"ojs"},{"id":999,"value":"getBitmask = (l) => (1 << l) - 1","pinned":false,"hidden":false,"mode":"ojs"},{"id":826,"value":"function getLeastSignificantIndex(bits, channels = 3) {\n return ((bits % channels) + channels - 1) % channels;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1641,"value":"#### Inputs","pinned":false,"hidden":false,"mode":"md"},{"id":1643,"value":"function applyOptimal() {\n console.log(\"A\");\n inputValueSetter(viewof selectedStrategy, \"EXPERIMENTAL\");\n inputValueSetter(viewof channelBitOffsetMode, \"NONE\");\n inputValueSetter(viewof jpgQuality, 100);\n return null;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1066,"value":"### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":345,"value":"import {\n canvasFromImage,\n canvasFromImageData,\n canvasChannelsFromImageData,\n integerBufferToColorizedImageData\n} from \"@mxfh/2d-canvas-utils\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1511,"value":"import {\n createColorScaleLUTBuffers,\n createEnum,\n getPermutationsOfArrayItems,\n inputValueSetter\n} from \"@mxfh/utils\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":791,"value":"import { noise2DOctaves } from \"@bensimonds/perlin-noise\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1118,"value":"jpeg = import(\"https://cdn.skypack.dev/jpeg-js@0.4.4?min\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":585,"value":"### Todo\n- validate usable range of YCbCr Color Space (safeRGBchannelDepthsForYCbCr)\n- Documentation / Write Up\n- benchmark with web image codecs\n- implement baseline comparisons vs LERC, Terraria PNG\n- get error metrics\n- reusable codec methods\n \n #### further options\n- split data in freqency domain (FFT)\n- explore chroma subsampling strategies\n- restructure loops for speed\n- webgl implementation\n- dynamic function to move least significant bit channel into chroma for any target bit resolition\n- encode ocean bathymetry at 10cm resolution in the value range 9000m to 10000m\n- bit flags for designating waterbodies\n- use of remainder bits or alpha channnel to encode landcover information to enable PBR shading","pinned":false,"hidden":false,"mode":"md"}],"files":"$W33","stdlib":"1","titling":"auto"}