Converting Lat/Long coords to Cartesian X/Y, then calculating shortest distance between point & line segment

I'm having an issue with accuracy when converting Lat/Long coordinates to X,Y and then finding the shortest distance from a Point to a Line with said coordinates.

The distance is off by around 40-50% of actual, which is unaccceptable for use.

First I convert the coordinates (which are in decimal format) to radians, and then X/Y (hope you guys don't mind some C# code):

private const double EarthRadius = 6367 / 1.61 * 5280;

private static double GetXCoord(double lat, double lon)
{
}

private static double GetYCoord(double lat, double lon)
{
}

public static double ToRadians(this double valueInDegrees)
{
return (Math.PI / 180) * valueInDegrees;
}


Using these methods, I get my point (x0, y0) and my line segment (x1, y1) (x2, y2).

I then perform the Point-Line Distance calculation found here: http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html

    var lateralPointCalc1 = Math.Abs(((x2 - x1) * (y1 - y0)) - ((x1 - x0) * (y2 - y1)));
var lateralPointCalc2 = Math.Pow((x2 - x1), 2) + Math.Pow((y2 - y1), 2);

lateralPointCalc2 = Math.Sqrt(lateralPointCalc2);

lateralPointDistance = lateralPointCalc1 / lateralPointCalc2;
return lateralPointDistance;


For test coords:

28.503946 / -99.453589 (x0, y0) (my point)

28.485044 / -99.453709 (x1, y1)

28.49823 / -99.46834 (x2, y2)

I would expect the shortest distance to be around ~4930 feet, but my method returns ~2969 (off by 2000 feet, which is a huge margin of error.)

What am I doing wrong here? Any help would be greatly appreciated, thanks!

-
I'm not sure what you are trying to do. Surely you need to compute (x,y,z) coordinates, not just (x,y)? You need to compute z, something like: EarthRadius * Math.Sin(lat.ToRadians()) – copper.hat Apr 16 '12 at 20:38
The linear approximation (using a tangent plane) is not very good for larger distances. Use spherical great arc distance. This is also only an approximation, since the earth is closer to being an oblate spheroid, but it should give you much better accuracy than linear! – bgins Apr 16 '12 at 20:46
@bgins: The result is much worse that it would be for a linear approximation on a tangent plane. The problem is rather, as copper.hat pointed out, that these are just the $x$ and $y$ coordinates, so this is not projected onto a tangent plane but onto the equatorial plane, which is completely wrong. A linear approximation on a tangent plane may well have sufficient accuracy, depending on the application. – joriki Apr 16 '12 at 20:54
Okay, I've seen info on getting the Z coordinate, but I'm not sure what equation to use with (X, Y, Z) to find the shortest distance from a Point to a Line Segment. Would this be what I'm looking for? mathworld.wolfram.com/Point-LineDistance3-Dimensional.html -- As far as accuracy, I'll be working with very short distances, so I imagine a "linear approximation on a tangent plane" would be sufficient. – ParaSwarm Apr 16 '12 at 21:00
@joriki: yes, thanks, i saw that. So I was trying to type up a start of an explanation of how the third dimension fits in. – bgins Apr 16 '12 at 21:28

The simplest general-purpose mathematical model would be to use spherical coordinates with standard latitude $\phi$ (north of equator) and longitude $\theta$ (say east of Greenwich): $$X(\theta,\phi) =\left[\matrix{x\\y\\z\\}\right] =\left[ \matrix{ r\cos\phi\cos\theta\\ r\cos\phi\sin\theta\\ r\sin\phi }\right] =r\left[ \matrix{ \cos\phi\cos\theta\\ \cos\phi\sin\theta\\ \sin\phi }\right]$$ From an arbitrary point $X$ as above, you can compute the east and north directions on the tangent plane from the partial derivatives $X_\theta$ and $X_\phi$, respectively (useful for local comparisons with compatible map projections): $$X_\theta =\frac{\partial}{\partial\theta}X(\theta,\phi) =r\left[ \matrix{ -\cos\phi\sin\theta\\ \cos\phi\cos\theta\\ 0 }\right] \qquad\text{(east)}$$

$$X_\phi =\frac{\partial}{\partial\phi}X(\theta,\phi) =r\left[ \matrix{ -\sin\phi\cos\theta\\ -\sin\phi\sin\theta\\ \cos\phi }\right] \qquad\text{(north)}$$ The zenith (outward radial vector direction) is of course just $$X_r =\frac{\partial}{\partial r}X(\theta,\phi) =\left[ \matrix{ -\sin\phi\cos\theta\\ -\sin\phi\sin\theta\\ \cos\phi }\right] \qquad\text{(zenith).}$$ Now given two points $X_1(\theta_1,\phi_1)$ and $X_2(\theta_2,\phi_2)$, the great arc between them is on a circle intersecting the (assumed) spherical surface of the earth (with assumed constant radius $r$, which is a simplification) and the plane through the center of the earth (the origin of our coordinate system) and the two points $X_1$ & $X_2$. The normal to this plane can therefore be found by taking the cross product: $$N = X_1 \times X_2$$ And the distance along this great arc is just the arc length from the angle $\alpha$ between the radial vectors of the two points, which is given by the arc cosine of their dot product (a caveat, however, is that it is ill-conditioned for angles near $0$ or $\pi\text{ rad}=180^\circ$): $$s=r\alpha\qquad\text{for}\qquad\cos\alpha=X_1 \cdot X_2$$ If you want the bearing (starting direction on the map) from $X_1$ to travel on this great arc path, well that, too, can be easily found. Take the tangent vector $V$ at $X_1$ in the direction of $X_2$ by crossing the tangent plane's normal, $N$, with $X_1$ (possibly up to sign): $$V = N \times X_1 = (X_1 \times X_2) \times X_1$$ Plugging in the above definitions or using some vector and trigonometric identities will give you explicit formulas for each of these in terms of $\theta_i$ and $\phi_i$.

So now suppose you have the point $X_0$ and the line (or great arc) through points $X_1$ and $X_2$ on the surface of the sphere. To find the distance from $X_0$ to this great circle, we essentially need to project (or better, rotate) $X_0$ onto the plane of the great arc (or onto the great arc itself). First we need to find the closest point, $F$ (for 'foot'), to $X_0$ on that arc. This will be the normalized version of $$M = N \times \left(X_0 \times N\right)$$ scaled by the radius $r$, $$F = r\,\widehat{M} = r\,\frac{M}{||M||}$$ and thus our distance $d=r\beta$, from $X_0$ to the "line" through $X_1,X_2$, can be obtained from the great arc angle $\beta$ from $X_0$ to $F$, by taking the arc (inverse) cosine of $$\cos\beta=X_0\cdot\widehat{M}=\frac{X_0\cdot M}{||M||}=X_0\cdot\frac{F}{r}\,.$$ Another explanation can be found at this stackoverflow post.

You can find more background on this, with good diagrams, under wikipedia's spherical coordinates, geometry, trigonometry, law of cosines, half-side formula and haversine formula pages, along with some pages to get a sense for the physical realism of this model such as the earth's radius, among others.

Playing around in sage, I got a distance of $4841.93165384$ feet, or $1613.97721795$ yards, or $1475.72976066$ meters:

var('x,y,z,lat,lon')
rkm = 6371 # km
rmi = 3959 # mi
d2r = pi / 180
ll = Matrix(RDF, 3, 2, [\
28.503946, -99.453589,\
28.485044, -99.453709,\
28.498230, -99.468340])
xyzll = lambda lat,lon: [\
(cos(d2r*lat)*cos(d2r*lon)).n(),\
(cos(d2r*lat)*sin(d2r*lon)).n(),\
(sin(d2r*lat)).n()]
xyz = Matrix(RDF, 3, 3)
for i in range(3):
xyz[i] = xyzll(ll[i,0],ll[i,1])
print xyz
v0 = vector(xyz[0])
v1 = vector(xyz[1])
v2 = vector(xyz[2])
vn = v1.cross_product(v2); vn=vn/vn.norm(); print 'vn:', vn, vn.norm()
vm = vn.cross_product(v0.cross_product(vn)); vm=vm/vm.norm(); print 'vm:', vm, vm.norm()
a0 = arccos(vm.dot_product(v0)); print 'a0:', a0, (a0/d2r).n()
print 'distance:', a0*rmi, a0*rkm

# The OUTPUT:
[-0.144339114164  -0.86684945361  0.477219283871]
[-0.144366780752 -0.867004401598  0.476929345107]
[-0.144570113877 -0.866859220201  0.477131611326]
vn: (-0.760935788228, -0.21083794795, -0.613615584791) 1.0
vm: (-0.144515375391, -0.866898313757, 0.477077163445) 1.0
a0: 0.000231632359231 0.0132715565826116
distance: 0.917032510197 1.47572976066

-
Thanks! This definitely gives me something to work with. Will come back to upvote when I have the required rep. – ParaSwarm Apr 17 '12 at 3:42
Bgins, I've implemented this in code, but it appears to consider the "line segment" as not having an end. Occasionally this will find the shortest distance from the point to an invisible line that continues to extend past the endpoint of the line segment. Is there no way around this? – ParaSwarm Apr 19 '12 at 15:11
The foot $F$ is the closest point from $X_0$ to the great circle through the two other points. I think $F$ lies outside the arc between $X_1$ & $X_2$ unless the angles $\angle FOX_i$ have opposite sign and magnitudes summing to $\angle X_1OX_2 < \pi$ rad$=180^\circ$. There is surely a clever way to test for this, but I can't think of it right now. – bgins Apr 19 '12 at 22:10
Alright, then what if we approach it this way (forgive my poor mathspeak): We draw two lines, each originating on the endpoints of our original line segment. These new lines are angled 90 degrees from our original line segment, and eventually end up back where they originated. Then we do a check to see if our point x0 resides between these two lines. If so, then we perform the calculation from your answer. If not, then I do a simple point-to-point distance check (Haversine) from x0 to the closest point on the line segment. Any idea how to accomplish this, if it makes sense? Thanks again! – ParaSwarm Apr 20 '12 at 17:50
Scratch that. I'll just need to find the point given by your above formula, then check if it's lat/long is greater than/less than the points on the line segment. – ParaSwarm Apr 20 '12 at 19:35