Why is there a transformation shift/distortion of ~1.5km in OpenLayers 5.3.0 going from EPSG:4326 to EPSG:3857?

Effectively and contrary to what I intuited in the beginning, the coordinates of the points are well proyected.

QGIS is showing the line in another position, because the layer is assigned with the EPSG:4326 CRS and is reprojecting on the fly to EPSG: 3857.

This causes QGIS to calculate the curved (or projected) path from vertices of the geometry to the edges of the screen, and from there draw a straight line segment. Therefore, at each zoom level, the line is located in another place, since it changes the location of the edges of the screen.

The solution to see the straight line in its correct location (the straight union between the projected coordinates of the vertices of the polygon) in QGIS, is to export the layer to the EPSG:3857 CRS.

I created a system similar to Mike's job, but in QGIS. The great circle lines were approximated from a custom gnomonic projection, segmentized in 50 parts, symbolized in green. The straight triangle lines are the black ones. And a densify of 10 segments of each one was symbolized in red.

I attach a screencast of the behavior of QGIS rendering on the fly in this particular case:

Here's the great circle geometry with EPSG:4326

and the great circle geometry with EPSG:3857/3395

so it appears QGIS is producing an accurate representation of the EPSG:4326 straight line in an EPSG:3857 projection. Different to OpenLayers but both are inaccurate in both projections, although QGIS in consistently inaccurate.

This is the code I used to create the geodesic triangle. It requires https://api.mapbox.com/mapbox.js/plugins/arc.js/v0.1.0/arc.js Also included is code to reproduce the QGIS EPSG:4326 geometry in any projection behaviour

var triangleFeature = format.readFeature(triangle, {
            dataProjection: 'EPSG:4326',
            featureProjection: 'EPSG:4326' // Keep coordinates in EPSG:4326 needed to construct geodesic polygon
      });

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

var polyCoords = [];
for (var n=0; n<triangleCoords.length-1; n++) {
      // create an arc circle between the two locations
      var arcGenerator = new arc.GreatCircle(
            {x: triangleCoords[n][0], y: triangleCoords[n][1]},
            {x: triangleCoords[n+1][0], y: triangleCoords[n+1][1]});
      var arcLine = arcGenerator.Arc(100, {offset: 10});
      if (arcLine.geometries.length === 1) {
            polyCoords = polyCoords.slice(0,-1).concat(arcLine.geometries[0].coords);
      }
}

var newPoly = new ol.geom.Polygon([polyCoords]);
newPoly.transform(ol.proj.get('EPSG:4326'), ol.proj.get(initialProj));

var featureGeo = new ol.Feature({
      geometry: newPoly
});

var featureGeoLayer = new ol.layer.Vector({
      source: new ol.source.Vector({
            features: [featureGeo]
      })
});

var poly4326Coords = [];
for (var n=0; n<triangleCoords.length-1; n++) {
      var line = new ol.geom.LineString([triangleCoords[n],triangleCoords[n+1]]);
      poly4326Coords = poly4326Coords.slice(0,-1);
      for (var i=0; i<=100; i++) {
          poly4326Coords.push(line.getCoordinateAt(i/100));
      }
}

var newPoly4326 = new ol.geom.Polygon([poly4326Coords]);
newPoly4326.transform(ol.proj.get('EPSG:4326'), ol.proj.get(initialProj));

var feature4326 = new ol.Feature({
      geometry: newPoly4326
});

var feature4326Layer = new ol.layer.Vector({
      source: new ol.source.Vector({
            features: [feature4326]
      })
});

triangleFeature.getGeometry().transform('EPSG:4326', initialProj);  // Now convert triangle to view projection

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

Based on the input from @mike and @gabriel I came up with the following solution to handle this clientside when rendering the features (so not modifying the geometry of the original feature):

Basically I make use of the Style object and the geometry-function (which returns a geometry to render). Runtime it will check if the area of a polygon is over a certain size; If so - I create a new polygon with segmented lines: thus allowing the projected view to render a representation which by the human eye is perceived to be similar to the WGS 84 rendering.

Of course the area size which triggers segmentation, as well as number of segments, must be considered in each scenario.

A code snippet:

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<geomCoords.length-1; n++) {
                      var line = new ol.geom.LineString([geomCoords[n],geomCoords[n+1]]);
                      coords = coords.slice(0,-1);
                      for (var i=0; i<=numSegmentations; i++) {
                          coords.push(line.getCoordinateAt(i/numSegmentations));
                      }
                }
                geometry = new ol.geom.Polygon([coords]);
                geometry.transform('EPSG:4326',map.getView().getProjection());
                return geometry;
            }

        }
        return geometry;
    },
    ....
});