Possible bug in the east-north coordinate system?

[mchitre]

Consider two nearby points on Earth:

p1 = LatLon(1.214914, 103.851194)
p2 = LatLon(1.214986, 103.851527)

The distance between them is:

julia> norm(ECEF(p1) - ECEF(p2))
37.90659161709082

This is the tunneling distance, but given how close the points are, the geodesic distance won't be very different.

The documentation claims that the EastNorth is a local coordinate system that wraps around Earth. I would expect that distances, at least locally, in that system should be correct. Try:

julia> norm(EastNorth(p1) - EastNorth(p2))
2.3117759335122004e7

Clearly this isn't the distance between p1 and p2.

The problem doesn't seem to be in rounding off in floating point due to large east-north coordinates either:

julia> norm(EastNorth(LatLon{BigFloat}(p1)) - EastNorth(LatLon{BigFloat}(p2)))
2.311775933512200348236462687224363807320209393693279196902270161034214710284282e+07

P.S. To handle rounding off problems, it would be nice if the EastNorth conversion could be done with respect to any specified origin for the local coordinate system.

[ettersi]

The problem is this rule:

[...] you can add and subtract EastNorth to / from any Coordinate{2}. The EastNorth coordinate is then interpreted as a translation of the other coordinate.

julia> LatLon(1,2) + EastNorth(3,4)
LatLon{Float64}(1.0000361746684363, 2.0000269535362025)

This rule currently spuriously applies also when both coordinates are EastNorth, and it does so in a way which is not correct:

julia> EastNorth(p1)
EastNorth{Float64}(1.1558080570082083e7, 134338.4378920306)

julia> EastNorth(p2)
EastNorth{Float64}(1.1558117325222772e7, 134346.39927583127)

julia> EastNorth(p1) - EastNorth(p2)
EastNorth{Float64}(2.3116197895304855e7, 268684.8371678619)

This is clearly a bug which needs to be fixed, but I'll first take some time to think through the overall roadmap for coordinate arithmetic. My original intention with this package was only to provide a convenient API for coordinate conversion; I so far haven't thought much about arithmetic, and hence there are a couple of quirks in this area.

[ettersi]

P.S. To handle rounding off problems, it would be nice if the EastNorth conversion could be done with respect to any specified origin for the local coordinate system.

Rounding errors happen the moment you convert your latitude and longitude to Float64.

julia> (1.214915 - big"1.214915") / big"1.214915" |> Float32
-2.7020058f-17

After this, the conversion to EastNorth does not further increase the rounding error in the big-O sense.

julia> f64 = EastNorth(LatLon(1.214914, 103.851194))[]
       bigf = EastNorth(LatLon(big"1.214914", big"103.851194"))[]
       (f64 .- bigf) ./ bigf .|> Float32
(1.1388555f-16, 0.0f0)

(The North error is zero here because the Lat to North conversion internally rounds to Float64.)

I'm therefore not sure whether there's anything we can do here. But also, there's actually not much reason for concern. The maximal rounding error is the world's circumference ( $O(10^7)$ ) times the machine epsilon ( $O(10^{-16})$ ), i.e. on the order of nanometers. Problems only arise if you try to do coordinate arithmetic in Float32, but in that case MapMaths can't really help you as pointed out earlier.