# 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

a = (sin(dlat / 2) * sin(dlat / 2) +
sin(dlon / 2) * sin(dlon / 2))
c = 2 * atan(sqrt(a), sqrt(1 - a))
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
end