Re-project Web Mercator tiles to arbitrary projection with D3?

This answer builds upon:

These three resources also build on one another. Understanding these three examples will help in terms of understanding what occurs in my example below.

The answer also uses, as a base, my slow moving, ongoing attempt to build a tile library.

The goal of this answer is not to present a finalized resource, but a rough demonstration of how one might be put together along with associated information. The answer will also evolve as I further my thought on the problem.

Web Mercator Tiles

A Mercator map that stretches over 360 degrees of longitude and ~170 degrees of latitude (+/- 85 degrees) will fill a square (Going beyond 85 degrees in latitude causes distortion to get out of hand, and inclusion of the poles isn't advisable as the poles are at +/-infinity on the projected plane).

This square of most of the world is, for a web mapping service (with mercator tiles), zoom level 0. The map is 2^0 squares across and 2^0 squares high.

If we divide that square into a grid of two squares by two squares we have zoom level 1. The map is 2^1 by 2^1 squares.

Consequently, the zoom level dictates how many squares across and high the map is: 2^zoomLevel. If each square is the same size in pixels, then each increase of the zoom level by one increases the pixel width of the world by 2x.

Luckily for us, there is no land north of ~85 degrees and it isn't often we want to show Antarctica, so this square is suitable for most web mapping applications. However, it means if we are reprojecting web mercator tiles to anything that shows above these latitudes, we'll have a bald spot:

Wᴇʙ-Mᴇʀᴄᴀᴛᴏʀ ᴛɪʟᴇs ʀᴇᴘʀᴏᴊᴇᴄᴛᴇᴅ ғᴏʀ ᴀ Mᴇʀᴄᴀᴛᴏʀ ᴘʀᴏᴊᴇᴄᴛɪᴏɴ ᴛʜᴀᴛ ʜᴀs ʙᴇᴇɴ ʀᴏᴛᴀᴛᴇᴅ ᴛᴏ sʜᴏᴡ ᴛʜᴇ Nᴏʀᴛʜ Pᴏʟᴇ.

Lastly, web mercator tiles are rendered in a projected space that is predictable and regular relative to the tiles. If we re-project the tiles, we may be enlarging or shrinking the projected space for each tile, we should be mindful of this. In the image above, the tiles around the North Pole are reprojected as much smaller than those further south. Tiles will not necessarily be uniformly sized after projection.

Reprojection and Resampling Plate Carree

The biggest challenge to reprojecting web service tiles is time - and not just time spent understanding projections and reading answers like this.

Projection functions are complex time consuming operations that must be done on each and every pixel that is rendered. All d3 examples I have seen use the process (or a close variant) as seen here for the actual reprojection and resampling. This example will only work if the original image is projected with Plate Carree. The process is as follows:

Create a blank new image.
For each pixel in the new image, take its location in pixels and invert it (using the desired projection) to get a latitude and longitude.
Determine which pixel in the original image overlaps that longitude and latitude.
Take the information from that pixel in the original image and assign it to the appropriate pixel in the new image (the pixel in step 2)

When the original image uses a Plate Carree projection we don't need a d3-geoProjection, the relationship is linear between projected and unprojected coordinates. For example: if the image is 180 pixels high, each pixel represents 1 degree of latitude. This means step 3 doesn't take very long as compared with step 2 and projection.invert(). Here's Mike's function for step 3:

var q = ((90 - φ) / 180 * dy | 0) * dx + ((180 + λ) / 360 * dx | 0) << 2;

The time required for step 2 is related to the projection used for the reprojected image. All examples I've seen use d3.geoProjection.invert() for step two in the above list - taking a pixel location in a new image and finding out its latitude and longitude. Not all projections are born equal. Cylindrical projections will generally outperform conical projections, and conical projections will generally out perform azimuthal projections. I've also seen some odd speed discrepencies between d3v4 and d3v5 in terms of projection.invert() time:

Tɪᴍᴇ ᴛᴏ ᴜɴᴘʀᴏᴊᴇᴄᴛ ᴀ ᴘᴏɪɴᴛ ᴏɴ ᴀ ᴍᴀᴘ ᴡɪᴛʜ D3 (ᴄᴏɴᴠᴇʀᴛ ғʀᴏᴍ ᴘɪxᴇʟs ᴛᴏ ʟᴀᴛ/ʟᴏɴɢ). Iᴛ ɪs ᴜɴᴄʟᴇᴀʀ ᴡʜʏ D3ᴠ4 ᴡᴏᴜʟᴅ ʙᴇ ғᴀsᴛᴇʀ.

And for completeness, here's a wider range of projections found in d3-geo-projections:

Cᴏᴍᴘᴀʀɪɴɢ ᴛɪᴍᴇ ғᴏʀ ᴘʀᴏᴊᴇᴄᴛɪᴏɴ.ɪɴᴠᴇʀᴛ ғᴏʀ ᴀᴅᴅɪᴛɪᴏɴᴀʟ ᴘʀᴏᴊᴇᴄᴛɪᴏɴs

Notes

This approach can run into problems described here, these problems can be addressed through the answers there - but different solutions will cost additional time to process.
This approach uses a nearest neighbour approach - which might lead to quality issues. More advanced sampling such as bilinear or cubic will add time to the process but may result in more desirable image.
If the base image has text, text may be rotated or otherwise manipulated to make it less or unreadable.
Mike's example is for a single image, for tiles, the process is altered to some degree, we are now creating multiple images, and this requires knowing each original tile's bounds and each reprojected tile's bounds in degrees as well as in tile units for the former and pixels for the latter - minor details.

Reprojection and Resampling Web Mercator

When I started looking at this question I looked at Alan McConchie's solution/example as a reference. It took a while to notice, but step 3 in this example (and I believe Jason Davies' work too) doesn't account for a web Mercator tile in the resampling - only in determining the tile bounds. But, the relationship between pixels on the y axis is no longer linear as it was in the Plate Carree.

