Use linear interpolation between two LLA coordinate system points

General Usage
, ,
1

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
        lat::Float64
        lon::Float64
        alt::Float64
end
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) )
end

Any comments are greatly appreciated.

2

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.

5 Likes
3

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 GeographicLib.inverse or Proj4.geod_inverse
  • Interpolate along this geodesic from the first point toward the second using GeographicLib.forward or Proj4.geod_destination.

@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.

5 Likes
4

@NiclasMattsson, Is there a great circle route?

5

GMT.jl can do it.

2 Likes
6

Thanks for giving that explanation, @c42f. 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.

2 Likes
7

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 delta and 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.

1 Like
8

Here’s a solution which uses the original GeographicLib C++ library via Cxx.jl:

using Cxx
using Libdl

# Workaround weird issue with dynamic linking
Libdl.dlopen("/usr/lib/x86_64-linux-gnu/libGeographic.so", Libdl.RTLD_GLOBAL)

cxx"""
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/GeodesicLine.hpp>

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)
    lat[], lon[]
end

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))
2 Likes
9

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)

using GMT

julia> project(start_pt=(0,0), end_pt=(45,45), step=5000, km=true)
1-element Array{GMT.GMTdataset,1}:
 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)
1 Like