[GIS] How to interpolate a lat/long point between two others at a given distance

distanceinterpolationpoint

I'm trying to solve the problem of interpolating a point at a given distance (in meters) between two other lat/long points. (e.g. P1 (50,10), P2(50.01, 9.990), I want to interploate P3 at a distance "X" on the shortest route between P1 and P2, assuming that X is always smaller than the distance between P1 and P2)

Whilst at first this seemed kind of easy I couldn't for the life of me find an easy to implementation for this (e.g. a few lines of code and a few external constants, as to be easy to debug and quick to run).

Currently I'm using the formula:

val x : Double = RadiusOfEarth*Math.toRadians(lon)*CentralLatCos
val y : Double = RadiusOfEarth*Math.toRadians(lat)

to project to 2d cartesian, I interpolate the point in cartesian coordinates (since for that there are loads of formulas that I actually understand and can use). And then I go back to lon/lat using:

val lat : Double = Math.toDegrees(y/RadiusOfEarth)
val lon : Double = Math.toDegrees( x/(RadiusOfEarth*CentralLatCos) )

However, this seems to be causing quite large errors so I'm quite sure the way I project my points is flawed.

Could you help me by suggesting either and easy algorithm to project back/forth from sphere to cartesian (that keep heading&distance) or a way to interpolate a point between two other lon/lat points at a given distance ?

Also, please do provide some source material that "proves" you solution so that I can better understand it, I asked this question somewhere else and some formulas with no context where thrown at me (which doesn't really help, I can find loads of those using the internet).

Do note, I'm trying to do this for small distance, the largest of which might be something around the lines of 2-3k meters (distance between the tow original points I know). So I really have no need for a formula that's accurate on a global scale, if it has a small error but its easy to implement I'm more than fine with it.

(Note: I'm working on a map between lat 45 – 55 and lon 7-13)

Best Answer

The name of the curve you are looking to interpolate along is called either a great circle for a sphere, or a geodesic for an earth-shaped ellipsoid of revolution.

Using GeographicLib, you can create an inverse geodesic between P1 and P2, and interpolate points in between, using distances in metres. There are bindings to the most popular programming languages, including Java. You can follow this example to see how a InverseLine object is created and used interpolate points.

GeodesicLine line = geod.InverseLine(50.0, 10.0, 50.01, 9.990,
                                     GeodesicMask.DISTANCE_IN |
                                     GeodesicMask.LATITUDE |
                                     GeodesicMask.LONGITUDE);

// Here is the half-way point in the middle of the arc
GeodesicData g = line.Position(line.Distance() / 2.0,
                               GeodesicMask.LATITUDE |
                               GeodesicMask.LONGITUDE);
System.out.println("Mid-point: " + g.lat2 + " " + g.lon2);
Related Question