From 15a39b6e42d3f1202ce8dae81dfe96768c92858d Mon Sep 17 00:00:00 2001 From: Alex Kovac <5686693+kovaca@users.noreply.github.com> Date: Sun, 27 Oct 2024 00:25:25 -0600 Subject: [PATCH] add wip gdaldem histojams --- package-lock.json | 13 +- package.json | 3 +- src/data/dicopaldata.json.js | 20 + src/data/histexample.json | 780 +++++++++++++++++++++++++++++++++++ src/histojam.md | 225 ++++++++++ 5 files changed, 1039 insertions(+), 2 deletions(-) create mode 100644 src/data/dicopaldata.json.js create mode 100644 src/data/histexample.json create mode 100644 src/histojam.md diff --git a/package-lock.json b/package-lock.json index 940fca5..45d6fbd 100644 --- a/package-lock.json +++ b/package-lock.json @@ -8,7 +8,8 @@ "@observablehq/framework": "^1.10.1", "d3": "^7.9.0", "d3-dsv": "^3.0.1", - "d3-time-format": "^4.1.0" + "d3-time-format": "^4.1.0", + "dicopal": "^0.8.1" }, "devDependencies": { "rimraf": "^5.0.5" @@ -1652,6 +1653,16 @@ "npm": "1.2.8000 || >= 1.4.16" } }, + "node_modules/dicopal": { + "version": "0.8.1", + "resolved": "https://registry.npmjs.org/dicopal/-/dicopal-0.8.1.tgz", + "integrity": "sha512-wt5O5plia2sCmfOQuu8b8yxkXsGSEvGPeLNrZSC4JDVtxcWv2BD/TSV+Uk1xuUm+p86kGokHIWXIXNMtOvB/4Q==", + "dependencies": { + "d3-array": "^3.2.4", + "d3-interpolate": "^3.0.1", + "d3-scale": "^4.0.2" + } + }, "node_modules/eastasianwidth": { "version": "0.2.0", "resolved": "https://registry.npmjs.org/eastasianwidth/-/eastasianwidth-0.2.0.tgz", diff --git a/package.json b/package.json index a7feeb1..862d981 100644 --- a/package.json +++ b/package.json @@ -12,7 +12,8 @@ "@observablehq/framework": "^1.10.1", "d3": "^7.9.0", "d3-dsv": "^3.0.1", - "d3-time-format": "^4.1.0" + "d3-time-format": "^4.1.0", + "dicopal": "^0.8.1" }, "devDependencies": { "rimraf": "^5.0.5" diff --git a/src/data/dicopaldata.json.js b/src/data/dicopaldata.json.js new file mode 100644 index 0000000..608ec3b --- /dev/null +++ b/src/data/dicopaldata.json.js @@ -0,0 +1,20 @@ +import * as dicopal from "dicopal"; + +const div = dicopal.getPalettes({ type: "diverging", number: 7 }); +const seq = dicopal.getPalettes({ type: "sequential", number: 7 }); +const custom = [ + {"type":"sequential","colors":["#E4EEC6","#C3E0A7","#8DCB81","#64AD66","#378C4D","#296634","#1C401F"]}, + {"type":"sequential","colors":["#D9E9F5","#BACCE8","#8D9CCD","#6678B8","#3C59A6","#214080","#001C5E"]}, + {"type":"sequential","colors":["#E1F0DC","#BCE1D0","#94D2C4","#5EB4B4","#0098A6","#03707D","#054957"]}, + {"type":"sequential","colors":["#FEE1D6","#FEC0BF","#F496A0","#E16E87","#C64974","#A32A64","#77184F"]}, + {"type":"sequential","colors":["#FAECB7","#FDD881","#FCB817","#F49D21","#EA8026","#BA4F28","#8B1B26"]}, + {"type":"diverging","colors":["#8B1B26","#BA4F28","#F08E23","#FCDD8E","#FBF7EA","#B6DED2","#6CAFAC","#07717D","#054957"]}, + {"type":"diverging","colors":["#5A0C0C","#A03035","#E1746F","#F9C5BC","#FBF7EA","#CFDBE8","#90A2D1","#3A56A3","#001C5E"]}, + {"type":"diverging","colors":["#5C1339","#A33166","#D783A0","#FCD5D2","#FBF7EA","#E5E8B3","#84C77B","#2F8749","#1C401F"]} +] + +const cmaps = [...custom, ...seq, ...div] + + + +process.stdout.write(JSON.stringify(cmaps)); \ No newline at end of file diff --git a/src/data/histexample.json b/src/data/histexample.json new file mode 100644 index 0000000..ce28166 --- /dev/null +++ b/src/data/histexample.json @@ -0,0 +1,780 @@ +{ + "description":"example.tif", + "driverShortName":"GTiff", + "driverLongName":"GeoTIFF", + "files":[ + "example.tif", + "example.aux.xml" + ], + "size":[ + 12000, + 6000 + ], + "coordinateSystem":{ + "wkt":"GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]", + "dataAxisToSRSAxisMapping":[ + 2, + 1 + ] + }, + "geoTransform":[ + -180.0, + 0.03, + 0.0, + 90.0, + 0.0, + -0.03 + ], + "metadata":{ + "":{ + "AREA_OR_POINT":"Area" + }, + "IMAGE_STRUCTURE":{ + "COMPRESSION":"LZW", + "INTERLEAVE":"BAND" + } + }, + "cornerCoordinates":{ + "upperLeft":[ + -180.0, + 90.0 + ], + "lowerLeft":[ + -180.0, + -90.0 + ], + "lowerRight":[ + 180.0, + -90.0 + ], + "upperRight":[ + 180.0, + 90.0 + ], + "center":[ + -0.0, + 0.0 + ] + }, + "wgs84Extent":{ + "type":"Polygon", + "coordinates":[ + [ + [ + -180.0, + 90.0 + ], + [ + -180.0, + -90.0 + ], + [ + 180.0, + -90.0 + ], + [ + 180.0, + 90.0 + ], + [ + -180.0, + 90.0 + ] + ] + ] + }, + "bands":[ + { + "band":1, + "block":[ + 12000, + 1 + ], + "type":"Float32", + "colorInterpretation":"Gray", + "min":0.005, + "max":23.087, + "minimum":0.005, + "maximum":23.087, + "mean":2.277, + "stdDev":1.433, + "histogram":{ + "count":256, + "min":-0.040118246095474978, + "max":23.132044400310761, + "buckets":[ + 402, + 2045, + 3115, + 3826, + 4943, + 7565, + 12445, + 19043, + 25057, + 30610, + 34246, + 36798, + 39375, + 40312, + 37219, + 32719, + 28820, + 27129, + 26130, + 25316, + 24571, + 24647, + 23949, + 23648, + 22949, + 21096, + 19483, + 18079, + 16486, + 15395, + 14843, + 14102, + 13736, + 13208, + 12981, + 12842, + 12869, + 12521, + 12426, + 12088, + 12015, + 11760, + 11410, + 10841, + 10485, + 9868, + 9425, + 8690, + 8069, + 7441, + 6843, + 6245, + 5798, + 5109, + 4664, + 4168, + 3715, + 3486, + 3139, + 2773, + 2639, + 2342, + 1999, + 1888, + 1695, + 1591, + 1429, + 1317, + 1171, + 1018, + 926, + 882, + 733, + 743, + 613, + 569, + 487, + 473, + 399, + 437, + 373, + 344, + 325, + 276, + 264, + 249, + 235, + 198, + 213, + 210, + 193, + 158, + 150, + 151, + 134, + 114, + 116, + 100, + 109, + 109, + 108, + 91, + 85, + 69, + 84, + 63, + 70, + 54, + 53, + 54, + 59, + 44, + 41, + 29, + 35, + 37, + 35, + 44, + 34, + 26, + 32, + 37, + 26, + 21, + 16, + 24, + 24, + 16, + 26, + 12, + 17, + 17, + 20, + 14, + 16, + 17, + 14, + 12, + 9, + 10, + 7, + 11, + 9, + 7, + 10, + 9, + 12, + 6, + 10, + 7, + 4, + 6, + 4, + 6, + 5, + 9, + 4, + 9, + 8, + 5, + 4, + 2, + 4, + 1, + 2, + 2, + 3, + 4, + 3, + 3, + 4, + 2, + 2, + 4, + 2, + 2, + 3, + 4, + 1, + 0, + 1, + 1, + 3, + 0, + 1, + 2, + 1, + 1, + 2, + 1, + 1, + 2, + 1, + 0, + 1, + 0, + 0, + 1, + 3, + 0, + 0, + 0, + 1, + 0, + 2, + 1, + 1, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 1, + 1, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1 + ] + }, + "noDataValue":-3.4e+38, + "metadata":{ + "":{ + "STATISTICS_COVARIANCES":"2.054343908110655", + "STATISTICS_MAXIMUM":"23.086786270142", + "STATISTICS_MEAN":"2.2774180973565", + "STATISTICS_MINIMUM":"0.0051398840732872", + "STATISTICS_SKIPFACTORX":"1", + "STATISTICS_SKIPFACTORY":"1", + "STATISTICS_STDDEV":"1.4332982620902" + } + } + } + ], + "stac":{ + "proj:shape":[ + 6000, + 12000 + ], + "proj:wkt2":"GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]", + "proj:epsg":4326, + "proj:projjson":{ + "$schema":"https://proj.org/schemas/v0.7/projjson.schema.json", + "type":"GeographicCRS", + "name":"WGS 84", + "datum_ensemble":{ + "name":"World Geodetic System 1984 ensemble", + "members":[ + { + "name":"World Geodetic System 1984 (Transit)", + "id":{ + "authority":"EPSG", + "code":1166 + } + }, + { + "name":"World Geodetic System 1984 (G730)", + "id":{ + "authority":"EPSG", + "code":1152 + } + }, + { + "name":"World Geodetic System 1984 (G873)", + "id":{ + "authority":"EPSG", + "code":1153 + } + }, + { + "name":"World Geodetic System 1984 (G1150)", + "id":{ + "authority":"EPSG", + "code":1154 + } + }, + { + "name":"World Geodetic System 1984 (G1674)", + "id":{ + "authority":"EPSG", + "code":1155 + } + }, + { + "name":"World Geodetic System 1984 (G1762)", + "id":{ + "authority":"EPSG", + "code":1156 + } + }, + { + "name":"World Geodetic System 1984 (G2139)", + "id":{ + "authority":"EPSG", + "code":1309 + } + }, + { + "name":"World Geodetic System 1984 (G2296)", + "id":{ + "authority":"EPSG", + "code":1383 + } + } + ], + "ellipsoid":{ + "name":"WGS 84", + "semi_major_axis":6378137, + "inverse_flattening":298.257223563 + }, + "accuracy":"2.0", + "id":{ + "authority":"EPSG", + "code":6326 + } + }, + "coordinate_system":{ + "subtype":"ellipsoidal", + "axis":[ + { + "name":"Geodetic latitude", + "abbreviation":"Lat", + "direction":"north", + "unit":"degree" + }, + { + "name":"Geodetic longitude", + "abbreviation":"Lon", + "direction":"east", + "unit":"degree" + } + ] + }, + "scope":"Horizontal component of 3D system.", + "area":"World.", + "bbox":{ + "south_latitude":-90, + "west_longitude":-180, + "north_latitude":90, + "east_longitude":180 + }, + "id":{ + "authority":"EPSG", + "code":4326 + } + }, + "proj:transform":[ + -180.0, + 0.03, + 0.0, + 90.0, + 0.0, + -0.03 + ], + "raster:bands":[ + { + "data_type":"float32", + "stats":{ + "minimum":0.005, + "maximum":23.087, + "mean":2.277, + "stddev":1.433 + }, + "histogram":{ + "count":256, + "min":-0.040118246095474978, + "max":23.132044400310761, + "buckets":[ + 402, + 2045, + 3115, + 3826, + 4943, + 7565, + 12445, + 19043, + 25057, + 30610, + 34246, + 36798, + 39375, + 40312, + 37219, + 32719, + 28820, + 27129, + 26130, + 25316, + 24571, + 24647, + 23949, + 23648, + 22949, + 21096, + 19483, + 18079, + 16486, + 15395, + 14843, + 14102, + 13736, + 13208, + 12981, + 12842, + 12869, + 12521, + 12426, + 12088, + 12015, + 11760, + 11410, + 10841, + 10485, + 9868, + 9425, + 8690, + 8069, + 7441, + 6843, + 6245, + 5798, + 5109, + 4664, + 4168, + 3715, + 3486, + 3139, + 2773, + 2639, + 2342, + 1999, + 1888, + 1695, + 1591, + 1429, + 1317, + 1171, + 1018, + 926, + 882, + 733, + 743, + 613, + 569, + 487, + 473, + 399, + 437, + 373, + 344, + 325, + 276, + 264, + 249, + 235, + 198, + 213, + 210, + 193, + 158, + 150, + 151, + 134, + 114, + 116, + 100, + 109, + 109, + 108, + 91, + 85, + 69, + 84, + 63, + 70, + 54, + 53, + 54, + 59, + 44, + 41, + 29, + 35, + 37, + 35, + 44, + 34, + 26, + 32, + 37, + 26, + 21, + 16, + 24, + 24, + 16, + 26, + 12, + 17, + 17, + 20, + 14, + 16, + 17, + 14, + 12, + 9, + 10, + 7, + 11, + 9, + 7, + 10, + 9, + 12, + 6, + 10, + 7, + 4, + 6, + 4, + 6, + 5, + 9, + 4, + 9, + 8, + 5, + 4, + 2, + 4, + 1, + 2, + 2, + 3, + 4, + 3, + 3, + 4, + 2, + 2, + 4, + 2, + 2, + 3, + 4, + 1, + 0, + 1, + 1, + 3, + 0, + 1, + 2, + 1, + 1, + 2, + 1, + 1, + 2, + 1, + 0, + 1, + 0, + 0, + 1, + 3, + 0, + 0, + 0, + 1, + 0, + 2, + 1, + 1, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 1, + 1, + 0, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1 + ] + }, + "nodata":null + } + ], + "eo:bands":[ + { + "name":"b1", + "description":"Gray" + } + ] + } +} diff --git a/src/histojam.md b/src/histojam.md new file mode 100644 index 0000000..bf57f2f --- /dev/null +++ b/src/histojam.md @@ -0,0 +1,225 @@ +--- +title: Gdaldem histojam +--- +# gdaldem histojam + +Making `gdaldem color-relief` color configuration files by hand can be a drag. Jam on some histogrammatic colormap file creation workflows instead. + +Calculate a histogram passing the `-hist` option to `gdalinfo` for your geotiff, write the output to a new json file, and upload it with the file selection option below. +```md +gdalinfo -json -hist file.tif > info.json +``` +Pick your colors and adjust the options to suit your needs, then copy the color configuration results to a text file for your next g'damn gdaldem color-relief attempt. +```js +const histFile = FileAttachment("data/histexample.json").json(); +``` +```js +const dicopalFile = FileAttachment("data/dicopaldata.json").json(); +``` + +```js +const jsonfile = view(Inputs.file({label: "gdalinfo JSON file", accept: ".json"})); +``` + +```js +const histData = jsonfile != null ? jsonfile.json() : histFile +``` +```js +const bandOpts = histData.bands.map((d, i)=> [d.band, i]) +``` +```js +const band = view(Inputs.select( + new Map(bandOpts), {label: "Band selection"} +)) +``` + +```js +const dataMin = histData.bands[band].min +const dataMax = histData.bands[band].max +const dataMean = histData.bands[band].mean +const histBinMin = histData.bands[band].histogram.min +const histBinMax = histData.bands[band].histogram.max +const histBinCount = histData.bands[band].histogram.count +const histBuckets = histData.bands[band].histogram.buckets +const minMax = [Math.floor(dataMin), Math.ceil(dataMax)] + +``` +```js +const cIn = Inputs.textarea({label: "Colors", rows: 2, value: 'purple,#ED2A24'}); +const cView = view(cIn) +``` + +```js +const colorVals = cView.split(',') +``` + +```js +const unoReverse = view(Inputs.button("Flip", {label:"Reverse colors",reduce:prestoChangeo})) + +``` +```js +const linspace = function (start, stop, nsteps) { + const delta = (stop - start) / (nsteps - 1); + return d3.range(nsteps).map((i) => start + i * delta); +} + +function prestoChangeo() { + cIn.value = colorVals.reverse().join(','); + cIn.dispatchEvent(new Event("input")); +}; + +function ramp(palette, n = 128) { + let color = d3.piecewise(d3.interpolateLab,palette); + const canvas = document.createElement("canvas"); + d3.select(canvas).attr("width", n).attr("height", 1); + const context = canvas.getContext("2d"); + canvas.style.margin = "2px"; + canvas.style.width = "30px"; + canvas.style.height = "30px"; + canvas.style.imageRendering = "-moz-crisp-edges"; + canvas.style.imageRendering = "pixelated"; + canvas.style.cursor = "pointer"; + canvas.onclick = (e) => { + if (!e.bubbles) return; + cIn.value = palette; + cIn.dispatchEvent(new Event("input")); + } + for (let i = 0; i < n; ++i) { + context.fillStyle = color(i / (n - 1)); + context.fillRect(i, 0, 1, 1); + } + return canvas; +} + +``` + +
+ Common colormaps + + +```js +const cmapFilter = view(Inputs.checkbox(["sequential","diverging"],{label:"Filter colormaps",value:["sequential","diverging"]})) +``` +```js +const selection = dicopalFile.filter(d => cmapFilter.includes(d.type)) +``` +```js +{ + const m = document.createElement("div") + selection.forEach(d => m.appendChild(ramp(d.colors))) + display(m); +} +``` + +
+ + +```js +const colorScheme = colorVals +``` +```js +const colorMin = view(Inputs.range(minMax, {label: "Color min", value: minMax[0], step: 0.01})); +const colorMax = view(Inputs.range(minMax, {label: "Color max", value: minMax[1], step: 0.01})); + +const interpolator = view(Inputs.select( + new Map([ + ["rgb", d3.interpolateRgb], + ["lab", d3.interpolateLab], + ["hcl", d3.interpolateHcl], + ["hcl-long", d3.interpolateHclLong] + ]), + { + label: "Color interpolator", + value: d3.interpolateHcl + })); + +``` +```js +const color = ({ + label: "min: " + colorMin + ", max: " + colorMax, + range: colorScheme, + domain: linspace(colorMin, colorMax, colorScheme.length), + legend: true, + type: "linear", + interpolate: interpolator, + clamp: true +}) +``` + +```js +const binW = (histBinMax - histBinMin) / histBinCount + +const binnedData = histBuckets.map((d,i)=>({ + value:d, + i:i, + x1: histBinMin + i * binW, + x2: histBinMin + (i+1) * binW + })) +``` + +```js +Plot.plot({ + y: { + clamp: false, + label:"count" + }, + color: { ...color, legend: false }, + marks: [ + Plot.rectY( + binnedData,{x1: "x1", x2: "x2", y2: "value", inset: 0.2, fill:"x1"} + ), + Plot.text( + [ + [colorMin, "color min"], + [colorMax, "color max"], + [dataMean ,"x̄"], + [dataMin,"min"], + [dataMax,"max"] + ], + { + x: "0", + y: 0, + text: (d) => d[1] + ": " + d[0], + textAnchor: "start", + dx: 5, + dy: -150 + } + ), + Plot.ruleX([colorMin, colorMax], { strokeDasharray: "2,2" }), + Plot.ruleX([dataMin, dataMax,dataMean], { strokeDasharray: "2,2", strokeOpacity:0.5 }) + ], + height: 200, + width: width, + marginRight: 100, + marginLeft: 50 +}) + +``` + + +```js +const steps = view(Inputs.range([2,256], {label: "Steps",transform: Math.sqrt, value: 11, step: 1})); +const alpha = view(Inputs.toggle({label:"Transparent nodata", value:true})) +``` + + +```js +const outVals = linspace(colorMin, colorMax,steps) +const outColors = d3.quantize(d3.piecewise(interpolator,colorVals),steps) + +let vrgba = outColors.map((d, i) => d3.format(".4")(outVals[i]) + ', ' + d.replace("rgb(", "").replace(")", ", 255")) +if (alpha) { + if (!vrgba.includes("nv, 0, 0, 0, 0")) { + vrgba.unshift("nv, 0, 0, 0, 0") + } +} else { + vrgba = vrgba.filter( d => d!=="nv, 0, 0, 0, 0") +} + +let formattedCmap = vrgba.join('\n') +document.getElementById("txtout").innerHTML = formattedCmap +``` + +

+
+Special thanks to the [gdal `gdaldem color-relief`](https://gdal.org/en/latest/programs/gdaldem.html#color-relief) function for bringing color to world, and to [dicopal.js](https://github.com/riatelab/dicopal.js) for making it easy to pull in a broad collection of color palettes.