9

Problem:


Given two locations:

L1 = (latitude1, longitude1, timestamp1), L2 = (latitude2, longitude2, timestamp2),

and a configurable, but constant, movement speed:

v = 1.39 meters per second (for instance).

How can we interpolate between these two locations to estimate a users location as he travels from L1 to L2?


I have been searching for solutions to this problem and so far I have found, that for small distances (away from the poles) linear interpolation can be used. So, I looked up linear interpolation on Wikipedia and found this:

// Imprecise method which does not guarantee v = v1 when t = 1,
// due to floating-point arithmetic error.
float lerp(float v0, float v1, float t) {
    return v0 + t*(v1-v0);
}

So I am thinking of using this lerp function to interpolate latitude and longitude between L1 and L2. That was the easy part. How do I calculate t? I'm guessing I have to calculate some time deltas, but how do I factor in the movement speed?


Edit:

I'm testing various methods for collecting GPS-locations. For this purpose, I am recording waypoint-locations throughout a walk. I need to interpolate between these waypoints, using the movement speed, to estimate my position along the walk. Then I can compare my results to the estimates to see how well they are faring.

Example:

map

transporter_room_3
  • 2,583
  • 4
  • 33
  • 51
  • 3
    so you are saying that the devices moves from point A to point B in a fixed amount of time, but that overall their speed is something else possible totally unrelated? – njzk2 Sep 28 '15 at 18:13
  • I have implemented different strategies for collecting gps locations, e.g. min-distance between updates, min-time between updates, accelerometer-enabled strategies, etc. I need to establish a 'ground truth' to compare the results against. I do this by collecting separate waypoints throughout an experiment and I need to interpolate between these waypoints (using the movement speed), so that I can compare them to the other results. I hope this makes sense, otherwise I'll try to explain better. – transporter_room_3 Sep 28 '15 at 18:22
  • 1
    There is a problem in your question: you say you want a constant configurable speed v. However you have measured timestamps and locations, which implicitly represent a speed. If the person was running at 10 km/h you cannot interpolate that with a given wrong speed. I dont see why you need the speed. you have the time stamps, as reference for your interpolation. You cannot interpolate at both. Decide what is your reference; Either it it is the pair (time and location) ; or it is the pair (speed and location) which is more complex. – AlexWien Sep 28 '15 at 19:56
  • Good point, AlexWien. I just realized this myself. I've got to figure out which approach to go with. I think the former is most relevant to my needs. – transporter_room_3 Sep 28 '15 at 20:53
  • The former approach is seen in my answer. Do you do this task for yourself, (app), university, or commercially? (get paid for) – AlexWien Sep 28 '15 at 21:32
  • I am developing this app for a course in university. The waypoint interpolation is a small but important part in the app. I appreciate the help you've provided. I implemented my own flavor of what you posted and it works like a charm. Thanks! – transporter_room_3 Sep 28 '15 at 21:50

5 Answers5

18

Take a close look at Calculate distance, bearing and more between Latitude/Longitude points

It contains several formulas and JavaScript examples that might help you. I know that it's NOT Java, but it should be simple enough to port the code over. Especially given the detailed description of the formula.

EDIT:

While it seems OK to use linear interpolation for shorter distances, it can in fact be quite off, especially as you get closer to the poles. Seeing from the example that you are in Hamburg, this will already have an effect that's noticable over a few hundred meters. See this answer for a good explanation.

The Problem: The distance between 1 degree in longitude varies greatly depending on your latitude.

This is because the earth is NOT flat, but a sphere - actually an ellipsoid. Therefore a straight line on a two dimensional map is NOT a straight line on the globe - and vice versa.

To get around this problem one can use the following approach:

  1. Get the bearing from the start coordinate (L1) to the end coordinate (L2)
  2. Calculate a new coordinate from the start coordinate (L1) along a great circle path, given the calculated bearing and a specified distance
  3. Repeat this process, but using the newly calculated coordinate as the starting coordinate

We can create a few simple functions that will do the trick for us:

double radius = 6371; // earth's mean radius in km

// Helper function to convert degrees to radians
double DegToRad(double deg) {
    return (deg * Math.PI / 180);
}

// Helper function to convert radians to degrees
double RadToDeg(double rad) {
    return (rad * 180 / Math.PI);
}

