The Geodesy.jl package is very useful for transforming between Latitude-Longitude and local coordinate systems. However, the only implemented local coordinate system is "ENU", i.e. East-North-Up. For the runway use-case, usually we want coordinates to be expressed in runway-centric coordinates "XYZ", i.e. "Alongtrack-Crosstrack-Up". This package extension implements that, taking inspiration from JuliaGeo/Geodesy.jl#71 implementing a North-East-Down coordinate system.
Usage should be similar to the Geodesy.jl package. For our use-case, probably we want something like
using GeodesyXYZExt; using Geodesy
DATUM=wgs84;
a, b, c = 27.1127, 109.3497, 5.0
origin_lla::LLA = LLA(a, b, c)
bearing::Float64 = deg2rad(45) # measured clockwise from north
# let's bind the transformations for easier handling
let XYZfromLLA_ = XYZfromLLA(origin_lla, bearing, DATUM),
LLAfromXYZ_ = LLAfromXYZ(origin_lla, bearing, DATUM)
origin::XYZ = XYZfromLLA_(LLA(a, b, c))
corners = XYZ[
origin + XYZ(0., 10., 0),
origin + XYZ(0., -10., 0)
]
LLAfromXYZ_.(corners)
@test LLAfromXYZ_(origin) ≈ origin_lla
endHowever, often we just want to fix one origin and bearing (and datum), which is also supported:
using GeodesyXYZExt
GeodesyXYZExt.fixdatum!(wgs84)
GeodesyXYZExt.fixorigin!(LLA(27.1127, 109.3497, 5.0))
GeodesyXYZExt.fixbearing!(deg2rad(-90))
enu = ENU(1.0, 0.0, 1.0)
xyz = XYZ(enu)
enu_ = ENU(xyz)
@assert enu_ ≈ enu # works :)
lla = LLA(xyz)
xyz_ = XYZ(lla)
@assert xyz_ ≈ xyz # works :)Notice that we have disabled adding vectors in different coordinate systems, as this is almost often wrong. Please convert to the same coordinate system.
@test_throws CoordinateSystemError XYZ(1., 2., 3.) + ENU(1., 2., 3.)Note that there's still many functions that are overloaded and will evaluate (wrongly). Due to Julia's dispatch system, it's not easy to "turn them all off". E.g.
XYZ(1., 2., 3.) == ENU(1., 2., 3.) # `true` (!)This package also provides an extension that should allow working with Unitful.jl quantities.
using Unitful; using Unitful.DefaultSymbols
Meters = typeof(1.0m)
XYZ(1.0m, 2.0m, 3.0m) == XYZ(1., 2., 3.) * 1m
ustrip.(uconvert.(mm, XYZ(1.0m, 2.0m, 3.0m))) == XYZ(1000., 2000., 3000.)
@assert typeof(XYZ(1.0m, 2.0m, 3.0m)) == XYZ{Meters}