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

compute the conjugate of the of the straight line angle

@joa-quim I’m not sure to understand what you mean about “conjugate”

In Mercator projection, bearing should simply be angle between segment North-Point1 and Point1-Point2.

Is there some Julia packages to show earth (for example) and see how trajectory between two points can be displayed on it and using different map projections?

PS : something like http://ressources.univ-lemans.fr/AccesLibre/UM/Pedago/physique/02/divers/ortholoxo.html but with Julia

Well, the angle between 2 points in a cartesian axis (theta) is the angle counterclockwise from the x-axis, but bearing (azimuth) is the angle clockwise from North. So azimuth = 90 - theta

You can do all sorts of mapping with GMT.jl See for example this example and the corresponding julia script.

And as a simpler example, let’s plot an orthodrome between two points

ortho = project(start_pt=[-50.0 10], end_pt=[50 60], step=10, km=true);
coast(region=[-180 180 -90 90], proj=:equidistCylindrical, frame=:g, res=:crude, land=:navy, B=:a)
plot!(ortho, lw=3, lc=:green, name="ortho_line.png", show=true)

(to plot an loxodrome one could interpolate linearly between two points in a Mercator, convert back to gepgraphic (using mapproject) and add that line to the plot below)

2 Likes

If you are only interested in the visualisation, this can be shortened to

coast(region=[-180 180 -90 90], proj=:equidistCylindrical, frame=:g, res=:crude, land=:navy, B=:a)
plot!([-50 10; 50 60], lw=3, lc=:green, name="ortho_line.png", show=true)

since GMT plots great circles between points by default.

2 Likes

:blush: Off course. I tried to see if there was an option to generate loxodromes that I forgot the obvious.

Thanks for your examples with GMT.jl
Being able to also plot loxodrome will be great!
I hope it will be possible to have such plots with GeoMakie.jl also. Pinging @asinghvi17

1 Like

It should be - if you have the geodesic line, all you need to do is get the coordinates in Mercator space and construct a LinRange for lat and lon - then, you can reproject into any coordinate system.

I haven’t been able to create an example yet, but probably will be able to in the next couple of days.

1 Like

Turns out that GMT can compute loxodromes directly with the module sample1d (sorry the syntax for this case was not yet expanded to the more comprehensible keywords)

loxo = sample1d([-50 10; 50 60], resamp="R+l", equi_space="10k");
ortho = project(start_pt=[-50.0 10], end_pt=[50 60], step=10, km=true);
coast(region=:global, land=:navy, res=:crude)
plot!(ortho, lw=2, lc=:green)
plot!(loxo, lw=2, lc=:red, name="loxo_ortho.png", show=true)

2 Likes

Hey, do you know how to get the geodesic line data (x,y,z)[lon,lat,alt] between two points coordinates?, I would like to have the data and then plot it with whatever library I wish, as you mentioned before.

Here is a translation to Julia of a Matlab code I wrote years ago. It is approximate in that it assumes a spherical earth model, but I found it completely adequate for plotting purposes.

using LinearAlgebra: ⋅, norm, ×

"""
    gcinterp(p1, p2, dθ)

Compute list of points on a great circle interpolation assuming a spherical earth model.

### Input Arguments

- `p1`, `p2`:  Iterables of length 2 containing the longitude λ and lattitude L 
               in degrees of the initial and terminal points of the desired arc.

- `dθ`: The desired angular increment in degrees for point spacing along 
        the interpolated great circle arc.

### Return Value
A vector of ordered pairs `(λᵢ, Lᵢ)` (longitudes and lattitudes in degrees) of
points on the great circle path joining `p1` and `p2`.  The arc length spacing of the
points will be equal to or slightly less than `dθ`.
"""
function gcinterp(p1, p2, dθ)
    λ₁, L₁  = float.(p1[1:2])
    λ₂, L₂  = float.(p2[1:2])
    rg1 = ll2rect(λ₁, L₁)
    rg2 = ll2rect(λ₂, L₂)
    ct0 = rg1 ⋅ rg2
    st0 = norm(rg1 × rg2)
    rg1 /= st0
    rg2 /= st0
    θ0 = atand(st0,ct0) # great circle arc length
    n = 1+ceil(Int, θ0/dθ)
    points = [(λ₁,L₁) for _ in 1:n]
    for (i,θ) in enumerate(range(0, θ0, length=n))
        rvec = rg1 * sind(θ0-θ) + rg2 * sind(θ)
        points[i] = rect2ll(rvec)
    end
    points[1] = p1[1], p1[2] # Eliminate rounding errors
    points[n] = p2[1], p2[2]
    points
end

"""
    ll2rect(λ, L)

Convert longitude `λ` and lattitude `L` (both in degreees) to a 3-vector of rectangular coordinates.
"""
function ll2rect(λ, L)
    sinλ, cosλ = sincosd(λ)
    sinL, cosL = sincosd(L)
    [cosL * cosλ, cosL * sinλ, sinL]
end

"""
    rect2ll(rvec)

Convert rectangular coordinates (a 3-vector of unit length) to a (longitude,lattitude) tuple.
Angular units on output are degrees.
"""
function rect2ll(rvec)
    λ = atand(rvec[2], rvec[1])
    L = asind(rvec[3])
    (λ, L)
end

Hope this helps.

2 Likes

Thanks, this works just fine for plotting. (Assuming we live in an spherical earth :smiley: )

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.