javascript - 将Web Mercator磁贴重新投影到D3的任意投影？

Question

杰森·戴维斯(Jason Davies)用Reprojected Raster Tiles引爆我们已经过去了数年，因为Mapbox阻止了他的网站， map 停止工作了，但是Mollweide Watercolour和Interrupted Goode Raster仍然是不错的演示。

现在在Observable HQ上，我看到了最新的d3-geo-projection和d3-tile的文档，但是没有关于如何执行Jason所做的工作的现代示例:重新投影标准Mercator贴图集。

我如何才能使d3平铺变形到新的投影？

Answer 1

该答案基于:

Mike Bostock's Raster Reprojection

Jason Davies' Quad Tiles和

Alan McConchie's Raster Reprojection

这三种资源也相互补充。了解这三个示例将有助于理解下面的示例中发生的情况。

答案还以我缓慢的，持续的尝试来构建tile library为基础。

该答案的目的不是提供最终的资源，而是粗略演示如何将其与相关信息放在一起。当我进一步思考该问题时，答案也将不断演变。

网络墨卡托瓷砖

延伸超过360度经度和〜170度纬度(+/- 85度)的Mercator map 将填满一个正方形(超过85度纬度会导致变形失控，并且不包含极点)由于极点在投影平面上为+/-无穷大，因此建议使用。

对于Web映射服务(带有墨卡托磁贴)，这个世界上大多数地方的缩放级别为0。 map 的宽度为2 ^ 0平方，高度为2 ^ 0平方。

如果将那个正方形分成两个正方形乘以两个正方形的网格，则缩放级别为1。 map 为2 ^ 1 x 2 ^ 1正方形。

因此，缩放级别决定了 map 上横跨多少个正方形:2 ^ zoomLevel。如果每个正方形的像素大小相同，则缩放级别每增加1，世界的像素宽度就会增加2倍。

对我们来说幸运的是，在〜85度以北没有土地，而且我们并不经常想要展示南极洲，因此这个正方形适合大多数网络 map 应用。但是，这意味着如果将网络墨卡托磁贴重新投影到这些纬度以上显示的任何内容，我们将有一个秃头:



Wᴇʙ-MᴇʀᴄᴀᴛᴏʀᴛɪʟᴇsᴀᴇʀᴄᴀᴛᴏʀʜᴀᴛʜᴇᴛʜᴇNᴛʜᴇPᴏʟᴇ。

最后，网络墨卡托瓷砖在一个相对于瓷砖而言可预测且规则的投影空间中进行渲染。如果我们重新投影图块，则可能会放大或缩小每个图块的投影空间，我们应该注意这一点。在上图中，北极周围的瓷砖被重新投影，比南端的瓷砖要小得多。投影后，瓷砖的尺寸不一定会统一。

重投影和重采样平板支架

重新投影Web服务磁贴的最大挑战是时间-而不仅仅是花在理解投影和阅读答案上的时间。

投影功能是耗时的复杂操作，必须对每个渲染的像素都必须完成。我看到的所有d3示例都使用here所示的过程(或近似变体)进行实际的重新投影和重采样。仅当原始图像使用Plate Carree投影时，此示例才有效。过程如下:

创建空白新图像。

对于新图像中的每个像素，以像素为单位定位并反转(使用所需的投影)以获取经度和纬度。

确定原始图像中的哪个像素与该经度和纬度重叠。

从原始图像中的该像素获取信息，并将其分配给新图像中的适当像素(步骤2中的像素)。

当原始图像使用Plate Carree投影时，我们不需要d3-geoProjection，投影坐标和非投影坐标之间的关系是线性的。例如:如果图像的高度为180像素，则每个像素代表1个纬度。这意味着与步骤2和projection.invert()相比，步骤3不需要花费很长时间。这是步骤3的Mike的功能:

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

步骤2所需的时间与用于重投影图像的投影有关。我见过的所有示例在上面的列表的第二步中都使用d3.geoProjection.invert()-在新图像中获取像素位置并找出其经度和纬度。并非所有的预测都是天生的。圆柱投影通常会胜过圆锥投影，而圆锥投影通常会表现出方位 Angular 投影。我还看到了在projection.invert()时间方面，d3v4和d3v5之间存在一些奇怪的速度差异:



ᴍᴇᴛʜᴛʜᴛʜᴛʜD3(ᴄᴏɴᴠᴇʀᴛᴇʟᴇʟʟᴏɴɢ/ʟᴏɴɢ)。 ᴛD3ᴠ4ғᴀsᴛᴇʀ。

为了完整起见，在d3-geo投影中可以找到更广泛的投影:



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

笔记

这种方法可能会遇到here描述的问题，可以通过那里的答案解决这些问题-但是不同的解决方案将花费额外的处理时间。

此方法使用最近邻居方法-可能会导致质量问题。更高级的采样(例如双线性或三次采样)会增加处理时间，但可能会产生更理想的图像。

如果基础图像包含文本，则可以旋转文本或通过其他方式处理文本以使其变得更少或不可读。

Mike的示例是针对单个图像的，对于图块，此过程已进行了一定程度的更改，我们现在正在创建多个图像，这需要了解每个原始图块的边界以及每个重新投影的图块的边界(以度为单位)以及图块单位前者和后者的像素-次要细节。

重新投影和重新采样Web墨卡托

