8

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)

3
  • Welcome to GIS Stack Exchange! Do you know the path between P1 and P2? I mean, are you managing some raster/shapefile or you are talking about a theoretical (math) problem?
    – mgri
    Dec 18, 2016 at 15:39
  • Well, no, that's one of my problems. I'm using something called graphhoper which simply gives me lists of points and not the shape of the route between them.
    – George
    Dec 18, 2016 at 15:45
  • As far as I understood though the points/nodes are joined by the shortest possible path (i.e. what can be thought of as "straight" lines since the points are very close)
    – George
    Dec 18, 2016 at 15:46

2 Answers 2

5

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

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.

3
  • While I do appreciate your answer I can't exactly inject a huge python script into my current works (done in java/Scala), hence why i asked for a formula as opposed to a library or software or web-service.
    – George
    Dec 18, 2016 at 22:20
  • The formula for UTM doesn't seem to be particularly easy to implement from scratch and I have a hard time finding a Java lib with examples and documentation :/
    – George
    Dec 18, 2016 at 22:21
  • OK how about this... which I found as the first answer to searching for "convert UTM to geographic JAVA". http://stackoverflow.com/questions/176137/java-convert-lat-lon-to-utm
    – Mr Purple
    Dec 18, 2016 at 22:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.