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

Is the following the correct way to calculate the straight line distance between two points / locations by longitude and latitude?

``````julia> using Geodesy
INFO: Precompiling module Geodesy.

julia> x_lla = LLA(-27.468937, 153.023628, 0.0)
LLA(lat=-27.468937°, lon=153.023628°, alt=0.0)

julia> y_lla = LLA(-27.465933, 153.025900, 0.0)
LLA(lat=-27.465933°, lon=153.0259°, alt=0.0)

julia> distance(x_lla, y_lla)
401.5431022017651
``````

Is there a simple / easy way to transfer the longitude and latitude into the x value and y value (cartesian coordinate system) of a point in a planar map?

A map in what projection?

If you want great circle distance, and are happy to think of the earth as a sphere, then it should just be `arccos(dot(v,w))` where v,w are 3-vectors like `v=[sin(lat), cos(lat)*sin(long), cos(lat)*cos(long)]`… up to units. But perhaps a package called Geodesy will already have this somewhere?

@improbable22, I read the package description of Geodesy and the above example if directly from their git link. But I’m not sure I totally understand it.

1 Like

Sorry perhaps a reading failure on my part.

Now that I look at the readme, it sounds like “distance” is probably a straight-line distance in 3D. Is this what you want?

They say "Future work may focus on geodesics and related calculations " which is what I was aiming for. But the package is more sophisticated than assuming the world is a sphere, and thus working out the accurate great-circle distance (i.e. geodesic length) accurately would be harder work.

To get the cartesian (2d-) distance on a map, have a look at Proj4.jl which is a Julia wrapper around a relatively widely used C library. As @improbable22 already wrote, this depends on the map projection.

1 Like

Doesn’t harversine distance in the Distances package work?

``````using Distances
julia> l1 = (-27.468937, 153.023628)
julia> l2 = (-27.465933, 153.025900)
julia> haversine(l1, l2, 6372.8)
0.39054922275889736
``````

I used it to calculate distances in an agricultural plot.

7 Likes

@alejandromerchan, yes. It’s correct. See https://www.movable-type.co.uk/scripts/latlong.html

1 Like

Just a few words to tell you that you need to be aware that straight line is not the shortest path.

No one here mentioned the difference between loxodromic route ie constant bearing and orthodromic navigation aka great-circle navigation.

I wonder if there is a Julia package which calculate bearing between 2 points (using either loxodromic navigation or orthodromic navigation).

1 Like

@scelles, there are several options for calculating great circle paths between points in Julia, including bearings:

• Straightforward Haversine calculations, as explained above—works for spheres only, but easy to implement oneself;
• Vincenty’s algorithm—works for moderately elliptical ellipsoids but has problems with near-antipodal points. A versions is available in Geodesics.jl (code mostly originally from GreatCircle.jl).
• Charles Karney’s algorithms—work for even very eccentric bodies and all combinations of points, with no slowdown compared to (2). Implemented in `GeographicLib.inverse`, or `Proj4.geod_inverse` (via the Proj4 external library).

I plan to get around to moving the GeographicLib algorithms into Geodesy.jl, when time permits.

(I don’t know of any constant-bearing methods implemented in Julia currently.)

3 Likes

For constant bearing it should only be a question of converting to Mercator (or other conformal projection) and compute the conjugate of the of the straight line angle.

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

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 )