地图上的多边形交点显示偏移

Question

我正在尝试查找地图上绘制的 2 个多边形之间的交叉区域。 我有 TurfJS intersect 方法来查找 2 个多边形之间的交集。

它适用于较小的区域，但对于较大区域的多边形，它开始在交点处显示一些偏移量，偏移量随着距离的增加而增加。 这也仅适用于带有斜线的多边形（垂直和水平线的交叉点似乎效果很好）。

我为此创建了一个 JSFiddle：https://jsfiddle.net/cLe6yo9d/

我试图找到黑色和蓝色多边形之间的交点，我得到的显示为红色多边形，它看起来从它应该在的位置移动。

var mapLayer = L.map('mapid', {
 zoomAnimation: false
});
var pid = 'karan44.pdmio34k';
var at = 'pk.mapbox-access-token-goes-here';
L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token={accessToken}', {
  id: pid,
  accessToken: at
}).addTo(mapLayer);

var polygon1 = turf.polygon([
  [
    [3.405762, 51.395350],
    [5.009766, 53.340303],
    [7.141113, 53.653999],
    [5.822754, 51.037508],
    [3.405762, 51.395350]
  ]
], {
  "fill": "#00000F",
  "stroke": "#00000F",
  "stroke-width": 1
});

var polygon2 = turf.polygon([
  [
    [0.241699, 54.173488],
    [10.162354, 50.908012],
    [8.854980, 50.062208],
    [0.241699, 54.173488]
  ]
], {
  "fill": "#0000FF",
  "stroke": "#0000FF",
  "stroke-width": 1
});

var polygon = turf.intersect(polygon1, polygon2);
polygon.properties = {
  "fill": "#FF0000",
  "stroke": "#FF0000",
  "stroke-width": 1
};


L.mapbox.featureLayer().setGeoJSON(polygon1).addTo(mapLayer);
L.mapbox.featureLayer().setGeoJSON(polygon2).addTo(mapLayer);

L.mapbox.featureLayer().setGeoJSON(polygon).addTo(mapLayer);
mapLayer.setView([52.754260888947776, 5.72100021667583], 8);

这个小提琴是通过修改turf.intersect example来重现问题而创建的。

希望有人可以帮助我了解问题所在。

Answer 1

我不清楚 turf 在计算交点时在做什么（它们不像 Ivan 所怀疑的那样位于多边形顶点之间的大圆上），但可以通过沿大圆弧插值多边形来解决问题，然后在这些上使用turf.intersect。下面的函数将采用 GeoJSON 多边形（如使用 turf.polygon 生成的多边形）并输出新的多边形，每个边都插入到测地线弧中（可选地指定要使用的步骤数）：

//interpolates simple GeoJSON polygon features along geodesic arcs
function geodesify(input, steps) {
  if (typeof steps === 'undefined') {
    steps = 50; //interpolation steps on each segment
  }
  var tempLine = {
    "type": "Feature",
    "geometry": {
      "type": "LineString",
      "coordinates": []
    },
    "properties": {}
  }
  if (input.geometry.type === "Polygon") {
    tempLine.geometry.coordinates = input.geometry.coordinates[0];
    tempLine.properties = input.properties;
    tempLine.properties.geodesic = "true"; //tells Leaflet.Geodesic to interpolate this feature
    tempLine.properties.geodesic_steps = steps;
    var geoLine = L.geoJson(tempLine).toGeoJSON();//convert interpolated feature back to GeoJSON
    var output = {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": geoLine.features[0].geometry.coordinates
      },
      "properties": tempLine.properties
    };
    output.properties.geodesic = "false"; //to prevent a second interpolation
    var outLen = output.geometry.coordinates[0].length;
    output.geometry.coordinates[0][outLen-1] = output.geometry.coordinates[0][0];
    return output;
  }
  console.log("geodesify input geometry must be a GeoJSON Polygon Feature");
  return false;
}

它依赖于Leaflet.Geodesic plugin，它将插入具有名为geodesic 的属性设置为"true" 的GeoJSON LineString 特征。大部分代码都用于将 GeoJSON 从 Polygon 转换为 LineString 并返回到 Polygon，如果插件一开始就接受 Polygon 功能，这将是不必要的。但无论如何，在你的例子中，你会像这样使用它：

var polygon1a = geodesify(polygon1, 30);
var polygon2a = geodesify(polygon2, 30);

