4

I have two features in plain WKT format, the source coordinates being in standard EPSG:4326 projection.

I use these features to run some intersection operations, and I'm troubled by the fact that in WGS 84 the features intersect visually, however if I change projection and render the map in Google/EPSG:3857 I see a shift of around 1500 meters - and the features do not intersect anymore.

I am aware of datum shifts between projections, and as such I do expect there could be smaller distortions. But 1500 meters seems alot?? Also, if I try these features in QGis, they render perfectly - intersecting each other in both WGS84 and Google projections.

Do I need to cater for datum shifts more explicitly in OpenLayers? What else could be causing this?

This is how it looks in OpenLayers using EPSG:4326 enter image description here

And this is the map in OpenLayers in EPSG:3857 enter image description here

And this are the features in QGIS in EPSG:3857 enter image description here

The WKTs of my two features:

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

var feature = '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.45310000000000628, 7.18760000000000066 65.45449999999999591, 7.18920000000000048 65.45740000000003533, 7.19099999999999895 65.45839999999998327, 7.19359999999999999 65.45959999999999468, 7.19669999999999899 65.45959999999999468, 7.19890000000000008 65.46059999999999945, 7.19430000000000103 65.46250000000000568, 7.19409999999999972 65.46519999999998163, 7.19709999999999894 65.46690000000000964, 7.19979999999999976 65.46769999999997935, 7.20379999999999932 65.4678000000000111, 7.20600000000000041 65.46989999999999554, 7.20770000000000088 65.47069999999999368, 7.20989999999999931 65.4713000000000136, 7.20989999999999931 65.47230000000001837, 7.20699999999999985 65.47399999999998954, 7.2110000000000003 65.47520000000000095, 7.21119999999999983 65.47849999999999682, 7.21429999999999882 65.47950000000000159, 7.21839999999999904 65.47950000000000159, 7.21989999999999998 65.48359999999999559, 7.21980000000000022 65.48509999999998854, 7.22240000000000038 65.48689999999999145, 7.22469999999999946 65.48709999999996967, 7.22839999999999971 65.48599999999999, 7.22579999999999956 65.48409999999998377, 7.22769999999999868 65.48250000000001592, 7.22999999999999954 65.4818000000000211, 7.22660000000000124 65.47890000000003852, 7.2275999999999998 65.47719999999995366, 7.23029999999999973 65.47620000000000573, 7.23479999999999901 65.47649999999995885, 7.23630000000000084 65.47499999999999432, 7.23550000000000004 65.47259999999999991, 7.23109999999999964 65.47099999999997522, 7.22130000000000116 65.46989999999999554, 7.21729999999999983 65.46869999999998413, 7.2157 65.46600000000000819, 7.21630000000000038 65.46350000000003888, 7.21359999999999957 65.46289999999999054, 7.21199999999999974 65.45980000000000132, 7.21090000000000053 65.4567000000000121, 7.20970000000000066 65.45420000000001437, 7.2085000000000008 65.45279999999996789, 7.20810000000000084 65.45130000000000337, 7.20690000000000097 65.44900000000001228, 7.20709999999999962 65.44670000000002119, 7.20939999999999959 65.44589999999999463, 7.20779999999999976 65.44259999999999877, 7.20530000000000026 65.43970000000001619, 7.20360000000000067 65.43840000000000146, 7.20150000000000023 65.43619999999998527, 7.19669999999999899 65.43350000000000932, 7.1955000000000009 65.43129999999996471, 7.19740000000000002 65.43019999999998504, 7.19660000000000011 65.42899999999997362, 7.19669999999999899 65.42730000000000246, 7.19460000000000033 65.42500000000001137, 7.19470000000000098 65.4234000000000151, 7.19350000000000112 65.42130000000000223, 7.19040000000000035 65.42050000000000409, 7.19050000000000011 65.41919999999998936, 7.18890000000000029 65.41650000000001342, 7.18590000000000018 65.41509999999999536, 7.18119999999999958 65.411200000000008, 7.17210000000000036 65.40720000000001733, 7.16809999999999903 65.40579999999999927, 7.1650999999999998 65.40440000000000964, 7.1570999999999998 65.40280000000001337, 7.15259999999999962 65.40319999999999823, 7.14909999999999979 65.40229999999999677, 7.14689999999999959 65.40149999999997021, 7.14430000000000032 65.40039999999999054, 7.14090000000000025 65.39839999999998099, 7.13780000000000037 65.397199999999998, 7.13569999999999993 65.39549999999999841, 7.12920000000000087 65.39230000000000587))';

I read the wkt in a standard way:

var format = new ol.format.WKT();

var triangleFeature = format.readFeature(triangle, {
            dataProjection: 'EPSG:4326',
            featureProjection: initialProj
      });

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

      var feature = format.readFeature(feature, {
        dataProjection: 'EPSG:4326',
        featureProjection: initialProj
      });

      var featureLayer = new ol.layer.Vector({
        source: new ol.source.Vector({
          features: [feature]
        })
      })