This means tiles are placed in the right places, but the sampling treats the y axis as linear within each tile. This distortion was most visible when showing the whole world with a low tile zoom level (midway down/up the tile), and may be what Alan speaks of when he mentions odd compression.

The solution is to properly project the latitude for each latitude/longitude pair in step 3 above. This adds time, always more time - the function involves Math.atan and Math.exp, but the difference shouldn't be too bad. In both Alan and Jason's work this is done with a simple formula (but only used for tile bounds, not each pixel):

Math.atan(Math.exp(-y * Math.PI / 180)) * 360 / Math.PI - 90;

In my example below, I've just used d3.geoMercator() to make for a more clear scaling factor, using the projection includes one extra operation to transform the x coordinate.

Otherwise, the 4 step process remains the same.

Finding the right tiles

I've only seen one clean approach to finding what tiles to show, Jason Davies' d3.quadTile, seen here. I believe that Alan McConchie uses an unminified version that may be otherwise altered. There is also this github repository for another version of d3.quadTiles, which is very similar.

For McConchie/Davies, d3.quadTile will, given a projection with a clip extent (not clip angle) and a tile depth, pull out all the tiles that intersect the view extent.

In Alan McConchie's solution/example, the zoom level is based off of the projection scale - but this is not necessarily the wisest: each projection has different scaling factors, a scale of 100 on one scale will show a different extent than a scale of a 100 on another. Also, the relationship between scale values and map size in a cylindrical projection may be linear, while non-cylindrical projections may have non linear relationships between map size and scale.

I've modified this approach a bit - I use a scale factor to determine an initial tile depth, and then reduce that tile depth if the tile count returned by d3.quadTile exceeds a certain number:

