8

I'm working on implementing the OSM Mobile Binary format for map data, and because it uses a 16 bit delta for each point along a way, it cannot represent a way if there is a large distance between two points.

The documented solution is to inject new points into the way, to keep the gap smaller.

Here's the relevant section of my code:

// calculate the distance between this point and the previous point,
// multiplied by 1,000,000 so it can be represented as an integer.
int32_t realNodeLatDelta = ((point.latitude - lastPoint.latitude) * 1000000);
int32_t realNodeLonDelta = ((point.longitude - lastPoint.longitude) * 1000000);

// cast as a 16 bit int (reduces filesize by about half)
int16_t nodeLatDelta = realNodeLatDiff;
int16_t nodeLonDelta = realNodeLonDiff;

// check if the cast caused corruption
if (realNodeLatDelta != nodeLatDelta || realNodeLatDelta != nodeLatDelta) {
  // if this point is reached, we need to add new points in the way,
  // which can be represented in a 16 bit delta
}

I'm not sure how to calculate the position of new points along the way? I assume it needs to be a great circle?

Points far enough apart to cause this issue are rare, but when they do happen it's usually a political boundary or similar, which tend to be quite long, and need to be accurately represented.

6
  • I suggest asking this question on OSM-dev mailing list (lists.openstreetmap.org/listinfo/dev), it's a bit too specific for a general GIS forum.
    – Igor Brejc
    Jan 8, 2012 at 8:28
  • 1
    Surely the math for finding a point half way between two other points isn't very specific? Jan 8, 2012 at 8:49
  • It isn't, but maybe there are other ways of dealing with this issue, which is format-specific.
    – Igor Brejc
    Jan 8, 2012 at 11:53
  • There aren't any alternative methods to deal with it. A 16 bit integer is the only format allowed for points in a way, and (at 6 decimal places of accuracy) a 16 bit integer is only capable of representing a difference of about 3.5km between two points. The only alternative is not to use that format (which is my interim solution). But I want to use something standard if possible. Jan 8, 2012 at 16:54
  • 1
    I think you may be off by an order of magnitude: the longest distances you can find on the earth are about 20,000 km. With 16 bits, that translates to 300m precision.
    – whuber
    Jan 12, 2012 at 0:11

3 Answers 3

15

The best accuracy is obtained with ellipsoidal models. In the interests of simplicity you want to avoid those when you have to code distances yourself. We pay a price: given that the earth's flattening is about 1/300, using a purely spherical model can potentially introduce a relative distance error of up to 1/300 for very long routes: about 3000 parts per million. This is worth exploring.

First, the spherical formula: simply convert (lat, lon) to cartesian coordinates, average the two cartesian points, then convert the average back to spherical coordinates. Here is pseudocode:

function cartesian(f,l) { // f = latitude, l = longitude
    return (cos(f)*cos(l), cos(f)*sin(l), sin(f))
}
function spherical(x,y,z) {
    r = sqrt(x^2 + y^2)
    if (r == 0) {
        if (z > 0) return (90, 0) 
        elseif (z < 0) return (-90, 0)
        else return Undefined // (x,y,z) == (0,0,0)
    } else {
        return (atan2(r, z), atan2(x, y)) // atan2 must return *degrees*
    }
}
function midpoint(f0,l0, f1,l1) { 
    return spherical((cartesian(f0,l0) + cartesian(f1,l1))/2)
}

