Skip to content

Instantly share code, notes, and snippets.

@mukhtyar

mukhtyar/.block

Forked from rasmuse/.block
Created October 26, 2017 19:26
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mukhtyar/ebb54bcbd571888c26a9957c16feb4bb to your computer and use it in GitHub Desktop.
Save mukhtyar/ebb54bcbd571888c26a9957c16feb4bb to your computer and use it in GitHub Desktop.
Download ZIP
A plugin for raster reprojection in d3
Raw
.block
height: 700
scrolling: yes
border: no
Raw
README.md

A plugin for raster reprojection in d3

Here is a demo of a simple d3 plugin to reproject (warp) a raster image using d3.

I was motivated to make it because I could not get gdalwarp to produce nice reprojections of rasters when the target projection had the antimeridian inside the picture, e.g. in a rotated orhographic projection (showing Earth's "back side"). The gdalwarp utility typically leaves a thin sliver unpainted close to the antimeridian, which shows up as a transparent line along the antimeridian.

How to use the plugin

The plugin behind this demo (see below, or in the GitHub repo) is only a proof of concept, not ready for production. But it illustrates that we can easily create a pretty neat API for raster reprojection, something like this:

var dstCanvas = d3.select('#warped').node();
var srcProj = d3.geoEquirectangular(),
    dstProj = d3.geoOrthographic(),

var world = {type: 'Sphere'},

// Fit the whole world inside the destination canvas
dstProj.fitSize([dstCanvas.width, dstCanvas.height]], world);

var warp = d3.geoWarp()
    .srcProj(srcProj)
    .dstProj(dstProj)
    .dstContext(dstContext);

var srcImg = d3.select('#src').node();

srcImg.onload = function () {
    // The source image is known to contain the whole world.
    // Of course, this can be hard coded if the image size is known.
    srcProj.fitSize([srcImg.width, srcImg.height], world);

    warp(srcImg);
}

srcImg.src = 'path/to/source/dataset.png';

Notes on implementation

The projection.fitExtent method (new in d3-geo v1.2.0) makes it easy to fit a d3 geoProjection to an image containing a map. The fitted projection object can then be used to sample the source image at some geographical coordinates without working out the conversion between image pixels and geographical coordinates.

No big assumptions are made about the source and destination projections. I only require that the destination projection is invertible. The source projection and the destination projection's inverse are called for each pixel in the destination raster. This is not fast, but it's simple :) As you can probably see from the poor animation frame rate on this page, the algorithm is not really useful for animated on-the-fly reprojection.

In this example, the source image is in equirectangular projection which is pretty typical for geographical raster datasets.

So far I only implemented nearest neighbor resampling, but other algorithms should be relatively straightforward to add.

