North-bearing between two points in Cartesian space

Question

I have two points $P$ and $Q$ given in earth-centered-earth-fixed Cartesian space $(XYZ)$

How can I compute the compass-bearing of $P$ towards $Q$ without going into Geodetic space (Lat, Lon)?

Answer 1

I'll assume what I think is the standard convention for ECEF coordinates: right-hand $x,$ $y,$ $z$ coordinates with the positive $x$ axis going through $0$ latitude $0$ longitude, the positive $y$ axis at $0$ latitude $90$ degrees east, and the positive $z$ axis at $90$ degrees north (the conventional north pole).

Let $p$ be the vector from the origin to $P$.

Let $\hat e$ be a unit vector in the equatorial plane orthogonal to $p.$ There are two such vectors; you want the one that is “east” of $P.$ One way to find $\hat e$ is to project $p$ onto the equatorial plane, normalize it to unit length, and rotate it $90$ degrees eastward. If $p = (x_P,y_P,z_P)$ then $$ \hat e = \frac1{\sqrt{x_P^2+y_P^2}} \left( -y_P, x_P, 0 \right). $$

An alternative procedure giving the same result is to take the cross product with a vector in the positive $z$ direction and then normalize the result to a unit vector: $$ \hat e = \frac1{\lVert p\times(0,0,1)\rVert} p\times(0,0,1). $$

Let $\hat n$ be a unit vector orthogonal to $p$ and $\hat e,$ selecting the more “northward” vector from the two possible vectors. There are various ways to get such a vector; a simple formula (but not necessarily the most efficient computation) is the cross product $$ \hat n = \frac1{\lVert P\rVert} P \times \hat e. $$

Now take inner products \begin{align} N = q\cdot\hat n,\\ E = q\cdot\hat e. \end{align} where $q$ is the vector from the origin to $Q.$

The course angle you want is the angle $\theta$ such that $E=r \sin \theta$ and $N=r \cos \theta$ for some positive number $r.$ (This implies that $r = \sqrt{N^2 + E^2},$ but you should not need to compute $r.$) There are various ways to find $\theta$ using some version or another of the inverse tangent function. I like the function $\mathrm{atan2}(E,N)$ found in many computer environments, which gives you the angle directly in radians.