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)
{
return (EarthRadius * Math.Cos(lat.ToRadians()) * Math.Cos(lon.ToRadians()));
}
private static double GetYCoord(double lat, double lon)
{
return (EarthRadius * Math.Cos(lat.ToRadians()) * Math.Sin(lon.ToRadians()));
}
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!
Best Answer
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: