Subtraction of LatLon

[mchitre]

Currently, if I subtract two LatLon coordinates, I get a LatLon with the differences:

julia> LatLon(1.214914, 103.851194) - LatLon(1.214986, 103.851527)
LatLon{Float64}(-7.199999999984996e-5, -0.00033299999999769625)

While this may seem reasonable, it is incorrect. The differences between latitudes is not a latitude, and the differences between longitude is not a longitude. This becomes important when the latitude-longitude differences later need to be converted to east-north distances:

julia> en = EastNorth(LatLon(1.214914, 103.851194) - LatLon(1.214986, 103.851527))
EastNorth{Float64}(-37.069390433874574, -7.961347859138245)

julia> norm(en)
37.91467745970203

which is quite different from:

julia> norm(ECEF(LatLon(1.214914, 103.851194)) - ECEF(LatLon(1.214986, 103.851527)))
37.90659161709082

The essential problem is that the actual latitude-longitude is lost in taking the difference, and then the east-north conversion wrongly treats the differences as being located at the origin of the latitude-longitude coordinate system.

Perhaps we need a different data type to represent differences in latitude-longitude, rather than returning them as a latitude-longitude pair?

[ettersi]

The fundamental problem here is that there isn't just one way to do coordinate arithmetic. For example, I'd argue that the following operation can have two very different results, both of which are perfectly reasonable in their own way.

x = LatLon(0,0) + 2*(LatLon(0,180) - LatLon(0,0))

• x = LatLon(0,360): LatLon(0,180) - LatLon(0,0) represents a rotation halfway around the earth, so two halfway-rotations should make a full rotation around the earth.
• x = (LatLon(0,180), Alt(1.3e7)): ECEF(LatLon(0,180)) - ECEF(LatLon(0,0)) = ECEF(1.3e7, 0, 0), so continuing in this direction should take us out to space.

Things get even trickier if we take differences at one latitude and then apply them at a different latitude. For example, what should the result of this expression be?

y = LatLon(45,0) + 2*(LatLon(0,180) - LatLon(0,0))

I'd say all of the following are sensible answers:

• LatLon(45, 360): We're still doing a full rotation, just at a different latitude.
• (LatLon(0 < lat < 45, 180), Alt(>0)): The result if we were doing the arithmetic in ECEF.
• LatLon(45, >360): A trip around the earth is shorter at latitude 45 than at latitude 0, so 2*(LatLon(0,180) - LatLon(0,0)) should be more than a full rotation.

Finally, what if even the difference is not taken at a single latitude?

z = LatLon(45,0) + 2*(LatLon(30,180) - LatLon(0,0))

Should LatLon(30,180) - LatLon(0,0) represent "30 degrees north and 180 degrees east"? If so, in what order? (Depending on how exactly you define "going north" and "going east", the two operations may not commute.) Alternatively, LatLon(30,180) - LatLon(0,0) could represent "go X km in the direction of the great circle which leaves the starting point at heading Y", which once again would give you a different answer...

The rationale behind the current API is that while its results may not always be particularly useful, at least it is very easy to explain what is going on: coordinate arithmetic is always executed elementwise; you as the user are responsible for choosing the coordinates such that elementwise arithmetic does what you want!

We actually have already seen an instance of this: in the current form of MapMaths, we have

julia> LatLon(45,0)[] .+ 2 .* (LatLon(0,180)[] .- LatLon(0,0)[])
(45.0, 360.0)

If instead of the rotation you want cartesian arithmetic, well then convert your arguments to ECEF beforehand:

julia> ECEF(LatLon(45,0))[] .+ 2 .* (ECEF(LatLon(0,180))[] .- ECEF(LatLon(0,0))[])
(-2.0994957121151067e7, 0.0, 4.48734840886592e6)

Finally, if instead of counting rotations you want to count distance, then convert your arguments to EastNorth:

julia> EastNorth(LatLon(45,0))[] .+ 2 .* (EastNorth(LatLon(0,180))[] .- EastNorth(LatLon(0,0))[]) |> EastNorth |> LatLon
LatLon{Float64}(44.999999999999986, 508.2641127930218)

All this is not to say that the current system is perfect. For example, I was quite shocked to see this:

julia> norm(ECEF(LatLon(1.214914, 103.851194)) - ECEF(LatLon(1.214986, 103.851527)))
37.90659161709082

julia> norm(EastNorth(LatLon(1.214914, 103.851194))[] .- EastNorth(LatLon(1.214986, 103.851527))[])
37.607499239327225 

Clearly, the difference between tunnel and surface distance cannot be 30cm over a range of 37m, and even more importantly, the surface distance cannot be less than the tunnel distance. What is happening here is that one point is about 8m further north than the other, and since East measures the distance to the prime meridian at the given latitude, this has a significant impact on the East component. This clearly needs to be fixed somehow. Maybe there should be a conversion to East relative to a meridian other than the prime meridian? I'm open to suggestions...

Apart from "hard" / "numeric" bugs like the one described just now, I guess there is also plenty of scope for more syntactic sugar. For example, it might be nice to have a macro

@coordinate(ECEF, x + 2*(y - z))

which expands to

ECEF(x) + 2*(ECEF(y) - ECEF(z)))

Finally, broken / missing methods like https://github.com/subnero1/MapMaths.jl/issues/3 and *(::Number, ::Coordinate) should be implemented...

[mchitre]

I see the primary use for math on lat-lon is to stay on the surface or earth. In most cases, a locally flat coordinate system around a given lat-lon is sufficient for this. In cases, where the distances are large, accounting for the earth's curvature is nice.

Thing's I'd typically want to do:

• measure distances between lat-lon pairs
• move from a lat-lon by certain distance along a particular bearing and get lat-lon at that location
• convert lat-lon to a locally flat ENU coordinate system in meters, and convert back to lat-lon

[ettersi]

I believe all of these objectives can be achieved without obstructing any other use case by defining coordinate arithmetic as follows.

• Coordinates themselves don't allow any arithmetic, i.e. LatLon(1,2) - LatLon(3,4) is a method error.
• We add functions cstrip(CoordType, coord) -> SVector and cwrap(CoordType, svec) -> Coord which convert locations to and from SVectors of plain numbers representing the coordinates in the given coordinate system.

The three actions that you mentioned could then be implemented as follows:

• measure distances between lat-lon pairs: norm(cstrip(ECEF, x) - cstrip(ECEF, y))
• move from a lat-lon by certain distance along a particular bearing and get lat-lon at that location: LatLon(cwrap(ENU(x), cstrip(ENU(x)) + r * SVector(cos(alpha), sin(alpha), 0))).
• convert lat-lon to a locally flat ENU coordinate system in meters, and convert back to lat-lon: cstrip(ENU(x), y), cwrap(ENU(x), d).

I hope I can get CoordRefSystems to adopt these conventions...