# Geodesy: how to calculate the straight line distance between two locations, which is represented by longitude and latitude?

**URL:** <https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984>\
**Category:** Geo\
**Tags:** diffeq, gmt, geodesy\
**Created:** [January 23, 2019, 2:03pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984 "2019-01-23T14:03:11Z")\
**Posts on this page:** 10\
**Page:** 2

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 19, 2021, 11:12pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/21 "2021-06-19T23:12:18Z")

</div>

For reference [this shows](https://www.generic-mapping-tools.org/GMT.jl/dev/gallery/buffergeo/buffer/) how to easily compute geodesic orthodroms … and buffers with GMT

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 20, 2021, 6:14pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/22 "2021-06-20T18:14:04Z")

</div>

> [@lazarusA](#):
>
> do you know how to get the geodesic line data (x,y,z)

@lazarus, please see code below that computes the geodesic line over the GRS80 ellipsoid given the (LAT, LONG) of 2 initial points.

This is a quick&dirty attempt to implement the differential equations in [Panou et al. (2013)](https://www.researchgate.net/publication/258797115_Solving_the_geodesics_on_the_ellipsoid_as_a_boundary_value_problem) with some simplifications. In particular, the iterative process to fine tune one of the boundary conditions was simplified, hopefully not too much.

In any case, the results seem to match the authors’ numeric tables for different points tested, in terms of geodesic length and Clairaut constants. The results seem to match also those from [this web calculator](https://geographiclib.sourceforge.io/cgi-bin/GeodSolve).

Not being an expert in any of these things, much less in Julia, any advice on how to improve the differential equations accuracy and the code are very welcome.

```julia
# Reference:
# Panou, G., Delikaraoglou, D. and Korakitis, R. (2013)
# Solving the geodesics on the ellipsoid as a boundary value problem.
# Journal of Geodetic Science, 3(1)

using DifferentialEquations

function geodesic!(du,u,p,λ)
    du[1] = u[2]
    du[2] = h1(u[1])*u[2]^2 + h2(u[1])
    du[3] = u[4]
    du[4] = (h3(u[1])*u[2]^2 + h4(u[1]))*u[3] + 2*h1(u[1])*u[2]*u[4]
end

trapez(t,A) = @views sum((A[1:end-1] + A[2:end]).*diff(t)/2) # trapezoidal integration

# GRS80 ellipsoid (Moritz, 1980)
a = 6378137; b = 6356752.3141 # m
e = sqrt(a^2 - b^2)/a; e2 = e^2

N(φ) = a / sqrt(1 - e2*sin(φ)^2) # radius of curvature in prime vertical normal section
E(φ) = a^2*(1 - e2)^2 / (1 - e2*sin(φ)^2)^3 # first fundamental coefficients E, F=0, G
G(φ) = a^2*cos(φ)^2 / (1 - e2*sin(φ)^2)

h1(φ) = -(2*(1 - e2)*tan(φ)+ 3*e2*sin(φ)*cos(φ)) / (1 - e2*sin(φ)^2)
h2(φ) = -(1 - e2*sin(φ)^2)*sin(φ)*cos(φ) / (1 - e2)
h3(φ) = -2*(1 - e2)*(1 - e2*sin(φ)^2 + 2*e2*sin(φ)^2 * cos(φ)^2) / ((1 - e2*sin(φ)^2)^2 * cos(φ)^2) +
      - 3*e2*((1 - e2*sin(φ)^2)*(cos(φ)^2 - sin(φ)^2) + 2*e2*sin(φ)^2 * cos(φ)^2) / (1 - e2*sin(φ)^2)^2
h4(φ) = -((1 - e2*sin(φ)^2)*(cos(φ)^2 - sin(φ)^2) - 2*e2*sin(φ)^2 * cos(φ)^2) / (1 - e2)

# INPUT geodetic endpoints (LAT, LONG), convert from degrees to radians
# NOTE: geodetic LAT: –π/2 < φ < π/2; geodetic LONG: –π < λ < π 
d2r = π/180
φ0 = 60*d2r ; λ0 = 0*d2r # λ0 < λ1
φ1 = -30*d2r ; λ1 = 120*d2r # λ0 = λ1 (meridian) excluded as trival case

# initial conditions
u0 = zeros(4)
u0[1] = φ0
u0[2] = (φ1 - φ0)/(λ1 - λ0) # need to iterate to find this initial condition so that geodetic ends at φ1(λ1)
u0[3] = 0.
u0[4] = 1.
λspan = (λ0, λ1) # longitude-span for ODE solving (takes the role of time variable)

prob = ODEProblem(geodesic!, u0, λspan) # define problem
sol = solve(prob, alg=RK4())
φ10 = sol[1,end]

n = 8000
λ = LinRange(λ0,λ1,n)
dϕ = 0.05*u0[2] # perturb dφ/dλ by 5% (empirical, can be improved)
edφ = +Inf
while abs(edφ) > 1e-12
    u0[2] += dϕ
    prob = ODEProblem(geodesic!, u0, λspan) # define problem
    sol = solve(prob, alg=RK4(),saveat=λ,reltol=1e-12)
    φ11 = sol[1,end]
    edφ = φ1 - φ11
    dφ1dφ0 = (φ11 - φ10)/dϕ # simplification of Panou et al., 2013, Eq.33)
    dϕ = edφ/dφ1dφ0
    φ10 = φ11
end
λ = sol.t
φ = sol[1,:]
dφ = sol[2,:]

# geodetic distance : s
A = @. sqrt(E(φ)*(dφ)^2 + G(φ))
s = trapez(λ, A) # geodetic distance

# Clairaut constant : C
C = mean(@. G(φ)/sqrt(E(φ)*(dφ)^2 + G(φ)))

# Plot computed geodesic line on ellipsoid using GMT
using Printf, GMT
str = @sprintf("%i m, over (LAT,LONG) = (%i,%i) -> (%i,%i) deg",s,φ0/d2r,λ0/d2r,φ1/d2r,λ1/d2r);
coast(region=:global, proj=(name=:ortho, center=((λ0+λ1)/d2r/2,(φ0+φ1)/d2r/2)), land=:lightblue, frame=:g,
      title="Geodesic L=$str", par=(FONT_TITLE="14p,Helvetica,black",))
GMT.plot!([λ/d2r φ/d2r], lw=1.5, lc=:red, show=true)

```

 ![GRS1980_Geodesic_DiffEqs_with_GMT_plot](https://global.discourse-cdn.com/julialang/original/3X/f/0/f0692c2d9b099f46219eaf9de7698df4582a32a9.png)

```julia
# plot computed geodesic line on ellipsoid using Plots.jl pyplot()
using Printf, Plots; pyplot()
n = 100; u = collect(LinRange(0,2*π,n)); v = collect(LinRange(0,π,n));
x = a*cos.(u)*sin.(v)'; y = a*sin.(u)*sin.(v)'; z = b*ones(n)*cos.(v)';
Plots.surface(x,y,z, colorbar=false)
x0 = @. N(φ)*cos(φ)*cos(λ);
y0 = @. N(φ)*cos(φ)*sin(λ);
z0 = @. N(φ)*(1 - e2)*sin(φ);
str = @sprintf("%i m, over (LAT,LONG) = (%i,%i) -> (%i,%i) deg",s,φ0/d2r,λ0/d2r,φ1/d2r,λ1/d2r);
plot!(x0,y0,z0, lc=:red, label="geodesic L=$str")

```

 ![GRS1980_Geodesic_DiffEqs_with_Plots_pyplot](https://global.discourse-cdn.com/julialang/original/3X/f/c/fc161ba34c7b4d8b12843b68330eb6b3da3aeff2.png)

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [June 20, 2021, 6:31pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/23 "2021-06-20T18:31:50Z")

</div>

thanks, I will tested it. I did find another solution for this case using Geodesy.jl I will compare and report back. But it was definitely much shorter.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 20, 2021, 7:20pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/24 "2021-06-20T19:20:21Z")

</div>

Like shown in the examples linked in my post above, you now do that in GMT with (for example). It uses PROJ for the calculations.

```julia
ortho = orthodrome([0 0; 70 60], step=1000, unit=:k);

```

hmm, better show also the result

```julia
using GMT
# Compute the great circle (orthodrome) line
ortho = orthodrome([0 0; 70 60], step=1000, unit=:k);
# Plot the orthodrome on an orthographic projection
coast(region=:global, proj=(name=:ortho, center=(0,45)), land=:peru, frame=:g)
plot!(ortho, lw=0.5, marker=:circ, ms=0.1, fill=:black, show=true)

```

 ![buffer1](https://global.discourse-cdn.com/julialang/original/3X/4/4/44420af82b5cdd4c5a552b579420d271b1c648b5.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 20, 2021, 7:35pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/25 "2021-06-20T19:35:10Z")

</div>

@joa-quim, an “orthodrome” is not the same thing as a geodesic on the ellipsoid, is it?

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [June 20, 2021, 7:49pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/26 "2021-06-20T19:49:30Z")

</div>

the problem with that approach, at least for me is that I usually like to have the data first and afterwards plotted. Data comes first. I don’t see how to get (x,y,z) from the above. Plus, here it seems like the line is effectively just the x,y projected.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 20, 2021, 7:53pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/27 "2021-06-20T19:53:50Z")

</div>

But you get the data (lon, lat positions on the selected ellipsoid)

```julia
ortho = orthodrome([0 0; 70 60], step=1000, unit=:k)
10×2 Matrix{Float64}:
  0.0 0.0
  4.32514 7.93354
  8.81972 15.8175
 13.676 23.5953
 19.1363 31.1948
 25.5301 38.5144
 33.3211 45.3992
 43.1465 51.6037
 55.7573 56.738
 70.0 60.0

```

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 20, 2021, 7:58pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/28 "2021-06-20T19:58:04Z")

</div>

If the computation is on the ellipsoid, as this one is, yes I think it is. Or maybe I’m wrong and gave it a misleading name. But this one implements the PROJ’s `geod` program and can work on the sphere or the ellipsoid.

I also implemented an `loxodrome` after [this reference](http://www.mygeodesy.id.au/documents/Loxodrome%20on%20Ellipsoid.pdf), which gives the same result as the GeographicLib loxodroms (PROJ does not have it) and the text clearly says `Loxodrome on the Ellipsoid`

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 20, 2021, 8:56pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/29 "2021-06-20T20:56:36Z")

</div>

@joa-quim, plotted your GMT “orthodrome” points over the geodesic solution computed via differential equations (code posted above) and for your info, they do match:

 ![Geodesic_and_GMT_orthodrome](https://global.discourse-cdn.com/julialang/original/3X/c/1/c16bd5e5d8525510f4a5c32265a01ab7284acccd.png)

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 20, 2021, 10:00pm UTC](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984/30 "2021-06-20T22:00:16Z")

</div>

Boa bravo Rafael, I had quite some work to implement this, and some companion functions that namely compute geodesic buffers which I think is a near unique feature that only ArcGIS has, but the real hard work here is done by PROJ & GDAL. But yours came out of your _engenho_.

[Previous page](https://discourse.julialang.org/t/geodesy-how-to-calculate-the-straight-line-distance-between-two-locations-which-is-represented-by-longitude-and-latitude/19984.md?page=1)