使用 d3-geo 在球体的一部分上绘制具有公制尺寸的多多边形

Posted 2023-03-23

【中文标题】使用 d3-geo 在球体的一部分上绘制具有公制尺寸的多多边形

【英文标题】：Draw a multi polygon with metric dimensions over part of a sphere, using d3-geo

【发布时间】：2018-06-26 16:52:34

【问题描述】：

我有以下需求：

给定一个具有公制尺寸的多面体，我需要将其绘制在球体的特定部分上，然后以 GeoJSON 格式输出。

此外，“球体的一部分”是一个由 4 个点组成的矩形，可以相对于赤道线旋转。

举个例子，想象一个 MultiPolygon，里面有一个矩形，像这样，其中的数字以米为单位：

[
  [
    [
      [0, 10],
      [0, 0],
      [20, 0],
      [20, 10]
    ]
  ]
]

如您所见，这样的多边形与轴成直角对齐。

然后，给定球体的任何矩形区域，例如下面的那个：

    leftTop: [ 13.377948, 52.521293 ],
    leftBottom: [ 13.377962, 52.521250 ],
    rightBottom: [ 13.378174, 52.521275 ],
    rightTop: [ 13.378161, 52.521318 ]

我想在这个区域平移、旋转和适应那个多面体。

因此，d3-geo 包似乎提供了适当的工具，例如 projection 结合了 scale/translate/rotate/fitExtent 和 invert 等。但我找不到合适的组合来处理旋转的矩形区域。

由于网络中 d3-geo 的绝大多数示例要么从 GeoJSON 到像素，要么使用 d3 作为 SVG 渲染引擎等。我很难将我的需求抽象为函数调用。

有没有人可以帮我解决这个问题？我猜几何/三角学中一定有一个合适的术语来描述这种计算。

我自己用过“在测地线区域上投影平面多边形”，但我想一定有一个更好、更具体的术语。至少如果我找到这样的术语，那将是一个好的开始。

更新：

我制作了一些插图来说明我正在努力实现的目标。请记住，尽管它们在这里是可视图像，但在我的情况下，我并不是在寻找可视化渲染工具，因为我只会将 Array/Objects 转换为 GeoJSON。

我要投影的多边形。请注意，浅灰色框显示了该多边形所在的虚拟区域，因此，由于多边形位于中间，它也应该投影到地理区域的中间。

大地水准面（本例中为地球）上的矩形扇区。注意指向参考角的红色标签，这应该表示旋转。

结果，多边形投影到给定区域（实际上是 GeoJSON 格式，具有相应的坐标）。

【问题讨论】：

如果我理解正确，第一个多边形（以米为单位）仅用于比例 - 您想拉伸/挤压和旋转它以适应某些投影特征吗？就像以像素为单位拉伸具有宽度/高度的卫星图像，以便将照片中的 x、y 位置转换为球体上的坐标或某些投影的地理坐标空间？

是的，正确的。我遇到的主要问题是第二个区域（投影特征）可以自行旋转，而不是与赤道线对齐。因此，假设一个适合该特征的“星形”多边形最终会一起旋转。因此，仅使用视口投影 x,y 是不够的。它需要投影在这 4 个角内。这就是我有限的数学打击我的地方：P

以正方形为例，将其与投影为正方形的地理特征进行匹配，应该如何确定旋转？可以使用四种不同的旋转，这有关系吗？是否有任何与非地理形状无关的球体旋转迹象？

旋转由第二个矩形特征决定。如果你仔细观察，就像这里一样 - bl.ocks.org/d/ac0cfaabcd5fcbd11814439203fe68b4 - 这不是一个正方形，它旋转了大约 10 到 15 度，但 4 个角暗示了旋转。在这种情况下，[0,10] 将被投影为 [ 13.377948, 52.521293 ] 和 [0, 0] 将是 [ 13.377962, 52.521250 ]（请注意，X、Y 都移动了，而在第一个方格中，X 始终为 0）。

另外，由于该特征是一个带有“leftTop”和“leftBottom”等键的对象，这也表示旋转，因为[0,0]应该投影到“leftBottom”中的任何坐标，而[0,10] 将被投影为与“leftTop”中的相同。因此，如果特征旋转 180 度，“leftBottom”将具有更高的“leftTop”，即使正方形仍然相同。这说明清楚了吗？

【参考方案1】：

我找到了一种实现我需要的方法，虽然我没有使用 d3-geo，而是使用了 Turf.js 和 的组合geodesy 图书馆。

在这里（搜索“projectPlainOverGeoArea”）：

https://gist.github.com/marinho/8a7e71b53c7da97f9579146742bb54de

代码执行以下操作：

计算距地理矩形（我要投影到的区域）以米为单位的宽度和高度。它使用大地测量球面/正弦函数

使用从多边形矩形到地理矩形的比例距离将多边形的每个路径的每个点投影到等效的地理坐标上。此处的结果是预期测量值与赤道线对齐的多边形。

使用 Turf 的 transformRotate 以左/下为轴心旋转投影多边形