geoTile.tileDepth = function(z) {
    // rough starting value, needs improvement:
    var a = [w/2-1,h/2]; // points in pixels
    var b = [w/2+1,h/2];
    var dx = d3.geoDistance(p.invert(a), p.invert(b)) ; // distance between in radians      
    var scale = 2/dx*tk;
    var z = Math.max(Math.log(scale) / Math.LN2 - 8, 2);
    z = Math.min(z,15) | 0;

    // Refine:
    var maxTiles = w*h/256/128;
    var e = p.clipExtent();
    p.clipExtent([[0,0],[w,h]])
    while(d3.quadTiles(p, z).length > maxTiles) {
        z--;
    }
    p.clipExtent(e);

    return z;
}

Then, using d3.quadTile I pull out the relevant tiles:

geoTile.tiles = function() {
    // Use Jason Davies' quad tree method to find out what tiles intercept the viewport:
    var z = geoTile.tileDepth();
    var e = p.clipExtent(); // store and put back after.

    p.clipExtent([[-1,-1],[w+1,h+1]]) // screen + 1 pixel margin on outside.
    var set = d3.quadTiles(p, Math.max(z0,Math.min(z,z1))); // Get array detailing tiles
    p.clipExtent(e);

    return set;
}

At first I thought that pulling tiles from multiple zoom depths (to account for discrepencies in size of reprojected tiles) would be ideal: but this will run into problems with things like line thickness in the raster as well as discontinuous annotations.

Adopting Jason and Alan's Work

I take the tile set generated above with geoTile.tiles() and run it through an enter/update/exit cycle using the tile coordinate (in tile coordinates, row, column, zoom depth) as a key, appending image elements to a parent g or svg. When loading images, once an image loads, we call a onload function to do the actual reprojection. This is largely unchanged from Jason and Alan, I've addressed the following challenges I saw in this code:

Resampling did not account for web mercator (mentioned above)
Tile depth was not selected well (mentioned above)
Tiles were reprojected as canvas-es placed in a div rather than an SVG - creating two parent containers, one for each type of feature: tile or vector.

I believe my example, with very minor tweaks, has addressed these. I've also added some more extensive comments for viewing:

    function onload(d, that) { // d is datum, that is image element.

        // Create and fill a canvas to work with.
        var mercatorCanvas = d3.create("canvas")
          .attr("width",tileWidth)
          .attr("height",tileHeight);
        var mercatorContext = mercatorCanvas.node().getContext("2d");           
        mercatorContext.drawImage(d.image, 0, 0, tileWidth, tileHeight); // move the source tile to a canvas.

        //
        var k = d.key; // the tile address.
        var tilesAcross = 1 << k[2]; // how many tiles is the map across at a given tile's zoom depth?

        // Reference projection:
        var webMercator = d3.geoMercator()
          .scale(tilesAcross/Math.PI/2) // reference projection fill square tilesAcross units wide/high.
          .translate([0,0])
          .center([0,0])

        // Reprojected tile boundaries in pixels.           
        var reprojectedTileBounds = path.bounds(d),
        x0 = reprojectedTileBounds[0][0] | 0,
        y0 = reprojectedTileBounds[0][1] | 0,
        x1 = (reprojectedTileBounds[1][0] + 1) | 0,
        y1 = (reprojectedTileBounds[1][1] + 1) | 0;

        // Get the tile bounds:
        // Tile bounds in latitude/longitude:
        var λ0 = k[0] / tilesAcross * 360 - 180,                     // left        
        λ1 = (k[0] + 1) / tilesAcross * 360 - 180,                   // right
        φ1 = webMercator.invert([0,(k[1] - tilesAcross/2) ])[1],     // top
        φ0 = webMercator.invert([0,(k[1] + 1 - tilesAcross/2) ])[1]; // bottom.             

        // Create a new canvas to hold the what will become the reprojected tile.
        var newCanvas = d3.create("canvas").node();

        newCanvas.width = x1 - x0,      // pixel width of reprojected tile.
        newCanvas.height = y1 - y0;     // pixel height of reprojected tile.
        var newContext = newCanvas.getContext("2d");    

        if (newCanvas.width && newCanvas.height) {
            var sourceData = mercatorContext.getImageData(0, 0, tileWidth, tileHeight).data,
                target = newContext.createImageData(newCanvas.width, newCanvas.height),
                targetData = target.data;

            // For every pixel in the reprojected tile's bounding box:
            for (var y = y0, i = -1; y < y1; ++y) {
              for (var x = x0; x < x1; ++x) {
                // Invert a pixel in the new tile to find out it's lat long
                var pt = p.invert([x, y]), λ = pt[0], φ = pt[1];

                // Make sure it falls in the bounds:
                if (λ > λ1 || λ < λ0 || φ > φ1 || φ < φ0) { i += 4; targetData[i] = 0; continue; }  
                    // Find out what pixel in the source tile matches the destination tile:
                    var top = (((tilesAcross + webMercator([0,φ])[1]) * tileHeight | 0) % 256  | 0) * tileWidth;
                    var q = (((λ - λ0) / (λ1 - λ0) * tileWidth | 0) + (top)) * 4;

                    // Take the data from a pixel in the source tile and assign it to a pixel in the new tile.
                    targetData[++i] = sourceData[q];
                    targetData[++i] = sourceData[++q];
                    targetData[++i] = sourceData[++q];
                    targetData[++i] = 255;
              }
            }
            // Draw the image.
            if(target) newContext.putImageData(target, 0, 0);
        }

        // Add the data to the image in the SVG:
        d3.select(that)
          .attr("xlink:href", newCanvas.toDataURL()) // convert to a dataURL so that we can embed within the SVG.
          .attr("x", x0)
          .attr("width", newCanvas.width)
          .attr("height",newCanvas.height)
          .attr("y", y0);
    }