Raw
d3-geo-warp.js
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-geo')) :
typeof define === 'function' && define.amd ? define(['exports', 'd3-geo'], factory) :
(factory((global.d3 = global.d3 || {}),global.d3));
}(this, function (exports,d3Geo) { 'use strict';
var NUM_COMPONENTS = 4; // RGBA
var SPHERE = {type: "Sphere"};
function linearCoords(img, col, row) {
return (row * img.width + col) * NUM_COMPONENTS;
}
function setPixel(img, col, row, val) {
var ix0 = linearCoords(img, col, row);
val.forEach(function (component, di) {
img.data[ix0 + di] = component;
});
}
function clamp(val, min, max) {
if (min > max) {
throw 'min > max'
}
return Math.max(min, Math.min(val, max));
}
function getPixel(imageData, col, row) {
var ix0 = linearCoords(imageData, col, row);
return imageData.data.slice(ix0, ix0 + NUM_COMPONENTS);
}
function isAlphaZero(imageData, col, row) {
return (imageData.data[linearCoords(imageData, col, row) + 3] == 0);
}
function interpolate(imageData, point) {
var x = point[0];
var y = point[1];
var col = clamp(Math.round(x), 0, imageData.width-1);
var row = clamp(Math.round(y), 0, imageData.height-1);
return getPixel(imageData, col, row);
}
function geoWarp() {
var srcProj, dstProj,
dstContext, dstCanvas,
maskObject;
function createCanvas(width, height) {
var canvas = document.createElement('canvas');
if (typeof width !== 'undefined') canvas.width = width;
if (typeof height !== 'undefined') canvas.height = height;
return canvas;
}
function warp(srcImage) {
var srcCanvas = createCanvas(srcImage.width, srcImage.height),
srcContext = srcCanvas.getContext('2d'),
maskData = makeMask();
srcContext.drawImage(srcImage, 0, 0);
srcImage = srcContext.getImageData(
0, 0, srcCanvas.width, srcCanvas.height);
var inverseProjection = function (point) {
return srcProj(dstProj.invert(point));
}
var dstImage = dstContext.getImageData(0, 0, dstCanvas.width, dstCanvas.height);
for (var row = 0; row < dstImage.height; row++) {
for (var col = 0; col < dstImage.width; col++) {
if (isAlphaZero(maskData, col, row)) continue;
var projectedPoint = [col+0.5, row+0.5],
inversePoint = inverseProjection(projectedPoint);
setPixel(dstImage, col, row, interpolate(srcImage, inversePoint));
}
}
dstContext.putImageData(dstImage, 0, 0);
}
warp.dstProj = function(_) {
if (!arguments.length) return dstProj;
dstProj = _;
return warp;
};
warp.srcProj = function(_) {
if (!arguments.length) return srcProj;
srcProj = _;
return warp;
};
warp.dstContext = function(_) {
if (!arguments.length) return dstContext;
dstContext = _;
dstCanvas = dstContext.canvas;
return warp;
};
function makeMask() {
var width = dstCanvas.width,
height = dstCanvas.height,
maskCanvas = createCanvas(),
context = maskCanvas.getContext('2d');
maskCanvas.width = width;
maskCanvas.height = height;
context.beginPath();
context.fillStyle = '#fff';
d3Geo.geoPath().projection(dstProj).context(context)(maskObject);
context.closePath();
context.fill();
return context.getImageData(0, 0, width, height);
}
warp.maskObject = function (_) {
if (!arguments.length) return maskObject;
maskObject = _;
return warp;
}
warp.createCanvas = function (_) {
if (!arguments.length) return createCanvas;
createCanvas = _;
return warp;
}
return warp.maskObject(SPHERE);
}
exports.geoWarp = geoWarp;
Object.defineProperty(exports, '__esModule', { value: true });
}));
Raw
index.html
<!DOCTYPE html>
<meta charset="utf-8">
<style>
body {
font-family: Arial;
}
</style>
<h2>Source image</h2>
<img id="src"/>
<h2>Reprojected image (with Natural Earth land polygons)</h2>
<canvas id="warped" alt="Reprojected map"></canvas>
<svg id="land"></svg>
<script src="//d3js.org/d3.v4.min.js"></script>
<script src="//d3js.org/topojson.v1.min.js"></script>
<script src="d3-geo-warp.js"></script>
<script>
var width = 500,
height = 500,
world = {type: 'Sphere'},
latSpeed = -0.02,
longSpeed = 0.01,
rasterWarper,
drawLand;
var dstProj = d3.geoOrthographic().fitExtent([[20, 20], [width-40, height-40]], world);
var dstCanvas = d3.select('#warped').node();
dstCanvas.width = width;
dstCanvas.height = height;
var dstContext = dstCanvas.getContext('2d');
dstContext.fillStyle = '#444'
dstContext.fillRect(0, 0, width, height);
path = d3.geoPath().projection(dstProj).context(dstContext);
var srcImg = d3.select('#src').node();
srcImg.onload = function () {
var srcProj = d3.geoEquirectangular()
.fitSize([srcImg.width, srcImg.height], world);
rasterWarper = d3.geoWarp()
.srcProj(srcProj)
.dstProj(dstProj)
.dstContext(dstContext);
d3.json("ne_110m_land.topojson", function(error, earth) {
if (error) return console.error(error);
var land = topojson.feature(earth, earth.objects.ne_110m_land);
drawLand = function() {
dstContext.beginPath();
dstContext.strokeStyle = 'rgb(200, 0, 0)';
dstContext.fillStyle = 'rgba(200, 0, 0, .2)';
path(land);
dstContext.closePath();
dstContext.stroke();
dstContext.fill();
}
});
}
srcImg.src = 'world.topo.bathy.200411.3x540x270.png'
d3.timer(function(elapsed) {
dstProj.rotate([longSpeed * elapsed, latSpeed * elapsed]);
if (!!rasterWarper) rasterWarper(srcImg);
if (!!drawLand) drawLand();
});
</script>
Display the source blob
Display the rendered blob
Raw
ne_110m_land.topojson
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment