2
$\begingroup$

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!

  • 0
    I think I have an answer you can use now to get a reliable estimate even over large distances on the sphere (without too much much extra computational work), assuming constant (e.g. mean) earth radius $r$. You only need to convert the points from from spherical $(\theta,\phi)=$(longitude,latitude) to Cartesian $(x,y,z)$ coordinates, compute some cross & dot products, divide by a vector norm, take the arc (inverse) cosine and multiply by $r$ to get your answer. Let me know if there are still questions. Also, you might want to encapsulate things like points in each coordinate system in a class.2012-04-16

1 Answers 1

3

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 
  • 0
    But what is it you need, then, if $X_0$ and its foot $F$ aren't inside this longitudinal sector of the (rotated) sphere?2012-04-20