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
```

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
```

© Daniel Marvin. Last modified: July 19, 2021.