¿Por qué hay un desplazamiento/distorsión de ~1,5 km en la transformación en OpenLayers 5.3.0 al cambiar de EPSG:4326 a EPSG:3857?

Question

Tengo dos características en formato WKT simple, las coordenadas fuente estando en proyección estándar EPSG:4326.

Utilizo estas características para realizar algunas operaciones de intersección, y me preocupa el hecho de que en WGS 84 las características se intersecan visualmente, sin embargo si cambio la proyección y renderizo el mapa en Google/EPSG:3857 veo un desplazamiento de alrededor de 1500 metros - y las características ya no se intersecan.

Soy consciente de los cambios de datum entre proyecciones, y como tal espero que pueda haber distorsiones más pequeñas. ¿Pero 1500 metros parece mucho? Además, si pruebo estas características en QGIS, se renderizan perfectamente, intersecándose en ambas proyecciones de WGS84 y Google.

¿Necesito considerar los cambios de datum de forma más explícita en OpenLayers? ¿Qué más podría estar causando esto?

Así es como se ve en OpenLayers usando EPSG:4326

Y este es el mapa en OpenLayers en EPSG:3857

Y estas son las características en QGIS en EPSG:3857

Los WKTs de mis dos características:

var triangulo = 'Polygon ((5.33973571664250635 66.68106075732242743, 10.45882119706011792 63.36917710864676678, 15.66616677196768492 67.09677082116945712, 5.33973571664250635 66.68106075732242743))';

var caracteristica = 'Polygon ((7.12920000000000087 65.39230000000000587, 7.12659999999999982 65.39210000000002765, 7.12289999999999957 65.39339999999998554, 7.1227999999999998 65.39440000000001874, 7.12539999999999996 65.39519999999998845, 7.13019999999999943 65.3972999999999729, 7.1344000000000003 65.40059999999999718, 7.13529999999999998 65.40160000000000196, 7.13740000000000041 65.40339999999997644, 7.13730000000000064 65.40430000000000632, 7.14029999999999898 65.40609999999998081, 7.14290000000000003 65.40690000000000737, 7.14369999999999994 65.40879999999998518, 7.14360000000000106 65.41050000000001319, 7.14520000000000088 65.41249999999999432, 7.14550000000000018 65.41440000000000055, 7.14860000000000007 65.41480000000004225, 7.14820000000000011 65.42059999999997899, 7.15309999999999935 65.42199999999999704, 7.15469999999999917 65.42430000000001655, 7.15589999999999904 65.42640000000002942, 7.1631999999999989 65.4304000000000201, 7.16530000000000111 65.43350000000000932, 7.16739999999999977 65.43559999999999377, 7.17300000000000093 65.43889999999998963, 7.17550000000000043 65.44089999999999918, 7.17879999999999985 65.44499999999999318, 7.1791999999999998 65.44570000000001642, 7.18089999999999939 65.44710000000000605, 7.18489999999999895 65.44889999999998054, 7.18480000000000008 65.45040000000000191, 7.18510000000000026 65.45179999999999154, 7.18590000000000018 65.453100000000

Answer 1

De manera efectiva y contrario a lo que intuí al principio, las coordenadas de los puntos están bien proyectadas.

QGIS está mostrando la línea en otra posición, porque la capa está asignada con el CRS EPSG:4326 y se está re-proyectando sobre la marcha a EPSG:3857.

Esto hace que QGIS calcule la ruta curva (o proyectada) desde los vértices de la geometría hasta los bordes de la pantalla, y a partir de ahí dibuja un segmento de línea recta. Por lo tanto, en cada nivel de zoom, la línea se encuentra en otro lugar, ya que cambia la ubicación de los bordes de la pantalla.

La solución para ver la línea recta en su ubicación correcta (la unión recta entre las coordenadas proyectadas de los vértices del polígono) en QGIS, es exportar la capa al CRS EPSG:3857.

He creado un sistema similar al trabajo de Mike, pero en QGIS. Las líneas de gran círculo fueron aproximadas desde una proyección gnomónica personalizada, segmentadas en 50 partes, simbolizadas en verde. Las líneas de triángulo recto son las negras. Y una densificación de 10 segmentos de cada una fue simbolizada en rojo.

Adjunto una captura de pantalla del comportamiento de la representación en tiempo real de QGIS en este caso particular:

Answer 2

Aquí está la geometría del círculo grande con EPSG:4326

y la geometría del círculo grande con EPSG:3857/3395

por lo tanto parece que QGIS está produciendo una representación exacta de la línea recta EPSG:4326 en una proyección EPSG:3857. Diferente a OpenLayers pero ambos son inexactos en ambas proyecciones, aunque QGIS es consistentemente inexacto.

Este es el código que utilicé para crear el triángulo geodésico. Requiere https://api.mapbox.com/mapbox.js/plugins/arc.js/v0.1.0/arc.js También se incluye el código para reproducir el comportamiento de la geometría QGIS EPSG:4326 en cualquier proyección.

var triangleFeature = format.readFeature(triangle, {
            dataProjection: 'EPSG:4326',
            featureProjection: 'EPSG:4326' // Mantener coordenadas en EPSG:4326 necesarias para construir polígono geodésico
      });

var triangleCoords = triangleFeature.getGeometry().getCoordinates()[0];

var polyCoords = [];
for (var n=0; n

`

http://mikenunn.16mb.com/demo/reprojected-lines..html http://mikenunn.16mb.com/demo/reprojected-lines-4326.html

`

Answer 3

Basado en la entrada de @mike y @gabriel, llegué a la siguiente solución para manejar esto en el lado del cliente al renderizar las características (sin modificar la geometría de la característica original):

Básicamente, hago uso del objeto Style y la función geometry (que devuelve una geometría para renderizar). En tiempo de ejecución, revisará si el área de un polígono es mayor a cierto tamaño; en ese caso, creo un nuevo polígono con líneas segmentadas, lo que permite que la vista proyectada renderice una representación que para el ojo humano es percibida como similar a la renderización de WGS 84.

Por supuesto, el tamaño del área que activa la segmentación, así como el número de segmentos, deben considerarse en cada escenario.

Un fragmento de código:

var style = new ol.style.Style({
    geometry: function (feature) {
        var geometry = feature.getGeometry();

        if(map.getView().getProjection().getCode() === 'EPSG:4326' ){
            return geometry;
        }

        if (geometry.getType() === 'Polygon'){
            var area = ol.sphere.getArea(geometry, {projection: map.getView().getProjection()} );

            if (area > 30000000000) {
                geometry = feature.getGeometry().clone();
                geometry.transform(map.getView().getProjection(),'EPSG:4326');

                var numSegmentations = 20;
                var coords = [];
                var geomCoords = geometry.getLinearRing(0).getCoordinates();
                for (var n=0; n