// Calculate the (initial) bearing between two points, in degrees
double CalculateBearing(Location startPoint, Location endPoint) {
    double lat1 = DegToRad(startPoint.latitude);
    double lat2 = DegToRad(endPoint.latitude);
    double deltaLon = DegToRad(endPoint.longitude - startPoint.longitude);

    double y = Math.sin(deltaLon) * Math.cos(lat2);
    double x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(deltaLon);
    double bearing = Math.atan2(y, x);

    // since atan2 returns a value between -180 and +180, we need to convert it to 0 - 360 degrees
    return (RadToDeg(bearing) + 360) % 360;
}

// Calculate the destination point from given point having travelled the given distance (in km), on the given initial bearing (bearing may vary before destination is reached)
Location CalculateDestinationLocation(Location point, double bearing, double distance) {

    distance = distance / radius; // convert to angular distance in radians
    bearing = DegToRad(bearing); // convert bearing in degrees to radians

    double lat1 = DegToRad(point.latitude);
    double lon1 = DegToRad(point.longitude);

    double lat2 = Math.asin(Math.sin(lat1) * Math.cos(distance) + Math.cos(lat1) * Math.sin(distance) * Math.cos(bearing));
    double lon2 = lon1 + Math.atan2(Math.sin(bearing) * Math.sin(distance) * Math.cos(lat1), Math.cos(distance) - Math.sin(lat1) * Math.sin(lat2));
    lon2 = (lon2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalize to -180 - + 180 degrees

    return new Location(RadToDeg(lat2), RadToDeg(lon2));
}

// Calculate the distance between two points in km
double CalculateDistanceBetweenLocations(Location startPoint, Location endPoint) {

    double lat1 = DegToRad(startPoint.latitude);
    double lon1 = DegToRad(startPoint.longitude);

    double lat2 = DegToRad(endPoint.latitude);
    double lon2 = DegToRad(endPoint.longitude);

    double deltaLat = lat2 - lat1;
    double deltaLon = lon2 - lon1;

    double a = Math.sin(deltaLat / 2) * Math.sin(deltaLat / 2) + Math.cos(lat1) * Math.cos(lat2) * Math.sin(deltaLon / 2) * Math.sin(deltaLon / 2);
    double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

    return (radius * c);
}

This uses a mean earth radius of 6371 km. See Wikipedia for an explanation of this number and its accuracy.

One can now calculate a new intermediary location between the two points, given a distance travelled (in km):

double bearing = CalculateBearing(startLocation, endLocation);

Location intermediaryLocation = CalculateDestinationLocation(startLocation, bearing, distanceTravelled);

Assuming a speed of v (e.g. 1.39) meters per second, one can now use a simple for loop to get points 1 second apart:

List<Location> locations = new ArrayList<Location>();

// assuming duration in full seconds
for (int i = 0; i < duration; i++){
    double bearing = CalculateBearing(startLocation, endLocation);
    double distanceInKm = v / 1000;
    Location intermediaryLocation = CalculateDestinationLocation(startLocation, bearing, distanceInKm);

    // add intermediary location to list
    locations.add(intermediaryLocation);

    // set intermediary location as new starting location
    startLocation = intermediaryLocation;
}

As an added bonus, you can even determin the time required to travel between any two points:

double distanceBetweenPoints = CalculateDistanceBetweenLocations(startPoint, endPoint) * 1000; // multiply by 1000 to get meters instead of km

double timeRequired = distanceBetweenPoints / v;

This will result in greater accuracy over any distance than a simple linear interpolation using just the delta of the coordinates. Although this approach is not perfect, it will have an error of generally 0.3% or less, which is quite acceptable. If you need a better solution, you might want to look into the Vincenty formula.

davethebrave
  • 823
  • 9
  • 21
  • I don't see anything which answers his question. – AlexWien Sep 28 '15 at 18:42
  • No, just posting a link without understanding the content make not any sense. The formular "distance" on that page does what it is named for: It gives the distance in meters between two coordinates. He does not need that. He must interpolate between the two measures. – AlexWien Sep 28 '15 at 19:06
  • You're right. That alone doesn't do it. However, using the information on that page one can come up with the correct answer and with good accuracy. Unfortunately calculations with latitude and longitude require a bit more math than simply adding and subtracting the points. – davethebrave Sep 28 '15 at 20:07
  • 1
    Link-only answers are bad answers: http://meta.stackexchange.com/questions/8231/are-answers-that-just-contain-links-elsewhere-really-good-answers – Juan Lopes Sep 29 '15 at 01:34
  • 2
    This was a very helpful answer! My only note would be in your method `CalculateDistanceBetweenLocations` you have `double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - 1));` This should instead be: `double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));` – user3798780 Apr 16 '17 at 04:19
  • 1
    I have submitted an edit to correct the CalculateDistanceBetweenLocations error reported by user3798780. The original posting should now be correctly displaying `Math.sqrt(1 - a)` – hsoi Aug 23 '19 at 21:29
  • 1
    Great answer, but can the OP be edited to spell `longitude` correctly? Line 7 in the `CalculateDestinationLocation` function spells it `logintude`. – lolololol ol Nov 17 '21 at 17:45
  • Thanks @lololololol. I didn't notice this before. I fixed the spelling. – davethebrave Nov 21 '21 at 10:00
6

Calculations like these are actually very simple if you first convert your lat/longs to n-vectors (https://en.wikipedia.org/wiki/N-vector). After converting you can use standard interpolation, and you will also avoid any problems with long distances, the poles or the dateline.

If you check "External links" on the Wikipedia page, there is a page (http://www.navlab.net/nvector/) where ten problems are solved, and problem 6 on that page (interpolated position) should be the same as your question. As you can see, that solution is exact for any distance and also works at any earth-positions, like the poles.

MathGuy
  • 61
  • 1
  • 4
3

I'm guessing I have to calculate some time deltas, but how do I factor in the movement speed?

At linear interpolation, in your casem you iterate between two time points, using the iteration variable t which runs from start time t1 to end time t2, with a predefined step. Asume step = 1 second, which is quite usable for your application.

long t1 = location1.getTimeStamp(); // in milliseconds;
long t2 = location2.getTimeStamp();
double deltaLat = location2.latitude - location1.latitude;
doule deltaLon =  location2.longitude- location1.longtude;
// remove this line if you don't have measured speed:
double deltaSpeed =  location2.speed - location1.speed;

long step = 1 * 1000; // 1 second in millis 
for (long t = t1; t1 < t2; t+= step) {

   // t0_1 shall run from 0.0 to (nearly) 1.0 in that loop
  double t0_1 = (t - t1) / (t2 - t1);
  double latInter = lat1 + deltaLat  * t0_1;
  double lonInter = lon1 + deltaLon  * t0_1;
  // remove the line below if you dont have speed
  double speedInter = speed1 + deltaSpeed  * t0_1;
  Location interPolLocation = new Location(latInter, lonInter, speedInter);
  // add interPolLocation to list or plot.
}
AlexWien
  • 28,470
  • 6
  • 53
  • 83
  • 2
    This code assumes that the latitude and longitude are on a flat Cartesian plane. Earth isn't exactly that, Earth is a sphere. – Ali Ashraf Jan 27 '17 at 14:48
2

There are some other interpolation strategies that perform better than linear interpolation, includind kinematic interpolation, which takes as input the initial and final speed of anchor points. As an example, see this comparison from a recent paper (Long JA (2015) Kinematic interpolation of movement data. Int J Geogr Inf Sci 8816:1–15. doi: 10.1080/13658816.2015.1081909):

From Long JA (2015) Kinematic interpolation of movement data. Int J Geogr Inf Sci 8816:1–15. doi: 10.1080/13658816.2015.1081909

There are R and Python implementations for the Kinematic Interpolation. It should be easy to write a Java version.

1

Yep. Linear interpolation.

L1 = (1, 2, 3)
L2 = (4, 5, 6)
desired_number_of_interpolation_points = 9
interpolation_points = []

lat_step = (L2[0] - L1[0]) / (desired_number_of_interpolation_points + 1)
lon_step = (L2[1] - L1[1]) / (desired_number_of_interpolation_points + 1)
time_step = (L2[2] - L1[2]) / (desired_number_of_interpolation_points + 1) 

for i in range(1, desired_number_of_interpolation_points + 1)
    interpolation_points.append((lat_step * i, lon_step * i, time_step * i))
Chad Kennedy
  • 1,716
  • 13
  • 16