A plugin for raster reprojection in d3

.block

height: 700
scrolling: yes
border: no

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.

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 });

}));

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>

ne_110m_land.topojson

preview.png

thumbnail.png

world.topo.bathy.200411.3x540x270.png