And I create the map according to standard and other examples:

var initialProj = ol.proj.get('EPSG:3857');
var initialCenter = ol.proj.fromLonLat([10, 65],initialProj);
var initialZoom = 5;      
var map = new ol.Map({
  controls: ol.control.defaults().extend([
          scaleLineControl
        ]),
  target: 'map',
  view: new ol.View({
    center: initialCenter,
    projection: initialProj,
    zoom: initialZoom
  })
});

The only difference between the two images is the value of the initialProj.

Can anyone guide me in the right direction here, or at least provide a good explanation as to why the shift is so huge in OL?

  • If you need to accurately show a route you need to plot the great circle line between the two points. For example a straight line on a map from London to Tokyo takes vastly different routes depending on the map projection. On most projections it appears to go over China just south of Mongolia. But the real direct route crosses Siberia well to the north of Mongolia. See openlayers.org/en/latest/examples/flight-animation.html – Mike Jan 16 at 18:18
  • 1
    It's possible that either OL or QGIS is using ellipsoid-based Mercator equations so that what they say is EPSG:3857 (spherical) is really EPSG:3395 (ellipsoidal). The Y values will differ. The first point of the triangle Y values: 3857 = 10065782.96 versus 3395 = 10026498.99 – mkennedy Jan 16 at 18:45
  • They don't intersect on EPSG:3395 either, the Y values do differ, but they differ consistently for both features. – Mike Jan 16 at 20:51
  • Beyond the spatial meaning of a straight line between two points defined by their latitudes and longitudes, as if they were orthogonal components of a plane coordinate system; or the enormous difference that can occur when confusing World Mercator with Pseudo Mercator (I do not think that QGIS or OL are confused in that sense); what I think the OP is asking is: why his code is not transforming the geometries accordingly (assuming that QGIS transforms them correctly)? I particularly think that the coordinates of their vertices are badly transformed. – Gabriel De Luca Jan 16 at 21:51
  • 1
    No system would make major errors in transforms between the two most widely used projections. If splitting an EPSG:4326 line into 100 segments to reproject it gives identical results to QGIS then QGIS must be doing something similar in non-4326 projections. I have now colour coded the lines and added a background layer to locate the land vertices. I've also done a version with the vectors in their native EPSG:4326 projection mikenunn.16mb.com/demo/reprojected-lines-4326.html where the red line is straight (as it should be) while in EPSG:3857 it isn't. It will be the same in QGIS. – Mike Jan 17 at 0:29
4

Here's the great circle geometry with EPSG:4326 enter image description here

and the great circle geometry with EPSG:3857/3395 enter image description here

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

  • Thank you mike for all effort. I'll consider this one the correct answer. – Eivind Rønnevik Jan 17 at 7:49
  • 1
    You question does raise a wider issue about how to display (or even how governments police), for example, international maritime boundaries. They are defined as a series of often widely spaced Lat/Lon points but do the international agreements specify if they are joined in flat WGS84 or great circles geometry? (it certainly won't be based on a projection conceived by Google). mikenunn.16mb.com/demo/uk-eez.html – Mike Jan 17 at 11:41
  • Starting from an absolute ignorance of international agreements, I suspect that the points should be joined with (curved) lines that keep their azimuth constant throughout their trajectory. Visible as straight lines in conformal projections (at the expense of an incorrect appreciation of their distances). It is not the (curved) shortest path between the points, it is the constant bearing path between them. Except for the reduction to the sphere, it would be more appropriate to draw them as straight lines in web mercator than in flat lon/lat, or great circles. – Gabriel De Luca Jan 17 at 16:49
5

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:

triangle

  • Thank you for the time spent looking into this! Really appreciated. I realize there's a lot more to it than I initially expected, but both your efforts have contributed to give me a better understanding of the challenges. – Eivind Rønnevik Jan 17 at 7:47
2

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;
    },
    ....
});
  • The area of the geometry isn't necessarily relevant, certainly not the spherical area as the problem get worse as latitude increases. Even large polygons are not affected if they are already well segmented, but long straight lines (e.g. in long thin polygons) can be a problem regardless of the polygon's area . If instead of defining a constant for segmentations calculating a suitable value (rounded up so the minimum is one) based on line.getLength() for each section of the linear ring would be more efficient, and in non-problem cases would simply clone the existing geometry. – Mike Jan 18 at 12:28
  • 1
    Thanks Mike, valuable input. In my scenario the polygons in questions are drawn by the user, hence I know I can expect very low-segmented polygons - and the area works ok to detect it. However as a more proper and generic implementation your suggestion is far better. Will have a go at it! – Eivind Rønnevik Jan 19 at 13:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.