[Math] Coordinates of interception point Y with XY being the shortest distance of X to AB on sphere

analytic geometryspherical coordinatesspherical-geometry

How would one calculate the interception point $Y$ with $\overleftrightarrow{XY}$ being the shortest distance of $X$ to $\overleftrightarrow{AB}$?

Cross track error coordinate

This answer to the question How to find the distance between a point and line joining two points on a sphere? describes how to calculate the distance, but not the coordinates.

What would be the most direct way to calculate the coordinates of interception point $Y$?

I since found out that this problem is also called the along track distance and have implemented it but have troubles checking my results.

For completeness here are the formulas and their code counterpart I have used and implemented:

$dxt = asin(\sin(d_{13}/R) * \sin( \phi_{13}−\phi_{12})) * r$

with

$d_{13}$ is distance from start point to third point
$\phi_{13}$ is (initial) bearing in degrees from start point to third point
$\phi_{12}$ is (initial) bearing in degrees from start point to end point
$r$ is the earth’s radius

in JavaScript:

var dXt = Math.asin(Math.sin(d13/R)*Math.sin(Math.PI / 180 * (brng13-brng12))) * R;

in Java:

final double dt = Math.asin(Math.sin(d13 / r) * Math.sin(Math.toRadians(b13 - b12))) * R;

then I calculate the Along Track Distance (following the formula on the site)

$dat = acos(\cos( d_{13}/r) / \cos(dxt/r)) * r$

where

$d13$ is distance from start point to third point
$dxt$ is cross-track distance
$r$ is the earth’s radius

in JavaScript:

var dAt = Math.acos(Math.cos(d13/R)/Math.cos(dXt/R)) * R;

in Java:

final double dAt = Math.acos( Math.cos(d13 / r) / Math.cos(dXt / r) ) * r;

As a test case I have

$A (0,0)$
$B (90,0)$
$X (45,45)$

and get the following results my with my implementation

$3335847.799336763$ for $\overleftrightarrow{XY}$
$6086322.17407125$ for $\overleftrightarrow{AB}$

and thus $(54.735, 0) as coordinates for Y

Using $r * \cos (lat)$ the radius on the given latitude is $4504977.303$, thus $1/8$ of the circumference is $3538200.9$.

Shouldn't $Y$ be $(0, 45)$? Surely the shortest distance from $X (45, 45)$ to $\overleftrightarrow{AB}$ is the great circle on 45° latitude ? Thus the cross track -distance should be $3538200.9$. Am I lacking imagination in three-dimensions?

Best Answer

