javascript - 使用 D3 将 Web Mercator 瓷砖重新投影到任意投影？

这个答案建立在：

这三种资源也相互依赖。理解这三个例子将有助于理解下面我的例子中发生的事情。

答案还以我缓慢、持续地尝试构建一个tile library为基础。

此答案的目标不是提供最终资源，而是粗略演示如何将资源与相关信息放在一起。随着我对这个问题的进一步思考，答案也会不断发展。

网页墨卡托瓷砖

延伸超过 360 度经度和约 170 度纬度（+/- 85 度）的墨卡托地图将填充一个正方形（纬度超过 85 度会导致失真失控，并且包含两极不是建议，因为极点在投影平面上处于 +/- 无穷大）。

对于网络地图服务（使用墨卡托瓦片），世界上大部分地区的这个正方形的缩放级别为 0。地图的宽度为 2^0 正方形，高度为 2^0 正方形。

如果我们将该正方形划分为两个正方形乘以两个正方形的网格，则缩放级别为 1。地图为 2^1 x 2^1 正方形。

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

对我们来说幸运的是，大约 85 度以北没有陆地，而且我们不经常想要展示南极洲，所以这个广场适合大多数网络地图应用程序。然而，这意味着如果我们将 web 墨卡托瓦片重新投影到显示在这些纬度之上的任何地方，我们将有一个光秃秃的地方：

wᴇʙ-mᴇʀᴄᴀᴛᴏʀsʀᴇᴘʀᴏᴊᴇᴄᴛᴇᴅғᴏʀᴀᴇʀᴄᴀᴛᴏʀᴘʀᴏᴊᴇᴄᴛɪᴏɴʜᴀʙᴇᴇɴʙᴇᴇɴʀᴏᴛᴀᴛᴇᴅʜᴏ

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

重投影和重采样板 Carree

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

投影函数是复杂的耗时操作，必须在渲染的每个像素上完成。我见过的所有 d3 示例都使用此处所示的过程（或近似变体）来进行实际的重投影和重采样。此示例仅适用于使用Plate Carree投影原始图像的情况。过程如下：

创建一个空白的新图像。
对于新图像中的每个像素，以像素为单位获取其位置并将其反转（使用所需的投影）以获得纬度和经度。
确定原始图像中的哪个像素与该经度和纬度重叠。
从原始图像中的该像素获取信息并将其分配给新图像中的适当像素（步骤 2 中的像素）

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

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

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

tɪᴍᴇᴛᴏᴘᴏɪɴᴛᴀᴍᴀᴘᴍᴀᴘᴍᴀᴘɪᴛʜɪᴛʜɪᴛʜɪᴛʜɪᴛʜɪᴛʜɪᴛʜɪᴛʜᴘɪᴘɪᴘɪᴘɪᴘɪᴘɪᴛᴏ/ʟᴏɴɢ/ʟᴏɴɢ）。我ᴛ ɪs ᴜɴᴄʟᴇᴀʀ ᴡʜʏ D3ᴠ4 ᴡᴏᴜʟᴅ ʙᴇ ғᴀsᴛᴇʀ。

为了完整起见，以下是 d3-geo-projections 中的更广泛的投影：

cᴏᴍᴘᴀʀɪɴɢᴛɪᴍᴇᴘʀᴏᴊᴇᴄᴛɪᴏɴɪɴᴠᴇʀᴛɪɴᴠᴇʀᴛᴀᴅᴅɪᴛɪᴏɴᴀʟ

笔记

这种方法可能会遇到这里描述的问题，这些问题可以通过那里的答案来解决 - 但不同的解决方案将花费额外的时间来处理。
这种方法使用最近邻方法——这可能会导致质量问题。更高级的采样（例如双线性或三次采样）会增加处理时间，但可能会产生更理想的图像。
如果基本图像具有文本，则可以旋转或以其他方式操纵文本以使其不那么可读或不可读。
Mike 的示例是针对单个图像，对于瓦片，该过程在某种程度上有所改变，我们现在正在创建多个图像，这需要知道每个原始瓦片的边界和每个重新投影瓦片的度数以及前者的瓦片单位和后者的像素 - 次要细节。

重投影和重采样 Web Mercator

当我开始查看这个问题时，我查看了 Alan McConchie 的解决方案/示例作为参考。花了一些时间才注意到，但本示例中的第 3 步（我相信 Jason Davies 的工作也是如此）并未考虑重采样中的 Web 墨卡托图块 - 仅用于确定图块边界。但是，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 坐标。

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

寻找合适的瓷砖

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

对于 McConchie/Davies，d3.quadTile 将在给定具有裁剪范围（不是裁剪角度）和切片深度的投影的情况下，拉出与视图范围相交的所有切片。

在 Alan McConchie 的解决方案/示例中，缩放级别基于投影比例 - 但这不一定是最明智的：每个投影都有不同的比例因子，一个比例上的 100 比例将显示与 a 比例不同的程度100 另一个。此外，圆柱投影中比例值和地图大小之间的关系可能是线性的，而非圆柱投影可能在地图大小和比例之间具有非线性关系。

我稍微修改了这种方法 - 我使用比例因子来确定初始瓦片深度，然后如果 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;
}

起初，我认为从多个缩放深度拉瓷砖（以考虑重新投影瓷砖的大小差异）是理想的：但这会遇到光栅中的线条粗细以及不连续注释等问题。

采用 Jason 和 Alan 的工作

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

重采样不考虑网络墨卡托（如上所述）
瓷砖深度没有选好（上面提到过）
瓷砖被重新投影为放置在 div 而不是 SVG 中的 canvas-es - 创建两个父容器，一个用于每种类型的功能：瓷砖或矢量。

我相信我的例子，经过非常细微的调整，已经解决了这些问题。我还添加了一些更广泛的评论以供查看：

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

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

具有重叠特征的常规瓦片地图有几个坐标系：

平铺单元 (3D)，标记每个平铺的列、行和缩放级别（分别为 x、y、z）
地理坐标 (3D)，用于标记三维球体上某个点的纬度和经度。
缩放单位 (3D)，跟踪缩放平移 (x,y) 和缩放比例 (k)。
投影单位 (2D)，将经纬度投影到的像素单位。

任何滑动地图的目标都是在可用系统中统一这些坐标。

当我们重新投影瓦片时，我们需要添加一个坐标空间：

(/a) 瓦片集投影。

我觉得这些示例在如何将所有坐标系联系在一起方面并不是特别清楚。因此，正如您可能已经看到的，我将上述方法放置在一个 geoTile 对象中，该对象取自个人项目的tile 库。这样做的目的是为了更顺畅地协调不同的单位。我不是想插上它，无论如何它仍在开发中（只是太忙而无法真正完成它）；但是，我会看看时间是否让我有机会用d3-tile装配一个示例。

前进的挑战

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

动画这样的地图可能是不可能的。

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

例子

不幸的是，代码对于一个片段来说太多了，所以我做了几个块来演示。

这些只进行了最少的测试，如果我设法完成了基础磁贴库，我将为此目的分叉它，同时它应该足以作为示例。代码的geoTile.tile()核心位于 d3-reprojectSlippy.js 文件中，该文件包含进入/更新/退出循环（相当基本）和上面描述的 onload 函数。当我在瓷砖上工作时，我会保持更新这个答案。

替代方案

重新投影瓷砖既麻烦又耗时。如果可能的话，另一种选择是在所需的投影中生成一个瓦片集。这已通过OSM 瓦片完成，但也很麻烦且耗时 - 仅适用于地图制作者，而不是浏览器。

TL;博士

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