Differences in angle measurements in Web Mercator

Question

I am trying to measure different angles on 2D map using lat/lon coordinates taken from Google Maps, and Java library JTS.

I have noticed a weird thing: my measurements are always distorted, except angles that are perpendicular to meridians and parallels. It gets as high as 25 degrees when the angle sides are at an angle of 45 degrees to the parallalels, example here:

The result is 68.97 and the library I use uses the law of cosines to measure the angle. I also double checked the result using the method provided here, and the difference was still there.

Here's the Java code:

double point2X = 2.35956;
double point2Y = 48.86054;

double point1X = 2.35922;
double point1Y = 48.86215;

double fixedX = 2.3582;
double fixedY = 48.86124;


double angle1 = Math.atan2(point1Y - fixedY, point1X - fixedX);
double angle2 = Math.atan2(point2Y - fixedY, point2X - fixedX);

return Math.toDegrees(angle1 - angle2);

Or as done using JTS:

Coordinate point1 = new Coordinate(2.35922, 48.86215);
Coordinate point2 = new Coordinate(2.35956, 48.86054);
Coordinate fixed = new Coordinate(2.3582, 48.86124);

return Math.toDegrees(Angle.angleBetween(point1, fixed, point2));

And here is a JS fiddle showing the calculations:

https://jsfiddle.net/hrg04v3n/

I found out that Web Mercator is non-conformal, thus it doesn't preserve angles, but does it get as bad for even such small areas as presented above?

Answer 1

Geographic (lat/lon) coordinates are already angular measures, so you can't compute a simple arctangent with the coordinates based on atan2(delta_y,delta_x) as if they were Cartesian values.

Using a port of the US Geodetic Survey FORTRAN code to solve the inverse (aka reverse) geodetic problem, and using the coordinate values in your Java code (assuming WGS84 datum), I got the following:

p1                     |p2                      |gcbearing
POINT (2.3582 48.86124)|POINT (2.35922 48.86215)|36.4845
POINT (2.3582 48.86124)|POINT (2.35956 48.86054)|127.9567

The difference between the bearings is 91.4722 degrees.

Reprojecting the points to WGS84 UTM Zone 31N and using Cartesian bearing, I got:

p1                            |p2                            |bearing
POINT (452926.994 5412229.237)|POINT (453002.658 5412329.763)|36.9681
POINT (452926.994 5412229.237)|POINT (453026.088 5412150.583)|128.4401

The difference here is 91.4720 degrees (since Paris is on the left side of the UTM cylinder the bearings are slightly different, as well)

Finally, reprojecting into a custom WGS84 datum flavor of the Europe Albers Equal Area Conic projection (central meridian 10.0, standard parallels 43.0,62.0), I got:

p1                            |p2                             |bearing
POINT ( -553711.962 2098386.8)|POINT (-553627.667 2098480.823)|41.8774
POINT ( -553711.962 2098386.8)|POINT (-553621.972 2098298.257)|134.5356

with a difference of 92.6582 degrees.

You cannot rely on a Web Mercator canvas for whether the intersection "looks like" a right angle because of the significant Y-axis distortion in Web Mercator, even at only 49N.

I don't have any experience in JTS, but it apparently only has Cartesian (not geodetic) functions, so you'd need to reproject to compute bearing changes (based on the above I'd recommend using UTM for relatively small areas).