将Web Mercator磁贴重新投影到D3的任意投影？
Ahm*_*sih 9 javascript gis map-projections d3.js
杰森·戴维斯（Jason Davies）用“ 重新投影的栅格图块”引爆我们已经过去了几年，该地图因Mapbox阻止了他的站点而停止工作，但是Mollweide水彩画和《中断的古德栅格》仍然是不错的演示。
现在在Observable HQ上，我看到了有关最新d3-geo-projection和d3-tile的文档，但是没有关于Jason所做的操作的现代示例：重新投影标准Mercator贴图集。
我如何才能使d3平铺扭曲到新的投影？
这个答案建立在：
这三种资源也是相互建立的。理解这三个示例将有助于理解我下面的示例中发生的事情。
答案还使用我缓慢移动、持续尝试构建tile library作为基础。
这个答案的目标不是提供最终的资源，而是粗略地展示如何将一个资源与相关信息放在一起。随着我对这个问题的进一步思考，答案也会不断发展。
网络墨卡托瓷砖
延伸超过 360 度经度和约 170 度纬度（+/- 85 度）的墨卡托地图将填充一个正方形（超过 85 度纬度会导致失真失控，并且不会包含极点）建议，因为极点在投影平面上的 +/- 无穷大）。
对于 Web 地图服务（使用墨卡托图块），这个世界大部分地区的缩放级别为 0。地图的宽度为 2^0 平方，高度为 2^0 平方。
如果我们将该正方形划分为两个正方形乘以两个正方形的网格，则缩放级别为 1。地图为 2^1 x 2^1 个正方形。
因此，缩放级别决定了地图的宽度和高度：2^zoomLevel。如果每个正方形的像素大小相同，那么缩放级别每增加 1 就会使世界的像素宽度增加 2 倍。
对我们来说幸运的是，北面没有约 85 度的陆地，而且我们并不经常想要展示南极洲，所以这个方块适合大多数网络地图应用程序。但是，这意味着如果我们将 Web 墨卡托图块重新投影到显示在这些纬度以上的任何内容，我们将有一个秃斑：
W??-M??????? ????s ????????????? ？？？? 米？？？？？？？？？？？？？？？？？？??s ???? ??????? ?? s？？？？？？？？？？P？？？。
最后，Web 墨卡托图块在相对于图块可预测且规则的投影空间中呈现。如果我们重新投影瓷砖，我们可能会扩大或缩小每个瓷砖的投影空间，我们应该注意这一点。在上图中，北极周围的瓦片被重新投影的比更靠南的瓦片小得多。投影后瓷砖的大小不一定一致。
重投影和重采样板卡里
重新投影 Web 服务磁贴的最大挑战是时间——而不仅仅是花时间理解投影和阅读这样的答案。
投影函数是复杂的耗时操作，必须在渲染的每个像素上完成。我见过的所有 d3 示例都使用此处所见的过程（或近似变体）进行实际的重投影和重采样。此示例仅适用于使用Plate Carree投影原始图像的情况。过程如下：
- 创建一个空白的新图像。
- 对于新图像中的每个像素，以像素为单位获取其位置并将其反转（使用所需的投影）以获得纬度和经度。
- 确定原始图像中的哪个像素与该经度和纬度重叠。
- 从原始图像中的那个像素中获取信息，并将其分配给新图像中的适当像素（步骤 2 中的像素）
当原始图像使用 Plate Carree 投影时，我们不需要 d3-geoProjection，投影和未投影坐标之间的关系是线性的。例如：如果图像高 180 像素，则每个像素代表 1 度纬度。这意味着与步骤 2 和 project.invert() 相比，步骤 3 不需要很长时间。这是 Mike 在第 3 步中的函数：
var q = ((90 - ?) / 180 * dy | 0) * dx + ((180 + ?) / 360 * dx | 0) << 2;
步骤 2 所需的时间与用于重投影图像的投影有关。我见过的所有示例都d3.geoProjection.invert()用于上面列表中的第二步 - 在新图像中获取像素位置并找出其纬度和经度。并非所有的预测都是生而平等的。圆柱投影通常优于圆锥投影，圆锥投影通常优于方位投影。我还看到 d3v4 和 d3v5 之间在投影.invert() 时间方面存在一些奇怪的速度差异：
T？？？?? ????????? ? ????? ?? ? ？？？？？？D3 (??????? ???? ??x??s ?? ???/????)。一世？?s ??????? ？？？D3？4？？？？?? ??s???。
为了完整起见，以下是 d3-geo-projections 中更广泛的投影：
C？？？？？？？？？？？？？？????????????.?????? ？？？？？？？？？？？？??????????s
笔记
这种方法可能会遇到这里描述的问题，这些问题可以通过那里的答案解决 - 但不同的解决方案将花费额外的时间来处理。
此方法使用最近邻方法 - 这可能会导致质量问题。更高级的采样（例如双线性或三次采样）会增加该过程的时间，但可能会产生更理想的图像。
如果基本图像有文本，则文本可能会被旋转或以其他方式处理以使其不那么可读或不可读。
Mike 的示例是针对单个图像，对于图块，该过程在某种程度上有所改变，我们现在正在创建多个图像，这需要知道每个原始图块的边界和每个重新投影图块的边界（以度为单位）以及前者的图块单位和后者的像素 - 次要细节。
重投影和重采样 Web 墨卡托
当我开始研究这个问题时，我将 Alan McConchie 的解决方案/示例作为参考。花了一段时间才注意到，但本示例中的第 3 步（我也相信 Jason Davies 的工作）并未考虑重采样中的 Web 墨卡托图块 - 仅用于确定图块边界。但是，y 轴上的像素之间的关系不再像 Plate Carree 中那样是线性的。
这意味着瓦片放置在正确的位置，但采样将 y 轴视为每个瓦片内的线性。这种失真在以较低的图块缩放级别（图块的中下方/上方）显示整个世界时最为明显，这可能就是 Alan 在提到奇数压缩时所说的。
解决方案是在上面的步骤 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 比例将显示不同的程度。 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 基本没有变化，我已经解决了我在这段代码中看到的以下挑战：
- 重采样没有考虑到 web 墨卡托（上面提到的）
- 瓷砖深度没有选择好（上面提到的）
- 瓷砖被重新投影为放置在 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 library的个人项目。这样做的目的是让不同单位之间的协调更加顺畅。我不是要插入它，无论如何它仍在开发中（只是太忙而无法真正完成它）；但是，我会看看时间是否能让我有机会用d3-tile装配一个示例。
前进的挑战
变焦速度和响应能力是我看到的最大挑战。为了解决这个问题，我将缩放功能设置为在缩放结束时触发 - 这在平移事件中最为明显，因为通常平移通过平移连续触发缩放功能，这可以通过转换现有图像来解决。然而，最可靠的使用方法是在静态地图上。为已经绘制的图像实现平移将是平移事件的理想选择，而不是像目前那样重新采样。
动画这样的地图可能是不可能的。
可能有优化将像素转换为纬度的计算的空间，但这可能很困难。
例子
不幸的是，代码对于一个片段来说太多了，所以我做了一些 bl.ocks 来演示。
这些只进行了最少的测试，如果我设法完成基础 tile 库，我会为此目的将其分叉，同时它应该足以作为示例。代码的内容可以在geoTile.tile()d3-reprojectSlippy.js 文件中找到，该文件包含进入/更新/退出循环（相当基本）和上述的 onload 函数。当我在瓷砖上工作时，我会不断更新这个答案。
替代方案
重新投影瓷砖既麻烦又费时。如果可能，另一种方法是在所需的投影中生成图块集。这已通过OSM 瓷砖完成，但也很麻烦且耗时 - 仅适用于地图制作者，而不是浏览器。
TL; 博士
重新投影的墨卡托瓷砖需要时间，您应该阅读以上内容。