Placing it within a larger structure.

A regular tile map with overlain features has a few coordinate systems:

Tile units (3D), which mark column, row, and zoom level of each tile (x,y,z respectively)
Geographic coordinates (3D), which mark latitude and longitude of a point on a three dimensional sphere.
Zoom units (3D), which keep track of zoom translate (x,y) and zoom scale (k).
projected units (2D), pixels units which latitude and longitude are projected to.

The goal of any slippy map is to unify these coordinates in a usable system.

When we reproject the tiles, we need to add a coordinate space:

The(/a) tile set projection.

I felt the examples were not particularily clear in how they tied together all the coordinate systems. So, I've placed the above methods, as you might have seen, in a geoTile object that is taken from a personal project for a tile library. The goal of this is for a bit smoother coordination of different units. I am not trying to plug it, it's still in development in any event (just too busy to really finish it); however, I will see if time affords me the opportunity to rig an example up with d3-tile.

Challenges moving forward

Zoom speed and responsiveness is the greatest challenge I see. To address this I've set the zoom function to trigger on zoom end - this is most noticable on pan events, as normally a pan triggers the zoom function continuously through the pan, this could be addressed by translating existing images. The most reliable way to use this, however, is on a static map. Implementing a translate for already drawn images would be ideal for pan events, rather than resampling as currently.

Animating such a map is probably not possible.

There is probably room for optimizing the calculations that turn pixel to latitude, but this might be difficult.

Examples

Unfortunately, the code is too much for a snippet, so I've made a few bl.ocks to demonstrate.

These have had only minimal testing, if I manage to finish the foundational tile library, I'll fork it for this purpose, in the meantime it should suffice as an example. The guts of the code is found in geoTile.tile() in the d3-reprojectSlippy.js file, which contains the enter/update/exit cycle (fairly basic) and the onload function described above. As I do work on the tiles a fair bit on the side I'll keep this answer updated.

The alternative

Reprojecting tiles is cumbersome and time consuming. An alternative, if possible, would be to produce a tile set in the desired projection. This has been done with OSM tiles, but is also cumbersome and time consuming - just for the map maker, not the browser.

TL;DR

Reprojected Mercator tiles take time, you should read the above.