Calculate angle/heading, clockwise relative to North, between ECEF points Let there be two data points, in ECEF, as follows (in meters):


*

*A1 (x1, y1, z1)

*A2 (x2, y2, z2)


I want to calculate the angle, Clock-wise relative to North (the Z-axis in ECEF), that is the heading from point 1 to point 2.
I can easily visualize this, and solve it, in 2D.
This question (https://stackoverflow.com/questions/642555/how-do-i-calculate-the-azimuth-angle-to-north-between-two-wgs84-coordinates) looks very similar to my question, yet the selected answer assumes you're on a flat surface.
This question (https://stackoverflow.com/questions/8123049/calculate-bearing-between-two-locations-lat-long) is in lat/lon, so again just 2D.
Can I just simply project my two points onto the Y-Z plane? For example:
theta = atan2(y1-y2, z1-z2) 
and then adjust for the different Earth "quadrants", if that makes sense.
I feel like solving this problem is much more complicated than this. I don't feel confident that my "projection" idea is correct... but I can't visualize it.
 A: What is “heading”? To me it means you establish a local coordinate system at point $A_1$, and giving the azimuth (and perhaps elevation) of $A_2$ with respect to this coordinate system. So what directions do you have in this coordinate system?


*

*Up is radial. The unit up vector is $A_1$, scaled to unit length.

*East is perpendicular to up and perpendicular to the axis of the earth. So compute the cross product of the vector $(0,0,1)$ and the up vector, normalize the result to unit length and make sure you got the correct sign.

*North is perpendicular to up and perpendicular to east, so take their cross product. The norm will already be one (since the other two are length one and perpendicular to one another), so just make sure you got the correct sign.


(Getting all signs correct means taking the cross product in the right order. There are many possible conventions. Instead of discussing this at length, just try out one generic example and if the signs are wrong, swap the order in the cross product.)
So now you have three basis vectors of an orthonormal basis, and you can compute the coordinates of your points with respect to this basis simply by computing the corresponding dot products. If I call these coordinates $u,e,n$ (in direction up, east resp. north), then the azimuth would be
$$\tan\theta=\frac en\qquad\theta=\operatorname{arctan2}(e,n)$$
with zero meaning north and positive angles going east from there. If you point a laser in the given direction at $A_1$ and adjust the elevation correctly, you get it to point at $A_2$ unless something (like the earth) gets in the way.
