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")
2 Likes
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.