var polygon = turf.intersect(polygon1a, polygon2a);
polygon.properties = {
  "fill": "#FF0000",
  "stroke": "#FF0000",
  "stroke-width": 1
};

L.mapbox.featureLayer().setGeoJSON(polygon1a).addTo(mapLayer);
L.mapbox.featureLayer().setGeoJSON(polygon2a).addTo(mapLayer);

L.mapbox.featureLayer().setGeoJSON(polygon).addTo(mapLayer);

这是一个在工作中展示它的小提琴：

https://jsfiddle.net/nathansnider/ycnno5df/

请务必注意，此功能仅在输入为单个 GeoJSON Polygon 功能时才有效，尽管它可以适用于其他类型。 Leaflet.Geodesic 的使用在这里可能有点过分，但它确实有效。我怀疑任何足够精细的插值实际上都可以正常工作（也就是说，turf 会计算在地图上绘制时沿着多边形边缘落下的交点），但这种方法具有地理上正确的优势。

Answer 2

这个问题又来了，终于找到了解决办法。分享以防万一它可能对其他人有帮助。

我对这个问题的理解是 turf.js 和 LeafletJs 使用不同的投影系统。所以我在传递给 turfjs 方法之前手动将坐标投影到 EPSG:3857，然后在传递给 Leaflet 之前将 turfjs 结果转换回 EPSG:4326。

https://jsfiddle.net/5Ls3vw8d/

var mapLayer = L.map('mapid', {
    zoomAnimation: false
});
var pid = 'karan44.pdmio34k';
var at = 'pk.mapbox-access-token-goes-here';
L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token={accessToken}', {
  id: pid,
  accessToken: at
}).addTo(mapLayer);
//Construct Truf projection polygons
var polygon1 = turf.polygon([[
    proj4("EPSG:4326","EPSG:3857",[3.405762, 51.395350]),
    proj4("EPSG:4326","EPSG:3857",[5.009766, 53.340303]),
    proj4("EPSG:4326","EPSG:3857",[7.141113, 53.653999]),
    proj4("EPSG:4326","EPSG:3857",[5.822754, 51.037508]),
    proj4("EPSG:4326","EPSG:3857",[3.405762, 51.395350])
]], {
  "fill": "#00000F",
  "stroke": "#00000F",
  "stroke-width": 1
});
var polygon2 = turf.polygon([[
    proj4("EPSG:4326","EPSG:3857",[0.241699, 54.173488]),
    proj4("EPSG:4326","EPSG:3857",[10.162354, 50.908012]),
    proj4("EPSG:4326","EPSG:3857",[8.854980, 50.062208]),
    proj4("EPSG:4326","EPSG:3857",[0.241699, 54.173488])
]], {
  "fill": "#0000FF",
  "stroke": "#0000FF",
  "stroke-width": 1
});
var polygon = turf.intersect(polygon1, polygon2);
polygon.properties = {
  "fill": "#FF0000",
  "stroke": "#FF0000",
  "stroke-width": 1
};
//Convert all polygons back to Leaflet projection
for(var i in polygon1.geometry.coordinates[0]) {
    polygon1.geometry.coordinates[0][i] = proj4("EPSG:3857","EPSG:4326",polygon1.geometry.coordinates[0][i]);
}
for(var i in polygon2.geometry.coordinates[0]) {
    polygon2.geometry.coordinates[0][i] = proj4("EPSG:3857","EPSG:4326",polygon2.geometry.coordinates[0][i]);
}
for(var i in polygon.geometry.coordinates[0]) {
    polygon.geometry.coordinates[0][i] = proj4("EPSG:3857","EPSG:4326",polygon.geometry.coordinates[0][i]);
}
L.mapbox.featureLayer().setGeoJSON(polygon1).addTo(mapLayer);
L.mapbox.featureLayer().setGeoJSON(polygon2).addTo(mapLayer);
L.mapbox.featureLayer().setGeoJSON(polygon).addTo(mapLayer);
mapLayer.setView([52.754260888947776, 5.72100021667583], 8);

Answer 3

这可能是因为（投影）地图上的直线不等于地球表面上的直线（或者换句话说，大地水准面表面上的great circles）。

请参阅some great circles 了解图形说明。

我鼓励你用大圆而不是线来绘制你的大多边形，看看交叉点是否更有意义。