[Math] Calculating a perpendicular distance to a line, when using coordinates (latitude & longitude)

algorithmstrigonometry

I'm trying to implement the Douglas-Peucker algorithm for simplifying a recorded GPS track (a list of coordinates). All implementations I can find assume a simple X/Y grid of squares, however ideally I would like to adapt it to work with Earth coordinates which are not square.

The algorithm relies on working out a points perpendicular distance from a line defined by two other points. As such, I need to create a function which can calculate the perpendicular distance to the line using latitude and longitude. Currently the function assumes a a simple X/Y grid.

I cannot work out how to adapt the function I have to use something like the Haversine formula for calculating distance between two points, as the existing method doesn't use the perpendicular point.

Is is possible to adapt this function to make it more accurate when applied to GPS coordinates? The easiest way to do this would be to be able to calculate the coordinates of the point on the line which is perpendicular to point p as I could then use my existing Haversine function to calculate the distance. How could I approach this if I don't know the coordinates of the perpendicular point? (or do I not need to?)

In my diagram, the points a and b define the line, and point p is the point where I need to determine its distance (in meters) from the line.

The current function is as follows:

public static double PerpendicularDistance(TrackPoint pointA, TrackPoint pointB, TrackPoint pointP)
    {
        // Area = |(1/2)(x1y2 + x2y3 + x3y1 - x2y1 - x3y2 - x1y3)|   *Area of triangle
        // Base = √((x1-x2)²+(x1-x2)²)                               *Base of Triangle*
        // Area = .5*Base*H                                          *Solve for height
        // Height = Area/.5/Base

        double area = Math.Abs(.5 * (pointA.Lat * pointB.Lon + pointB.Lat * pointP.Lon + 
            pointP.Lat * pointA.Lon - pointB.Lat * pointA.Lon - pointP.Lat * pointB.Lon - pointA.Lat * pointP.Lon));
        double bottom = Math.Sqrt(Math.Pow(pointA.Lat - pointB.Lat, 2) + Math.Pow(pointA.Lon - pointB.Lon, 2));
        double height = area / bottom * 2;
        return height;
    }

enter image description here

Best Answer

  1. I'll tell you how to do this in a moment.

  2. First, I'll tell you why it's a bad idea.

Your data, unless it comes from, say, an orbiting satellite with a huge velocity, consists of points whose distance is very small compared to the radius of the earth. The error introduced by treating lat/lon coordinates as $xy$-coords (after suitable scaling for the compression of longitude away from the equator) will generally be tiny compared to the error you make in simply applying the D-P algorithm. The exceptions are (1) the poles, (2) very fast motion, which probably should not be simplified because of the 3D analogue of "aliasing" in computer graphics, and (3) the international dateline. For cases 1 and 3, the best solution is probably to change to a different set of lat/long coords. For the dateline, for instance, you can add 180 to all longitudes, and then subtract 360 until the result is in the range -180 to 180. Most of your points will end up with "new longitudes" near zero. Once you've done the DP algorithm, you convert back.

For polar regions, you switch to a lat-lon system whose "poles" are two opposite points near what we call the equator, and whose "Greenwich meridian" is half of the thing we call the equator. This involves converting to rectangular coordinates, swapping y and z, and converting back.

OK. Now the "how to":

Given a lat-long point with latitude $u$ and longitude $v$ (which I'll write $[u, v]$) you can compute $xyz$ coordinates via the following: \begin{align} x &= \sin u \cos v \\ y&= \sin v\\ z &= \cos u \cos v \end{align}

I'll write the result as $(x, y, z)$.

Your problem is that you have points $A, B, C$ in $[u, v]$ coordinates, and want to project $C$ onto the arc between $A$ and $B$, and find the length of the projection arc.

To do this, I'm going to talk about $A$, $B$, and $C$ as if they've been converted to $xyz$ triples already, to avoid new notation. I'll use the dot-product of two triples, $$ (x_1, y_1, z_1) \cdot (x_2, y_2, z_2) = x_1 y_1 + x_2 y_2 + x_3 y_3 $$ and the "cross product," which I'll let you look up on Wikipedia for fear that I'll get a sign wrong.

OK.

a. Compute $n = A \times B$ (the cross product) and $$ N = n / \sqrt{n \cdot n} $$ Explanation: There's a plane passing through $A$, $B$, and the earth's center; the line from the earth's center to $N$ is perpendicular to this plane, and $N$ is on the surface of the earth.

b. Compute $$ s = 90^\circ- |\arccos (C \cdot N)|. $$

Explanation: this is the angular distance between a ray from the earth's center to $C$ and the plane described above.

You can also compute the "distance" between $A$ and $B$ as $$ s' = \arccos(A \cdot B) $$

Assuming that you're working in degrees, this will range from $0$ to $180$. If you want miles, you need to multiply this by the circumference of the earth, about 25,000 miles. (All of this assumes that the earth is spherical -- another source of error).

Anyhow, this gives you an alternative approach to the DP algorithm: simply use the formulas I've given you for $s$ and $s'$ as the "distances" that DP requires. They'll all be between 0 and 180 degrees (or $0$ and $\pi$ radians), but they'll be proportional (ignoring non-sphericity) to earth's-surface distances, which is all that you need.

And then be careful not to apply the algorithm in any situation in which a distance is greater than about, say, 15 degrees, because any assumptions of "flatness" in D-P is likely to be broken by the curvature of the earth at that scale.

Related Question