Implement "wrapping" local coordinates?

[ettersi]

There are a few occasions where it would be handy to have local coordinates that "wrap" around the surface of the earth. For example:

• When putting axis ticks in meters / kilometers on a WebMercator map:

  yticks = [
      convert(
          WebMercator, 
          cwrap(
              WrappingEastNorth(map_origin), 
              (0u"km", y), 
          )
      ).y
     for y in (-5:5)u"km"
  ]

• We work with devices that communicate acoustically under water. In shallow bodies of water, sound is constrained to travel along the surface of the earth rather than tunnel through the seabed, so in this case working in a wrapping coordinate system is more accurate.

I must admit that I am not aware of any reference for such a coordinate system, so this proposal is quite speculative. Curious to hear what other people might think.

[juliohm]

I believe that we need to distinguish between periodic and aperiodic dimensions as we do in Meshes.jl. If we are talking about LatLon coordinates, these can clearly wrap around with the appropriate algebraic structure. On the other hand, if we are talking about Projected coordinates, there is not really a way to wrap around. Instead, we should convert back to LatLon, work in this periodic space, and then project back.

Notice that Projected coordinates must satisfy bounds, catalogued in their EPSG codes. These bounds are a function of the LatLon bounds of the projection. Operating on the projected coordinates directly is risky, and we should avoid that.

[ettersi]

I believe I may have chosen my words poorly. By "wrapping coordinates", I mean coordinates that wrap around the surface of the earth like a paper wraps around a stone, not "wrapping" as in periodic boundaries. So what I'm asking for is not to change the behaviour of any existing coordinate type, but to add a new type such that LatLon(0,0) + WrappingEastNorth(1e6,0) gives me a point close to LatLon(0, 90) rather than somewhere out in space.

[juliohm]

How does it differ from #48 then? You mean that we should be able to consider this periodic behavior when adding the ENU shifts? That will be automatic given the constructor of LatLon.

[ettersi]

ENU coordinates are usually assumed to be cartesian, while what I'm proposing here is not. For illustration, assuming the earth is a sphere with radius 1:

Cartesian(1, 0, 0) + CartesianEastNorth(1, 0) -> Cartesian(1, 1, 0)
Cartesian(1, 0, 0) + WrappingEastNorth(1, 0) -> Cartesian(cos(1), sin(1), 0)

[juliohm]

So it is a matter of terminology then. We thought that ENU was the wrapping case where LatLon + ENU would still produce LatLon on the sphere. Why do we need a CartesianENU? Isn't Cartesian enough?

[ettersi]

Imagine I have a GPS, a compass, a plumb line and a ranging device. Then I can compute the location of any other visible object using something like

gps + convert(CartesianENU, Spherical(range, 90u"°" - slope, azimuth + 90u"°"))

(Not sure what the syntax for the Spherical-to-CartesianENU conversion should be, but thanks a ton for the Spherical type 😉 )

[juliohm]

Not clear yet the meaning of CartesianENU, can you provide a MWE with actual numbers to see if we are talking about the same ideas?

[ettersi]

Cartesian:

julia> using Geodesy
       x = LLA{Float64}(0,0,0)
       d = ENU{Float64}(10e3,0,0);

julia> LLAfromENU(x, wgs84)(d)
LLA(lat=0.0°, lon=0.08983145480493704°, alt=7.839274896367615) # Almost 8m off the ground

julia> ECEFfromENU(x, wgs84)(d)
ECEF(6.378137e6, 10000.0, 0.0) # y component is exactly 10km

Wrapping:

julia> using MapMaths
       x = LatLon(0,0)
       d = EastNorth(10e3,0);

julia> x + d
LatLon{Float64}(-1.43e-322, 0.08983152841195213) # Still exactly at sea level

julia> ECEF(x + d)
ECEF{Float64}(6.378129160721892e6, 9999.99590304674, -3.130123e-317) # y component is slightly less than 10km

[juliohm]

I believe that this issue is a dup of #48 in the sense that we will only need a ENU type that wraps around the ellipsoid. The CartesianENU is not different than a normal Cartesian.