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.

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.

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

1 Like

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)

1 Like

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.

1 Like

: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

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)