使用 D3 将 Web Mercator 磁贴重新投影到任意投影?

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

杰森·戴维斯 (Jason Davies) 以 Reprojected Raster Tiles—that map stopped working because Mapbox is blocking his site, but the Mollweide Watercolour and Interrupted Goode Raster 的出色演示让我们震惊已经过去几年了。

现在在 Observable HQ 上,我看到了最新 d3-geo-projection and d3-tile 的文档,但没有关于如何执行 Jason 所做的现代示例:重新投影标准墨卡托图块集。

如何让 d3-tile 扭曲到新的投影?

此答案基于:

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

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

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

Web 墨卡托拼贴

延伸超过 360 度经度和约 170 度纬度(+/- 85 度)的墨卡托地图将填满一个正方形(超过 85 度纬度会导致变形失控,并且包含极点是不可取的,因为极点在投影平面上处于 +/- 无穷远)。

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

如果我们 divide 将那个正方形分成两个正方形乘两个正方形的网格,我们的缩放级别为 1。地图是 2^1 x 2^1 正方形。

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

对我们来说幸运的是,~85 度以北没有陆地,而且我们通常不想显示南极洲,因此这个正方形适合大多数网络地图应用程序。然而,这意味着如果我们将 web 墨卡托图块重新投影到显示在这些纬度以上的任何地方,我们就会有一个秃点:

wᴇʙ-mᴇʀᴄᴀᴛᴏʀsʀᴇᴘʀᴏᴊᴇᴄᴛᴇᴅғᴏʀᴀᴀᴘʀᴏᴊᴇᴄᴛɪᴏɴʜᴀʜᴀSʙᴇᴇɴSʜᴏSʜᴏSʜᴏSʜᴏSʜᴏSʜᴏnᴏʀᴛʜpᴏʟᴇ。

最后,Web 墨卡托图块在投影 space 中呈现,相对于图块是可预测和规则的。如果我们 re-project 瓷砖,我们可能会放大或缩小每个瓷砖的投影 space,我们应该注意这一点。在上图中,北极周围的地块经过重新投影后比更南端的要小得多。 Tiles 投影后大小不一定统一

再投影和重采样板 Carree

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

投影函数是复杂的耗时操作,必须在渲染的每个像素上完成。我见过的所有 d3 示例都使用 here for the actual reprojection and resampling. This example will only work if the original image is projected with Plate Carree 中所见的过程(或接近的变体)。过程如下:

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

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

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

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

tɪᴍᴇᴛᴏᴀᴀᴍᴀᴘ Iᴛɪs ᴜɴᴄʟᴇᴀʀ ᴡʜʏ D3ᴠ4 ᴡᴏᴜʟᴅ ʙᴇ ͓ᴀsᴛᴇʀ.

为了完整起见,这里是 d3-geo-projections:

中的 wider 投影范围

cᴏᴍᴘᴀʀɪɴɢᴏᴍᴘᴀʀɪɴɢᴘʀᴏᴊᴇᴄᴛɪᴏɴᴘʀᴏᴊᴇᴄᴛɪᴏɴɪɴᴠᴇʀᴛғᴏʀS

备注

  • 这种方法可以运行转化为描述的问题,这些问题可以通过那里的答案来解决——但不同的解决方案将花费额外的时间来处理。

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

  • 如果基本图像有文本,文本可能会被旋转或以其他方式处理以使其变少或不可读。

  • Mike 的示例是针对单个图像的,对于图块,过程在一定程度上有所改变,我们现在正在创建多个图像,这需要知道每个原始图块的边界和每个重新投影的图块的边界度数以及前者的平铺单位和后者的像素 - 次要细节。

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

当我开始看这个问题时,我参考了 Alan McConchie 的 solution/example。花了一段时间才注意到,但此示例中的第 3 步(我相信 Jason Davies 的工作也是如此)在重采样中没有考虑 web 墨卡托图块 - 仅在确定图块边界时。但是,y 轴上的像素之间的关系不再像 Plate Carree 中那样是线性的。

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

解决方案是在上面的第 3 步中正确投影每个 latitude/longitude 对的纬度。这增加了时间,总是更多的时间——函数涉及 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 我只看到一种干净的方法来查找要显示的图块,我看到 here. I believe that Alan McConchie uses an unminified version that may be otherwise altered. There is also this github 另一个版本的 d3.quadTiles 的存储库,这非常棒相似的。

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

在 Alan McConchie 的 solution/example 中,缩放级别基于投影比例 - 但这不一定是最明智的:每个投影都有不同的比例因子,一个比例为 100 将显示一个与另一个 100 的等级不同。此外,在圆柱投影中比例值和地图大小之间的关系可能是线性的,而 non-cylindrical 投影在地图大小和比例之间可能具有非线性关系。

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

起初我认为从多个缩放深度拉出图块(以解决重新投影图块大小的差异)会 ideal:但这会 运行光栅中的线条粗细以及不连续的注释等问题。

采用 Jason 和 Alan 的作品

我将上面使用 geoTile.tiles() 和 运行 生成的图块集通过使用图块坐标(在图块坐标、行、列、缩放深度中)的 enter/update/exit 循环作为一个键,将 image 元素附加到父 gsvg。加载图像时,一旦图像加载完毕,我们就会调用 onload 函数来进行实际的重新投影。这与 Jason 和 Alan 在很大程度上没有变化,我已经解决了我在这段代码中看到的以下挑战:

  • 重采样未考虑 web 墨卡托(如上所述)
  • 图块深度没有选择好(上面提到)
  • 瓷砖被重新投影为 canvas-es 放置在 div 而不是 SVG - 创建两个父容器,一个用于每种类型的特征:瓷砖或矢量。

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

    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),标记纬度和经度f 三维球体上的一个点。
  • 缩放单位 (3D),跟踪缩放平移 (x,y) 和缩放比例 (k)。
  • 投影单位(2D),纬度和经度投影到的像素单位。

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

当我们重新投影瓦片时,我们需要添加一个坐标space:

  • (/a) 瓦片集投影。

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

来组装一个例子

前进的挑战

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

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

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

例子

不幸的是,代码对于一个片段来说太多了,所以我做了一些 bl.ocks 来演示。

这些只进行了最少的测试,如果我设法完成基础图块库,我会为此目的分叉它,同时它应该足以作为一个例子。代码的核心位于 d3-reprojectSlippy.js 文件的 geoTile.tile() 中,其中包含 enter/update/exit 循环(相当基本)和上述的 onload 函数。当我在 side 上处理瓷砖时,我会及时更新这个答案。

备选方案

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

TL;DR

重新投影的墨卡托图块需要时间,您应该阅读以上内容。