这似乎在我所做的测试中运行良好，但老实说我对此感到不满意，因为我觉得必须有一个更优雅的解决方案 d3-geo 或另一个图书馆。

【讨论】：

我不确定是否有很多优雅的方法可以做到这一点，我已经添加了我对使用 d3-geo 执行此操作的想法，但我认为匹配球形“矩形”总是不优雅平面矩形总是会导致一些不愉快。

【参考方案2】：

如果我正确理解问题，我相信我有解决方案。但是，这个问题很宽泛，没有任何代码，因此，我的回答会很宽泛，但应该清楚地概述实现预期结果的方法。

根据图像，问题似乎是：如何将形状适合地理边界框，其中边界框投影到平面空间并可能具有旋转。

如果这是准确的，那么我的回答应该很有用。为了清楚起见，当我使用形状一词时，我指的是非地理多边形。

解决的一般模式是：

获取参考多边形的方向 旋转地图 - 不是形状 使用 geoIdentity，使形状适合地理边界框

一个。完成与开始时相反的旋转，

b.旋转整个 svg/canvas 以对抗初始旋转；或者，

c。保留旋转后的地图。

---

获取地理矩形的方向

您需要建立地理特征相对于投影平面的旋转，这似乎是问题的症结（您说明了球体的旋转矩形区域，但区域只能是投影中的矩形空间，并且根据图像，这种解释似乎是正确的）。

D3 在这方面可能有点挑战，因为它不会沿笛卡尔坐标渲染路径，而是插入很大的圆距离。对于国家或大陆规模的地理特征，这将是一个更大的问题，但在城市规模上应该可以忽略不计。

如果从左上角 ([x1,y1]) 开始并使用右下角 ([x2,x1])，您应该能够确定这条线相对于它应该是什么的角度：一条垂直线从第一个坐标开始，到第二个坐标结束。鉴于缠绕顺序对 d3 很重要，如果您有一个点，则第二个坐标将始终对应于矩形的同一顶点。

获得这个角度的方法相当简单：

var p1 = [long,lat];   // geographic coordinate space for the two points
var p2 = [long,lat];

var x1 = projection(p1)[0];  // projected svg coordinate space for the two points
var x2 = projection(p2)[0];
var y1 = projection(p1)[1];
var y2 = projection(p2)[1];

var dx = x1 - x2;  // run
var dy = y1 - y2;  // rise

var angle = Math.atan(dx/dy);  // in radians, multiply by 180/π to get degrees

我们正在投影两个坐标，计算 x 和 y 坐标的投影差异，以像素为单位。通过 Math.atan() 运行它，我们就有了一个角度。

但是等等，还有更多。

旋转地图

我们用来计算角度的方法很好，但是我们需要修改它来旋转地图。如果我们旋转地图，它将围绕[0,0] (long,lat) 旋转，这是大多数投影的默认旋转中心。我们需要将地图居中于[-x,-y]，其中 x 和 y 表示用于计算方向的两个点的中点，在地理坐标空间中测量（未投影）。

在计算角度之前，我们需要通过旋转使地图居中，因为这会改变角度。为此，我们需要 d3.geoInterpolate 来计算从 p1 (0) 到 p2 (1) 的点，我们需要中间点，所以我们输入 0.5：

var pMid = d3.geoInterpolate(p1, p2)(0.5);

现在我们可以在计算角度之前将其应用于投影：

projection.rotate([-pMid[0],-pMid[1]]);

为什么是负值？我们在我们脚下移动地球

现在我们可以进行角度计算了。一旦我们有了角度，我们就可以应用它：

projection.rotate([-pMid[0],-pMid[1],-angle])  // angle in degrees

我们已经完成了一半，使用下图，我们使用 A 和 B 的地理坐标来确定地理中心 C。然后使用 A 和 B 的投影坐标确定角度 α，然后使用它来旋转投影坐标，使 AB 线在地图上垂直。

使形状适合地理边界框

所以我们解决了一半的问题，现在我们需要投影形状。我们将为形状使用第二个投影，一个普通的 geoIdentity 就可以了，它允许我们在投影坐标时使用 fitSize 或 fitExtent 方法而不进行变换。请注意，您可能希望在 y 轴上翻转此功能：svg y 值从顶部 0 开始，更标准的笛卡尔 y 值从底部开始。

我们将要使用 fitExtent，它允许我们为形状设置一个边界矩形。 proejction.fitExtent([[x,y],[x,y]],feature) 采用一个包含边界框左上角和右下角（在 svg 坐标中）的数组来保存特征（geojson 特征）。

请记住，我们纠正的投影地理特征的地理（未投影）尺寸越小，它的矩形就越长，更大的区域可能具有更少的正方形属性，尤其是相对于矩形的右侧.对于较大的地理区域，您可能需要修改旋转计算，但有时投影的矩形与球体上的“矩形”不对齐。

要获得 fitExtent 的边界框，我们可以使用 path.bounds()：

