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

Question

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)

Answer 1

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);

Answer 2

You can use a transformed projection such UTM to solve this easily. In QGIS for example. For Germany the UTM projection to use is zone 32North you can check this at: http://www.latlong.net/lat-long-utm.html

Here is the code to do it. Just start up QGIS (Free software) open the python console in QGIS (ctrl-alt-P) copy and paste the script below.

lat1=50
lon1=10
lat2=50.1
lon2=10.1
partial=1500 #distance in meters from point 1 towards point 2

from qgis import core
#Projections
crs_latlon = core.QgsCoordinateReferenceSystem(4326,core.QgsCoordinateReferenceSystem.EpsgCrsId)
crs_utm32U = core.QgsCoordinateReferenceSystem(32632,core.QgsCoordinateReferenceSystem.EpsgCrsId) 

#Transform
latlon2UTM = QgsCoordinateTransform(crs_latlon, crs_utm32U)          
UTM2latlon = QgsCoordinateTransform(crs_utm32U, crs_latlon)

#Apply transform        
UTMp1=latlon2UTM.transform(QgsPoint(lon1,lat1))
UTMp2=latlon2UTM.transform(QgsPoint(lon2,lat2))            
UTMx1=UTMp1.x()
UTMy1=UTMp1.y()
UTMx2=UTMp2.x()
UTMy2=UTMp2.y()

#linear interpolation
#Get world distance by pythagorus
world_distance=((UTMx2-UTMx1)**2+(UTMy2-UTMy1)**2)**0.5
UTMx3=UTMx1+(UTMx2-UTMx1)*partial/world_distance
UTMy3=UTMy1+(UTMy2-UTMy1)*partial/world_distance

latlonp3=UTM2latlon.transform(QgsPoint(UTMx3,UTMy3))
lon3=latlonp3.x()
lat3=latlonp3.y()

print lon3, lat3

You can prove it for your self in google earth by adding points, editing the lat and lon for each, then measure the distance from point 1 to point 3

If you dont want to use the QGIS transform functions I'm sure there are other online UTM converiosn functions you could use or you could formulate your own from the mathematical transform itself which is probably available in wikipedia.