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'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())2012-04-16
  • 0
    The linear approximation (using a tangent plane) is not very good for larger distances. Use spherical [great arc distance](http://en.wikipedia.org/wiki/Great-circle_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!2012-04-16
  • 0
    @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.2012-04-16
  • 0
    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? http://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.2012-04-16
  • 0
    @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.2012-04-16
  • 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
    Thanks! This definitely gives me something to work with. Will come back to upvote when I have the required rep.2012-04-17
  • 0
    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?2012-04-19
  • 0
    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.2012-04-19
  • 0
    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!2012-04-20
  • 0
    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.2012-04-20
  • 0
    I think you could check whether $X_0$ is inside the "rhombus" with vertices $X_1,X_2$ & $\pm r\widehat{N}$ and great arc edges. This rhombus is divided by the diagonal between $X_1,X_2$ into two right spherical triangles, so we can say [a lot](http://en.wikipedia.org/wiki/Spherical_trigonometry#Napier.27s_Pentagon) about it. This diagonal is like the arc of a rotated equator, and the points $\pm r\widehat{N}$ are like rotated poles. Within this rotated coordinate system, we're asking whether the longitude of $X_0$ lies between that of the other points, and we have to be careful of periodicity.2012-04-20
  • 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