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

For reference this shows how to easily compute geodesic orthodroms … and buffers with GMT

@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) 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.

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.

# 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)

# 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")

1 Like

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.

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.

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

hmm, better show also the result

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)

1 Like

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

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.

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

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
2 Likes

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, which gives the same result as the GeographicLib loxodroms (PROJ does not have it) and the text clearly says Loxodrome on the Ellipsoid

1 Like

@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:

1 Like

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.