(The arithmetic for midpoint involves a vector sum and a scalar division of that sum, so it's really hiding three sums and three divisions.)

That's the midpoint in spherical geometry. The calculation requires two cosines, two sines, a square root, two arctangents, and some multiplication and addition: reasonably fast and easy. There will be no problems near the poles or crossing the +-180 meridian. The result will be undefined when the two points are diametrically opposite.

One way to measure the error is to compute the increased distance traveled via the midpoint compared to the distance between the original points. If the increase is tiny compared to the original distance, we have little to complain about. I have computed these errors using the accurate ellipsoidal distance for the WGS84 ellipsoid. As a typical example of the results, here is a plot of relative errors when one of the endpoints is fixed at (lat, lon) = (45, 0):

Relative errors

The contours are on a logarithmic (base 10) scale: the -6 contours show points where the relative error is 10^(-6); that is, one part per million (ppm). The -5 contours (barely visible near (-45, 180), the diametrically opposite point) are 10 ppm. The -7, -8, etc. are fractions of a ppm: highly accurate.

Evidently, as long as we're not trying to compute midpoints of two points that are almost diametrically opposite, we will do fine. (Remember, the calculation is perfectly correct for the sphere; these errors are due to the flattening of the sphere.)

Given that 16 bit precision is about 16 ppm (a base-10 log equal to -4.8), it looks ok to use the spherical formula for midpoint-finding provided the two points are more than one degree away from being diametrically opposite.

What about the simpler linear formula? To study this, let's compare the distance between the linear midpoint (obtained by averaging the two latitudes and two longitudes) to the spherical midpoint, relative to the distance between the two endpoints. The next figure fixes one endpoint at (45, 180) and explores a relatively small region around it.

Relative errors 2

Most of these contours (base-10 logarithms once again) are near -2: that's one part per hundred (1%) error. For north-south directions there is no error, but for all other directions the error is unacceptable for many applications.

To see whether the linear approximation ever becomes OK, let's zoom into that previous map by a factor of 10. Now it is one degree wide (50 miles at this latitude) and one-half degree across (35 miles): we're looking at the scale of a large city or a small city and its suburbs.

Now the contours are around -3 to -4: that's 100 to 1000 parts per million (0.01% to 0.1%). Pretty crude and just barely noticeable on a high-res computer screen if you look closely.

Relative errors 3

Looking back, it is apparent that the slightly more complicated--but still easily implementable--spherical formula achieves greater accuracy across the globe than the simple linear formula does even in nearby locations. (I'm fudging a little, because I have used two different ways to measure errors, so they are not directly comparable.)

The bottom line:

--The linear formula will misposition the midpoint by a relative error of 0.01% to 0.1% at a city scale; over larger areas, the mispositioning can be grossly wrong (1% on up to hundreds of %).

--The spherical formula is absolutely correct for a spherical earth model. Compared to the more accurate ellipsoidal formula, it should still do well except for points that are almost diametrically opposite.

3
  • 1
    Wow, thanks! That's some really good info. I don't need perfect precision, it's just an iPhone app for personal GPS use, mainly only focused on streets. Since a highway isn't perfectly straight for a hundreds kilometers even linear precision is "good enough". But political boundaries are often straight and long (the right side of Western Australia is 1,800km). The spherical formula should be perfect for that. I wish I could give you a bounty for your answer, you deserve it. Jan 13, 2012 at 4:36
  • 1
    Note to others who want to implement the pseudocode above: trigonometric functions in most programming languages tend to take radians instead of degrees, and the argument order for atan2 tends to be the inverse of the order used in the pseudocode above. Here's an implementation in Swift: gist.github.com/ts95/bf113ba080c36d9ff645b322081360d6 Oct 21, 2020 at 10:53
  • Thank you @Toni: that's a good point about atan2; I should have noted that in this post.
    – whuber
    Oct 21, 2020 at 12:29
8

Your solution works for small distances, but it won't work for larger ones. The easiest way to see why is to look at the following map (taken from here):

enter image description here

It is is a world map in equirectangular projection - longitudes and latitudes are simply linearly projected into X and Y axes. Your math can be represented by the blue rectangle. However, the blue diagonal does not represent the shortest "line" between two points on Earth's surface. What you need is a great circle (represented by black curve).

Although you're using C++, I recommend taking a look at Mike Gavaghan's C# Geodesy Library, it's open source and the code is of high quality, you should be able to figure out how to do the calculation along the great circle. You can also take a look at Vincenty's formulae if you're interested in math behind the calculation.

1
  • Thanks! That's exactly what I'm looking for. Mike Gavaghan's code looks easy to port to C (i'm actually using Objective-C btw). I was pretty sure my linear math would be way out, but needed to have something. Jan 11, 2012 at 23:42
-1

Simple linear math seems to work fine:

double newLat = lastCoord.latitude + ((coord.latitude - lastCoord.latitude) / 2);
double newLon = lastCoord.longitude + ((coord.longitude - lastCoord.longitude) / 2);

I'm recursively inserting half way new points in each section of the way that's too large, until they're all small enough.

I'm not sure if more complex math is needed for larger distances (can anyone confirm this for me?), but it works fine over this one (about 5km):

example of long straight line between points

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.