当我开始研究这个问题时，我以Alan McConchie的solution/example为引用。花了一段时间才注意到，但是此示例中的第3步(我也相信Jason Davies的工作)在重新采样中并未考虑Web Mercator磁贴-仅用于确定磁贴边界。但是，y轴上像素之间的关系不再像Plate Carree中那样线性。

这意味着将图块放置在正确的位置，但是采样将每个图块中的y轴视为线性。当以低图块缩放级别(在图块的中/下)显示整个世界时，这种失真最明显，并且可能是艾伦在提到奇数压缩时所说的。

解决方案是在上述步骤3中为每个纬度/经度对正确投影纬度。这会增加时间，总会增加时间-该函数涉及Math.atan和Math.exp，但差别不应该太糟。在Alan和Jason的作品中，这都是通过一个简单的公式完成的(但仅用于图块范围，而不是每个像素):

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

在下面的示例中，我只是使用d3.geoMercator()来使比例因子更清晰，使用投影包括一项额外的操作来转换x坐标。

否则，四步过程保持不变。

寻找合适的瓷砖

我只看过一种干净的方法来查找要显示的图块，Jason Davies的d3.quadTile见here。我相信Alan McConchie使用的是未缩小的版本，否则可能会被更改。还有另一个版本d3.quadTiles的this github存储库，这非常相似。

对于McConchie/Davies，d3.quadTile将在给定具有裁剪范围(而非裁剪 Angular )和贴图深度的投影的情况下，拉出与 View 范围相交的所有贴图。

在艾伦·麦康琪(Alan McConchie)的solution/example中，缩放级别基于投影比例-但这不一定是最明智的:每个投影都有不同的比例因子，一个比例上的100比例将显示一个比100上的比例不同的范围。其他。而且，圆柱投影中的比例值和 map 尺寸之间的关系可以是线性的，而非圆柱投影可以在 map 尺寸和比例之间具有非线性的关系。

我已经对该方法进行了一些修改-我使用比例因子来确定初始图块深度，然后如果d3.quadTile返回的图块计数超过一定数量，则减小该图块深度:

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

然后，使用d3.quadTile提取相关的图块:

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

最初，我认为从多个缩放深度拖动图块(以解决重新投影的图块的大小差异)将是理想的:但这会遇到诸如栅格中的线宽以及不连续注释之类的问题。

采纳杰森和艾伦的作品

我将上面使用geoTile.tiles()生成的图块集并以图块坐标(在图块坐标，行，列，缩放深度中)作为键，通过进入/更新/退出循环运行它，并将image元素附加到父g或svg。加载图像时，一旦加载图像，我们就调用onload函数进行实际的重新投影。这与Jason和Alan基本上没有变化，我已经解决了我在这段代码中看到的以下挑战:

重采样未解决网络墨卡托操作(如上所述)

未正确选择图块深度(如上所述)

将图块重新投影为放置在div中而不是SVG中的 Canvas ，从而创建了两个父容器，每个父容器对应一种要素类型:图块或矢量。

我相信我的示例经过很小的调整就解决了这些问题。我还添加了一些更广泛的评论以供查看:

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

将其放置在较大的结构中。

具有覆盖特征的常规图块 map 具有一些坐标系:

图块单位(3D)，用于标记每个图块的列，行和缩放级别(分别为x，y，z)

地理坐标(3D)，用于标记三维球体上某个点的纬度和经度。

缩放单位(3D)，用于跟踪缩放平移(x，y)和缩放比例(k)。

投影单位(2D)，即投影到纬度和经度的像素单位。

任何草率的贴图的目标是在可用系统中统一这些坐标。

重新投影图块时，我们需要添加一个坐标空间:

(/a)瓷砖设置投影。

我觉得这些示例在如何将所有坐标系联系在一起方面并不是特别清楚。因此，正如您可能已经看到的，我将上述方法放置在geoTile对象中，该对象是从tile library的个人项目中获取的。这样做的目的是使不同单位之间的协调更加顺畅。我并不是要插入它，无论如何它仍在开发中(只是太忙而无法完成它)。但是，我将看看时间是否能让我有机会使用d3-tile搭建示例。

前进的挑战

变焦速度和响应速度是我看到的最大挑战。为了解决这个问题，我将缩放功能设置为在缩放结束时触发-在平移事件中最明显，因为通常平移会通过平移连续触发缩放功能，这可以通过转换现有图像来解决。但是，使用此方法最可靠的方法是在静态 map 上。为平移事件而不是当前重采样，为已经绘制的图像实现转换将是理想的选择。

动画化这样的 map 可能是不可能的。

可能存在优化将像素转换为纬度的计算的空间，但这可能很困难。

例子

不幸的是，对于一段代码来说，代码太多了，因此我做了一些bl.ocks的演示。

Selectable tile set

Basic Zoom/Pan with vector layer

Restricted pan with fitSize

Static Map

这些仅进行了很少的测试，如果我设法完成基础的图块库，我将为此目的对其进行 fork ，与此同时，它应该足以作为示例。可以在d3-reprojectSlippy.js文件的geoTile.tile()中找到该代码的内容，该文件包含上述的输入/更新/退出循环(相当基本)和onload函数。当我在边上进行平铺工作时，我将不断更新此答案。

替代品

重新投影图块既麻烦又费时。如果可能的话，一种替代方法是将瓷砖设置为所需的投影。这是使用OSM tiles完成的，但是又麻烦又耗时-仅对于 map 制作者，而不是浏览器。

TL; DR

重新投影的墨卡托瓷砖需要时间，您应该阅读以上内容。