返回投影的平面边界框（通常以像素为单位） 指定的 GeoJSON 对象。边界框由 二维数组：[[x₀, y₀], [x₁, y₁]]，其中 x₀ 是最小值 x坐标，y₀是最小值 y坐标，x₁是最大值 x 坐标，y₁ 是最大 y 坐标。 (API docs)

太好了，现在我们有两个应该具有相同点的边界框，我们使用：

projection2.fitExtent(path.bounds(geoRectangle));

现在我们已经在地理特征上覆盖了形状：

var feature =  "type": "Feature","properties": ,"geometry": "type": "Polygon","coordinates": [ [ [-81.62841796875,24.307053283225915],[-75.9375,21.88188980762927],[ -77.3876953125,18.8543103618898],[       -83.1884765625,21.268899719967695 ], [-81.62841796875,24.307053283225915]]];
var triangle = "type": "Polygon","coordinates": [[[0, 0], [10, 10], [20, 0], [0, 0]]];


var width = 500; var height = 300;
var svg = d3.select("body")
  .append("svg")
  .attr("width",width)
  .attr("height",height);
  
var projection = d3.geoMercator().scale(500/Math.PI/2).translate([width/2,height/2]);
var path = d3.geoPath().projection(projection);


d3.json("https://unpkg.com/world-atlas@1/world/110m.json", function(error, world) 
  if (error) throw error;

  var p1 = feature.geometry.coordinates[0][1]; // first point in geojson
  var p2 = feature.geometry.coordinates[0][2]; // second point in geojson
  
  var pMid = d3.geoInterpolate(p1, p2)(0.5);   // halfway point between both points
  
  projection.rotate([-pMid[0],-pMid[1]]);                          // rotate the projection to center on the mid point
  projection.fitExtent([[135,135],[width-135,height-135]],feature) // optional: scale the projected feature, may offer benefits for very small features
 
  var dx = projection(p1)[0] - projection(p2)[0]; // run, difference between projected points x values
  var dy = projection(p1)[1] - projection(p2)[1]; // rise, difference between projected points y values
  
  var a = Math.atan(dx/dy) * 180 / Math.PI; // get angle and convert to degrees
  
   projection.rotate([-pMid[0],-pMid[1],-a]);  // adjust rotation to straighten feature
   projection.fitExtent([[135,135],[width-135,height-135]],feature) // scale and translate the feature.
   
  // draw world map, draw feature
  svg.append("path")
    .attr("d",path(topojson.mesh(world)))
    .attr("fill","none")
    .attr("stroke","black")
  
  svg.append("path")
    .attr("d", path(feature))
    .attr("fill","none")
    .attr("stroke","steelblue");


  // set up the projection and path for the shape       
  var projection2 = d3.geoIdentity().reflectY(true);
  var path2 = d3.geoPath().projection(projection2);

  // scale the shape's projection for the shape using the bounds of the geographic feature
  projection2.fitExtent(path.bounds(feature),triangle);

  // draw the shape
  svg.append("path")
    .attr("d", path2(triangle));
  
);

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.11.0/d3.min.js"></script>
<script src="https://unpkg.com/topojson-client@3"></script>

这是一个手绘的多边形，我在显示器上放了一张名片，在上面描了画，所以瑕疵比一般的大

结束游戏

现在您已经绘制了一个形状并将其覆盖在某个地理特征之上，希望相当体面。怎么办？从技术上讲，这解决了关键问题，但如果您不想要倾斜的地图，我们已经创建了一个新问题。解决此问题的选项是旋转整个 svg/canvas 或获取形状中的每个点并重新投影它们并将它们转换为地理坐标，而不是无聊的旧笛卡尔坐标。还有其他的，但这两个似乎最直接。

如果处理图像，您可以逐个像素地重新投影图像，但不要指望这会很快，请参阅this block。这个answer 会查看图像并将它们拟合到投影（如果它们具有已知的边界），墨卡托应该在非常小的距离内作为投影正常工作。

重投影点和 GeoIdentity

如 cmets 所示，如果遵循以下模式，d3.geoIdentity 将无法进行重投影：var latlong = projection.invert(projection2([x,y]))，因为 geoIdentity 不返回函数，而只是返回一个对象。相反，我们可以定义行为（翻转 y，并使用恒等比例和变换）：

projection2.project = function([x,y]) 
  var s = this.scale();
  var t = this.translate();
  return [(x * s) + t[0] , -y * s + t[1]]

在粗略测试中似乎返回投影点，这是 geoIdentity 中不包含的 geoProjection 的标准行为。使用模式应如下所示：

projection.invert(projection2.project([x,y]))

【讨论】：

谢谢。我理解了这个概念，但最后一行 projection.invert(projection2([x,y]) 失败了，因为它抱怨 projection2 不是一个函数。文档在github.com/d3/d3-geo/blob/master/README.md#_projection 中说它应该是可调用的，并且在 3.0 版本中它也曾经是可调用的，所以，我仍在努力解决这个问题。

这令人惊讶，因为此功能看起来很有用、易于实现且符合预期。不过，我们可以自己定义功能，我已经更新了答案以反映这一点。