I think the best way to do this would be to look at the sphere as a subset of 3-dimensional space, find the coordinates of A, B, and X in that 3-dimensional space (that's just converting spherical coordinates to Cartesian), and then you can (orthogonally) project everything onto the plane containing A, B, and X. A, B, and X project to themselves, of course, and the arc between A and B will project to the line between A and B, and the arc between X and Y will project to the line between X and Y. Therefore, Y will project to the orthogonal projection of X onto the line AB, and then from that you can just project back onto the sphere and convert that to latitude and longitude. It's somewhat indirect, but it avoids having to work with distances in spherical coordinates.

Related Solutions

[Math] Calculating longitude degrees from distance

The problem is that you must have made a mistake while copying your formula in Mathematica. After extracting $\Delta\lambda$ from the formula of $d$, knowing $\phi$ constant, I obtain (when $\phi\neq\pi/2 + k\pi$): $$ \Delta\lambda = \pm 2 \arcsin\left\lvert \frac{\sin\frac{d}{2r}}{\cos\phi_1} \right\lvert $$

I must say I'm concerned that in this formula, the argument of arcsin might become greater than 1.. I did this rapidly but I might try to check how this works. Oh, and use radians with trig functions.

Edit: Here is your example with Matlab:

>> phi = 40 * pi/180; a = 10 * pi/180; b = 11*pi/180;
>> r = 6371;
>> d = 2*r*asin( sqrt( cos(phi)^2 * sin( (a-b)/2 )^2 ) )
d =  85.180
>> delta = 2*asin( abs( sin(d/(2*r))/cos(phi) ) )
delta =  0.017453
>> b-a
ans = 0.017453

[Math] Determine if one point lies between two other points on a sphere

For this answer I will not rely on Cartesian coordinates, but I will not work directly in terms of latitude and longitude either. Instead, it will be convenient to represent a point with spherical coordinates $(\phi,\theta)$ as $z=e^{i\phi}\tan\dfrac\theta2$, corresponding to a stereographic projection of the sphere into the complex plane. We can now apply Mobius transformations to move the complex numbers $\alpha$, $\beta$ (corresponding to $A$, $B$) to more helpful locations.

Specifically, we can use rotations to move $\alpha$ to zero, then map $\beta$ to the real line, and finally rotate $\{\alpha,\beta\}$ to $\{r,r^{-1}\}$ on the positive real line. This amounts to rotating our coordinates so that $A$ and $B$ have zero longitude and opposite latitudes. To make this simpler, we shall assume that $\alpha$ is real without loss of generality (i.e. we take the longitude of $B$ relative to $A$).

We now state these transforms explicitly. First, note that the subgroup of the Mobius group corresponding to rotations consists of mappings of the form $z\mapsto \dfrac{u z-\overline{v}}{v z+\overline{u}}$ for complex $u,v$. A mapping of this form which takes $\alpha$ to the origin is $z\mapsto \dfrac{z-\alpha}{\alpha z+1}$; this also maps $\beta\mapsto \dfrac{\beta-\alpha}{1+\beta\alpha}= R e^{i\psi}$ for some real $R,\psi$. A rotation $z\mapsto e^{-i\phi} z$ then brings $\beta$ to $R$.

For the last transformation, consider the mapping $z\mapsto \dfrac{1+z r}{r-z}$ for real $r>0$ which takes the real line to itself. Then $\{0,R\}\mapsto \left\{\dfrac{1}{r},\dfrac{1+r R}{r-R}\right\}$; since we want the images of $0$ and $R$ to be on equal and opposite sides of the 'equator', we require their product to be unity. From this we deduce $r=R+\sqrt{1+R^2}$. Composing the three transformation gives the net rotation

$$\mathcal{T}:\;z \mapsto \dfrac{z-\alpha}{1+\beta\alpha} \mapsto e^{i\phi}\left(\frac{z-\alpha}{1+\beta\alpha}\right) \mapsto \dfrac{1+e^{i\phi}\left(\frac{z-\alpha}{1+\beta\alpha}\right)r}{r-e^{i\phi}\left(\frac{z-\alpha}{1+\beta\alpha}\right)} =\dfrac{e^{-i\phi}(1+\beta\alpha)+(z-\alpha)r}{r e^{-i\phi}(1+\beta\alpha)-(z-\alpha)}$$

We now can give a criterion for the original question: a point $z$ lies between $\alpha$ and $\beta$ if $r^{-1}\leq | \mathcal{T}(z)|r$. In principle one can do inverse transformations to map this annulus back to the original coordinates, but I won't do so explicitly. Instead I'll just summarize the result:

Let $\alpha,\beta$ be the stereographic coordinates of two points $A,B$ on the unit sphere, and define $w=R e^{i\phi}:=\dfrac{\beta-\alpha}{1+\beta\alpha}$. Then a point $C$ on the unit sphere lies between $A$ and $B$ iff its sterographic coordinate $z$ satisfies $$\frac{1}{r}< \left|\dfrac{e^{-i\phi}(1+\beta\alpha)+(z-\alpha)r}{r e^{-i\phi}(1+\beta\alpha)-(z-\alpha)}\right|<r$$ where $r=R+\sqrt{1+R^2}$.

Best Answer

Related Solutions

[Math] Calculating longitude degrees from distance

[Math] Determine if one point lies between two other points on a sphere

Related Question