I have the following simple piece of code to find the position between two coordinate points using linear interpolation. I’m not sure if this is a good / feasible solution or not. So what do you think?
mutable struct LLA
LLA(lat,lon) = LLA(lat,lon,0.0)
function linearInterpolation ( startLoc::LLA, endLoc::LLA, fraction::Float )
startLat = startLoc.lat
startLon = startLoc.lon
endLat = endLoc.lat
endLon = endLoc.lon
return LLA( startLat + fraction * (endLat - startLat), startLon + fraction * (endLon - startLon) )
Any comments are greatly appreciated.
Linear interpolation of the coordinates will only be approximately correct for points that are near each other, and it fails completely if you cross the date line. For larger distances you need to account for the curvature of the earth by calculating the great circle distance - and for that you need to use spherical trigonometry (instead of ordinary planar). The standard approach for lat-lon coordinates is to use the haversine formula, or alternatively n-vectors. Then I guess you could just interpolate the two altitudes if that makes sense for your application.
The most accurate way to do this would be to compute a geodesic line via the inverse geodesic problem, and use this object to interpolate between the lat-lon endpoints. This is the method used in the GeographicLib waypoints example in python (see also the underlying C++ documentation).
I don’t think we have a registered julia package which can do this in the most natural way yet. The obvious candidates would be Proj4.jl (wraps the C proj.4 library) and Geodesy.jl (pure julia, with some algorithms ported from GeographicLib). This may eventually be in Geodesy.jl (see https://github.com/JuliaGeo/Geodesy.jl/issues/40) but in the meantime you could try out https://github.com/anowacki/GeographicLib.jl or
Proj4.jl with a similar two-step procedure:
- Compute the azimuth and great circle distance along the geodesic from the first to second point using
- Interpolate along this geodesic from the first point toward the second using
@visr @anowacki is there a better alternative?
All the above is for interpolating the lat-lon components. For the altitude you can probably just use linear interpolation based on the fractional distance along the geodesic.
@NiclasMattsson, Is there a great circle route?
Thanks for giving that explanation, @Chris_Foster. That is indeed how I would do this. GeographicLib.jl now also exposes the waypoint functionality of the original library with
waypoints, though this doesn’t include fractional distances. (Perhaps an additional option could be added.)
As Chris says, this should be added to Geodesy if I can get round to it at some point this coming semester.
The Haversine method suggested by Niclas is for great circles on the sphere, so if the accuracy you need is not high, or you are only interested in spheres, that will do the job fine. (Feel free to steal the implementation of
step in https://github.com/anowacki/assorted-julia-modules/blob/master/SphericalGeom.jl). If accuracy on the real Earth is required, the other routes explained above with Proj or GeographicLib are needed.
Here’s a solution which uses the original
GeographicLib C++ library via Cxx.jl:
# Workaround weird issue with dynamic linking
using namespace GeographicLib;
# Get `(lat,lon)` position `d` meters along a GeodesicLine
function geoline_position(line, d)
lat = Ref(0.0)
lon = Ref(0.0)
@cxx line->Position(d, lat, lon)
latlon1 = (40.1, 116.6)
latlon2 = (37.6, -122.4)
geod = icxx"Geodesic(Constants::WGS84_a(), Constants::WGS84_f());"
line = @cxx geod->InverseLine(latlon1..., latlon2...)
dist = @cxx line->Distance()
numpoints = ceil(Int, dist/1000e3)
positions = geoline_position.(Ref(line), LinRange(0, dist, numpoints))
I was not sure of this but checked, GMT only does this on the sphere. If that’s enough, with the master version, to interpolate at 5000 km steps along a line from long, lat (0,0) to (45,45)
julia> project(start_pt=(0,0), end_pt=(45,45), step=5000, km=true)
GMT.GMTdataset([-3.1805546814635168e-15 3.1805546814635168e-15 0.0; 29.970589500178583 35.24036543346766 5000.0; 44.99999999999997 44.999999999999986 6671.703118513764], String, "", String, "", "")
with GMT.jl 0.11 version the syntax will have to be a bit different
julia> project(start_pt=(0,0), end_pt=(45,45), G=5000, Q=true)