Computing Geodesics

I'd like to be able to provide two points given latitude and longitude and compute the "as the crow flies" distance between them.

locations_latlon = Dict{Symbol, LatLon}(
    :washington_dc => LatLon("38.9072N", "77.0369W"),
    :london_uk => LatLon("51.5074N", "0.1278W"),
    :paris_fr => LatLon("48.8555N", "2.3522E"),
    :lyon_fr => LatLon("45.7597N", "4.8422E"),
    :tokyo_jp => LatLon("35.6804N", "139.7690E"),
    :buenosaires_ar => LatLon("34.6037S", "58.3816W"),
    :capetown_sa => LatLon("33.9249S", "18.4241E")
)

As shown above, I need some way of storing the lat and lon data, so I created a struct called LatLon and wrote a constructor. Since Most of the lat and lon data I have comes in string form with a cardinal direction stamped on the end, I'll also write a helper function to convert the string form to a float. Any South and West addresses will be negative, so I'll scale the parsed value by either \(1\) or \(-1\) depending on if I need to add a sign or not.

Another approach would be to multiply the parsed value by \((-1)^k\) where the parity of \(k\) would influence whether the final term remained positive or became negative, but we only have two states to worry about and a multiplicative operation will be faster than an exponent operation.

function transformlatlon(latlon::String)
    sgn = (latlon[end] ∈ ['S', 'W']) ? -1 : 1
    sgn * parse(Float64, latlon[1:end-1])
end

struct LatLon{T<:Real}
    lat::T
    lon::T
end

function LatLon(lat::String, lon::String)
    (lat, lon) = map(transformlatlon, [lat, lon])
    LatLon(lat, lon)
end

I'd like to be able to decompose a LatLon object quickly, e.g.

(lat, lon) = LatLon("38.9072N", "77.0369W")

which requires a simple implementation of Base.iterate. Since there are only two states to worry about, it's pretty straight-forward.

function Base.iterate(LL::LatLon, state=1)
    state > 2 && return nothing
    state == 1 ? (LL.lat, 2) : (LL.lon, 3)
end

Finally, we'll implement the Haversine Distance formula to compute the geodesic or "great circle" distance. Since I'm storing my latitude and longitude values in a dictionary, I'll leverage Julia's method dispatching and write a second wrapper function where I can just pass in symbols with the names of locations.

Given points \((\rho, \phi_1, \theta_1)\) and \((\rho, \phi_2, \theta_2)\), where \(\rho = 6371\) kilometers, \(\Delta \phi = \phi_2 - \phi_1\), and \(\Delta \theta = \theta_2 - \theta_1\), the Haversine Distance formula is defined as,

\[2 \cdot \rho \cdot \tan^{-1}\Bigg(\sqrt{\frac{\sin^2\Big(\frac{\Delta \phi}{2}\Big) + \cos(\phi_1)\cos(\phi_2) \sin^2\Big(\frac{\Delta \theta}{2}\Big)}{1 - \sin^2\Big(\frac{\Delta \phi}{2}\Big) - \cos(\phi_1)\cos(\phi_2) \sin^2\Big(\frac{\Delta \theta}{2}\Big)}}\Bigg)\]

function haversine(p1::LatLon, p2::LatLon)
    lat1, lon1 = p1
    lat2, lon2 = p2
    radius = 6371 # in kilometers
    
    dlat = deg2rad(lat2 - lat1)
    dlon = deg2rad(lon2 - lon1)
    
    a = (sin(dlat / 2) * sin(dlat / 2) + 
        cos(deg2rad(lat1)) * cos(deg2rad(lat2)) *
        sin(dlon / 2) * sin(dlon / 2))
    c = 2 * atan(sqrt(a), sqrt(1 - a))
    radius * c
end

function haversine(p1::Symbol, p2::Symbol)
    haversine(locations_latlon[p1], locations_latlon[p2])
end

Now I can easily decompose locations.

(lat, lon) = locations_latlon[:washington_dc]

@show lat
> lat = 38.9072

@show lon
> lon = -77.0369

And I can compute distances by name.

@show haversine(:washington_dc, :london_uk)
> haversine(:washington_dc, :london_uk) = 5897.618855872552

@show haversine(:lyon_fr, :paris_fr)
> haversine(:lyon_fr, :paris_fr) = 392.0501333501885

Spherical Coordinates

Another common coordinate system for geographic locations is spherical coordinates. If we think of polar coordinates, which transform from cartesian coordinates with,

\[\begin{aligned} x &= r \cos(\theta) \\ y &= r \sin(\theta) \end{aligned}\,,\]

where \(r\) is the radius and \(\theta = \tan^{-1}(\frac{y}{x})\). We can expand to cylindrical coordinates by adding a \(z = z\) component.

Finally, spherical coordinates are transformed with,

\[\begin{aligned} x &= \rho \cos(\phi)\sin(\theta) \\ y &= \rho \sin(\phi)\sin(\theta) \\ z &= \rho \cos(\theta)\,. \end{aligned}\]

Latitude, \(\ell_1\), is presented from -90 to 90 degrees and longitude, \(\ell_2\), is presented from -180 to 180 degrees. Therefore converting to radians simply requires a small adjustment.

Let \(\xi = \frac{\pi}{180}\).

\[\begin{aligned} \rho &= \rho \\ \phi &= \frac{\pi}{2} - \ell_1\cdot\xi\\ \theta &= \ell_2\cdot\xi\,. \end{aligned}\]

struct SphericalCoordinates{T<:Real}
    ρ::T
    ϕ::T
    θ::T
end

function SphericalCoordinates(geo::LatLon)
    ρ = 6.371e3 # in kms
    ϕ = deg2rad(90-geo.lat)
    θ = deg2rad(geo.lon)
    SphericalCoordinates(ρ, ϕ, θ)
end