Пересечение многоугольника на карте показывает смещение

Question

Я пытаюсь найти область пересечения между двумя многоугольниками, нарисованными на карте. Я пересекаю метод TurfJS, чтобы найти пересечение между двумя многоугольниками.

Он работает хорошо для небольших областей, но для больших многоугольников он начинает показывать некоторое смещение в точках пересечения, смещение растет с увеличением расстояния. Это также видно только для полигонов с наклонными линиями (вертикальные и горизонтальные пересечения линий, похоже, хорошо работают).

Я создал 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.js и LeafletJs используют разные проекционные системы. Поэтому я вручную проецировал координаты в EPSG: 3857, прежде чем перейти к методам turfjs, а затем преобразовать результаты 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 2

Это может быть связано с тем, что прямые линии на (проецируемой) карте не равны прямым линиям на поверхности земли (или, другими словами, great circles на поверхности геоида).

См. some great circles для графического объяснения.

Я призываю вас нарисовать большие многоугольники с большими кругами вместо линий и посмотреть, имеют ли пересечения смысл.

Answer 3

Непонятно мне, что именно торф делает, когда он вычисляет пересечения (они не лежат на больших кругах между вершинами многоугольника, как подозревал Иван), но проблема может быть решена путем интерполирования многоугольников вдоль больших окружных дуг , затем используя 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, который будет интерполировать функции GeoJSON LineString, которые имеют свойство geodesic, установленное на "true". Большая часть этого кода потрачена на преобразование 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 может быть излишним здесь, но оно работает. Я подозреваю, что любая достаточно тонкая интерполяция на самом деле будет работать нормально (т. Е. Торф будет вычислять пересечения, которые падают вдоль краев многоугольника при наложении на карту), но этот метод имеет то преимущество, что он географически корректен.