2b:null 2c:["$","title",null,{"children":"Raster reprojection | Jake Low | Observable"}] 31:T5fc,md` The map above is a Waterman-style polyhedral projection, but the source data is an image I borrowed from Wikipedia (shown below) that uses the Molleweide projection. In order to reproject, I defined two projections — one that described the source image, and one that described my desired projection. Then, I used \`geoWarp\` (a function I wrote, see the source below) to translate between the projections. The technique is based on [some work that Mike Bostock did](https://bl.ocks.org/mbostock/4329423) a while back. Basically, I iterate over each pixel in the target image (a blank canvas where I'd like to render the Waterman map). For each pixel, I invert from canvas coordinates to spherical coordinates using the Waterman projection's \`invert\` function, and then project these spherical coordinates into coordinates on the Molleweide map. Whatever color I find at that pixel, I paint onto my output canvas. This is a very rough draft of what a module for doing raster reprojection might support. Performance isn't as good as it probably could be, and the only resampling algorithm available is "nearest-neighbor" which is fast but not very pretty. ~~It also looks pixelated on high-DPI screens.~~ Update: High-DPI mode now works, but only because we assume that the canvas context's \`scale\` is equal to the \`devicePixelRatio\`. This is true in most cases on Observable, because it's how \`DOM.context2d()\` works internally, but ideally it would be possible to get the current scale of the context. `32:Tfde,function geoWarp(targetProjection, context) { /* * geoWarp is designed to have a similar interface to geoPath -- it takes a * projection and a context, and returns a function. The function can be used * to paint an image onto a canvas context, much like how a path function can * be used to paint GeoJSON geometries. * * The code has been instrumented to log some timing information so I could * better understand what the performance bottlenecks are. You can ignore the * various `performance.now()` calls and updates to the `stats` object; they are * not part of the core algorithm. */ return function(sourceProjection, image) { const start = performance.now(); const stats = { pixels: { total: 0, sourced: 0, cheated: 0 }, setup: 0, cheat: 0, invert: 0, project: 0, copy: 0, paint: 0, unaccounted: 0, }; let checkpoint = performance.now(); const bounds = d3.geoPath(targetProjection).bounds(sphere); const sourceContext = DOM.canvas(image.width, image.height).getContext("2d"); sourceContext.drawImage(image, 0, 0, image.width, image.height); const source = sourceContext.getImageData(0, 0, image.width, image.height) const sourceData = source.data; const dpr = window.devicePixelRatio || 1; // TODO the dimensions of the target are hardcoded; can we get them from the context instead? const target = context.getImageData(0, 0, width * dpr, height * dpr); const targetData = target.data; // make a "cheat sheet" by projecting the sphere and filling it; we can use the resulting image data // to decide whether to invert or not. const cheatSheetContext = DOM.context2d(width, height); cheatSheetContext.beginPath(); d3.geoPath(targetProjection, cheatSheetContext)(sphere); cheatSheetContext.fill(); const cheatSheetImageData = cheatSheetContext.getImageData(0, 0, width * dpr, height * dpr); const cheatSheetData = cheatSheetImageData.data; stats.pixels.total = targetData.length / 4; stats.setup = performance.now() - checkpoint; for (var y = 0; y < target.height; y += 1) { for (var x = 0; x < target.width; x += 1) { checkpoint = performance.now(); let i = (y * target.width + x) * 4; // compute `i` early so we can use it for the cheat sheet const cheatPixelAlpha = cheatSheetData[i + 3]; stats.cheat += performance.now() - checkpoint; // if any cheat pixel isn't black, it hasn't been filled when we projected the sphere, so skip it. if (cheatPixelAlpha === 0) { stats.pixels.cheated += 1; continue; } checkpoint = performance.now(); const p = [x / dpr, y / dpr]; const s = targetProjection.invert(p); stats.invert += performance.now() - checkpoint; if (!s) continue; const [λ, φ] = s; if (λ > 180 || λ < -180 || φ > 90 || φ < -90) continue; stats.pixels.sourced += 1; checkpoint = performance.now(); const q = sourceProjection([λ, φ]); const u = Math.round(q[0]), v = Math.round(q[1]); let now = performance.now(); stats.project += now - checkpoint; checkpoint = now; // let i = (y * target.width + x) * 4; let j = (v * image.width + u) * 4; targetData[i] = sourceData[j]; targetData[++i] = sourceData[++j]; targetData[++i] = sourceData[++j]; targetData[++i] = 255; stats.copy += performance.now() - checkpoint; } } checkpoint = performance.now(); context.putImageData(target, 0, 0); const end = performance.now(); stats.paint += end - checkpoint; stats.unaccounted = (end - start) - (stats.setup + stats.invert + stats.project + stats.copy + stats.paint); mutable stats = stats; } }33:[] 2f:{"id":"982b422280f470c1","slug":"raster-reprojection","title":"Raster reprojection","thumbnail":"9fa0f026c50db9a21f7f5f1d356da46aa648b2314f1227e9b1455eeaa797a3a1","access_level":"public","create_time":"$D2018-10-18T02:32:26.061Z","update_time":"$D2019-11-01T00:58:36.195Z","view_time":null,"first_public_time":"$D2018-10-19T03:05:12.195Z","fork":false,"stars":13,"starred":false,"owner":{"type":"team","id":"cbb2b32473d3e648","login":"jake-low","name":"Jake Low","avatar_url":"https://avatars-next.observableusercontent.com/fcff7503af57e71a7006d2543ad38643918ba16eb09e01f116b42d94a497c0a0","tier":"starter_2024"},"version":593,"publish_version":593,"role":null,"authored":false,"authors":[{"type":"individual","id":"e30eca3d7ea0ab2c","login":"jake-low","name":"Jake Low","avatar_url":"https://avatars-next.observableusercontent.com/fcff7503af57e71a7006d2543ad38643918ba16eb09e01f116b42d94a497c0a0","tier":null}],"members":[{"type":"individual","id":"e30eca3d7ea0ab2c","login":"jake-low","name":"Jake Low","avatar_url":"https://avatars-next.observableusercontent.com/fcff7503af57e71a7006d2543ad38643918ba16eb09e01f116b42d94a497c0a0","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":11,"collection":null,"type":"notebook","body":{"title":"Raster reprojection","cells":[{"id":0,"value":"md`# Raster reprojection`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":384,"value":"md`\n_Raster reprojection_ is the act of transforming a map encoded in a raster image (such as a JPEG or PNG) to a different map projection. This is useful for showing satellite imagery, or for visualising interesting datasets (like the one used below) using a different projection than the original authors used when they created the data. \n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":50,"value":"target = {\n const context = DOM.context2d(width, height);\n \n const warp = geoWarp(waterman, context);\n warp(molleweide, image);\n \n const path = d3.geoPath(waterman, context);\n context.lineWidth = lineWidth;\n context.strokeStyle = \"white\";\n context.beginPath(), path(land), context.stroke();\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":438,"value":"md`This map shows sea surface temperatures. The content comes from a source image that uses an entirely different projection.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":397,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":76,"value":"source = {\n const context = DOM.context2d(image.width, image.height);\n const path = d3.geoPath(molleweide, context);\n \n context.drawImage(image, 0, 0, image.width, image.height);\n \n context.lineWidth = lineWidth;\n context.strokeStyle = \"white\";\n \n context.beginPath(), path(land), context.stroke();\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":429,"value":"md`The source image [from Wikipedia (by user Giorgiogp2, CC-BY-SA 3.0)](https://en.wikipedia.org/wiki/Sea_surface_temperature#/media/File:MODIS_sst.png) shows sea surface temperature on a Molleweide projection. The white landmass borders are drawn on top with D3, demonstrating that the source projection is aligned correctly.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":105,"value":"$32","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":548,"value":"mutable stats = null","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":322,"value":"waterman = d3.geoPolyhedralWaterman()\n .fitSize([width, height], sphere)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":65,"value":"molleweide = d3.geoMollweide()\n .fitExtent([[0.005 * image.width, 0], [0.995 * image.width, 0.906 * image.height]], sphere)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":516,"value":"width = image.width","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":389,"value":"height = width / 1.6","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":520,"value":"lineWidth = 0.4 + width / 1000","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":24,"value":"url = \"https://upload.wikimedia.org/wikipedia/commons/f/fd/MODIS_sst.png\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":22,"value":"land = topojson.feature(world, world.objects.land)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":20,"value":"sphere = ({ type: \"Sphere\" })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":4,"value":"image = d3.image(url, { crossOrigin: \"anonymous\" })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":14,"value":"world = d3.json(\"https://unpkg.com/world-atlas@1/world/110m.json\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":18,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":7,"value":"d3 = require(\"d3-fetch\", \"d3-geo\", \"d3-geo-projection\", \"d3-geo-polygon\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":518,"value":"stylesheet = html`\n(restyle images to be no wider than the main text)\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W33","stdlib":"1","titling":"auto"}} 30:"$2f:body"