Spherical or Astronomical Coordinate Transformation

Converting from Galactic coordinates to altitude and azimuth involves three separate coordinate transformations: Galactic longitude and latitude to equatorial right ascension and declination [(${\displaystyle l}$, ${\displaystyle b}$) ${\displaystyle \rightarrow }$ (${\displaystyle \alpha }$, ${\displaystyle \delta }$)]; to hour angle and declination [(${\displaystyle \alpha }$, ${\displaystyle \delta }$) ${\displaystyle \rightarrow }$ (${\displaystyle ha}$, ${\displaystyle \delta }$)]; to azimuth and altitude [(${\displaystyle ha}$, ${\displaystyle \delta }$) ${\displaystyle \rightarrow }$ (${\displaystyle az}$, ${\displaystyle alt}$)].

It is possible (and recommended) to implement these coordinate transformations using rotation matrices. In this method, you generate a vector in the original coordinate system; convert the vector to another coordinate system by rotating the coordinates using matrix multiplication; and convert the vector to the angles of the new coordinate system. The method is general and can be applied to any coordinate transformation.

1 Rotation Matrices: The Method

To restate the problem: we start with (${\displaystyle long}$, ${\displaystyle lat}$) in one coordinate system and want to convert to (${\displaystyle long^{\prime }}$, ${\displaystyle lat^{\prime }}$) in some other coordinate system. Here’s the prescription:

Step 1. First, convert the angles to rectangular coordinates. One would usually call these (${\displaystyle x}$, ${\displaystyle y}$, ${\displaystyle z}$); here, to emphasize the vector/matrix flavor, we call them (${\displaystyle x_{0}}$, ${\displaystyle x_{1}}$, ${\displaystyle x_{2}}$) and denote the 3-element vector ${\displaystyle {\vec {x}}}$. To accomplish this conversion in Python:

x = numpy.array([0.,0,0])
x[0] = numpy.cos(lat) * numpy.cos(long)
x[1] = numpy.cos(lat) * numpy.sin(long)
x[2] = nump.sin(lat)


Step 2. Apply the rotation matrix ${\displaystyle {R}}$ (we’ll discuss its definitions below):

${\displaystyle {{\vec {x}}^{\prime }}={R\cdot {\vec {x}}}\;.\,\!}$

In Python, we’ll use numpy arrays to represent ${\displaystyle {\vec {x}}}$ and ${\displaystyle {{\vec {x}}^{\prime }}}$, and you do matrix multiplication with numpy.dot:

xp = numpy.dot(R, x)


Step 3. Convert the primed rectangular coordinates to the new set of spherical coordinates ${\displaystyle (long^{\prime },lat^{\prime })}$:

{\displaystyle {\begin{aligned}long^{\prime }&=\tan ^{-1}\left({{{\vec {x}}^{\prime }}_{1} \over {{\vec {x}}^{\prime }}_{0}}\right)\;,\\lat^{\prime }&=\sin ^{-1}({{\vec {x}}^{\prime }}_{2})\;.\end{aligned}}\,\!}

To do these in Python, where we’ll write longp, latp:

longp = numpy.arctan2(xp[1], xp[0])
latp = numpy.arcsin(xp[2])


That’s it! If you want to “go backwards”, you just apply the matrix multiplications using the inverse matrices. For rotation matrices, the inverse is always equal to the transpose(!)[1]—symbolically for the matrix ${\displaystyle {R}}$, ${\displaystyle {R^{-1}}={R^{T}}}$. So to go from the primed system to the unprimed, you switch the primed and unprimed in equation (2) and use the inverse rotation matrix

${\displaystyle {\vec {x}}={R^{-1}\cdot x^{\prime }}={R^{T}\cdot x^{\prime }}\;.\,\!}$

and in Python…

iR = numpy.linalg.inv(R)
x = numpy.dot(iR, xp)


2 Rotation Matrices: Specific for our Problem

2.1 (Ra, Dec) to (Ha, Dec)

Converting from ${\displaystyle (\alpha ,\delta )\rightarrow (ha,\delta )}$ keeps the declination the same and uses the relationship ${\displaystyle ha=LST-\alpha }$. It is easiest to think of this in two steps, so we express

${\displaystyle {R}_{(\alpha ,\delta )\rightarrow (ha,\delta )}={R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),2}\cdot {R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),1}\;.\,\!}$

First, we rotate around the equatorial pole by an angle equal to the Local Sidereal Time (${\displaystyle LST}$), which does ${\displaystyle \alpha \rightarrow (\alpha -LST)}$:

{\displaystyle {\begin{aligned}{R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),1}=\left[{\begin{array}{ccc}\cos(LST)&\sin(LST)&0\\-\sin(LST)&\cos(LST)&0\\0&0&1\\\end{array}}\;\right]\;.\end{aligned}}\,\!}

Next, the ${\displaystyle ha}$ and ${\displaystyle \alpha }$ go in opposite directions, which is equivalent to converting from the original left-handed to a right-handed coordinate system, so the second step is just to perform this reversal:

{\displaystyle {\begin{aligned}{R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),2}=\left[{\begin{array}{ccc}1&0&0\\0&-1&0\\0&0&1\\\end{array}}\;\right]\;.\end{aligned}}\,\!}

The full rotation matrix is the matrix product ${\displaystyle {R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),2}{\cdot }{R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),1}}$. Note the order! Applying ${\displaystyle {R}_{(\alpha ,\delta )\rightarrow (ha,\delta ),2}}$ at the beginning in the matrix product means that it operates last on the vector ${\displaystyle {\vec {x}}}$, which is what we want. So we have as the product…

{\displaystyle {\begin{aligned}{R}_{(\alpha ,\delta )\rightarrow (ha,\delta )}=\left[{\begin{array}{ccc}\cos(LST)&\sin(LST)&0\\\sin(LST)&-\cos(LST)&0\\0&0&1\\\end{array}}\;\right]\;.\end{aligned}}\,\!}

2.2 (Ha,Dec) to (Azimuth, Altitude)

Because ${\displaystyle (ha,\delta )}$ are Earth-based coordinates, this conversion depends only on your terrestrial latitude ${\displaystyle \phi }$:

{\displaystyle {\begin{aligned}{R}_{(ha,\delta )\rightarrow (az,alt)}=\left[{\begin{array}{ccc}-\sin \phi &0&\cos \phi \\0&-1&0\\\cos \phi &0&\sin \phi \\\end{array}}\;\right]\;.\end{aligned}}\,\!}

2.3 Equatorial to Galactic

{\displaystyle {\begin{aligned}{R}_{(\alpha ,\delta )_{1950}\rightarrow (\ell ,b)}=\left[{\begin{array}{rrr}-0.066989&-0.872756&-0.483539\\0.492728&-0.450347&0.744585\\-0.867601&-0.188375&0.460200\\\end{array}}\;\right]\;.\end{aligned}}\,\!}

In truth, the precession of the equatorial coordinate system makes this matrix a function of time: the equatorial coordinates move around the sky, but the Galactic ones do not. Precession amounts to nearly an arcminute per year! In principle, you should derive the matrix for the current epoch. In practice, you may not need such high accuracy. Better than the 1950 version are the numbers for epoch 2000, which are in Green’s equation (14.55); epoch 2000 is lots closer to the present than is epoch 1950:

{\displaystyle {\begin{aligned}{R}_{(\alpha ,\delta )_{2000}\rightarrow (\ell ,b)}=\left[{\begin{array}{rrr}-0.054876&-0.873437&-0.483835\\0.494109&-0.444830&0.746982\\-0.867666&-0.198076&0.455984\\\end{array}}\;\right]\;.\end{aligned}}\,\!}

2.4 Precession—converting equatorial between epochs.

We won’t need these for the lab course, but we give you the info for the sake of completeness. Generating the rotation matrix for precession is a bit tedious and we won’t give the explicit formulae here. They are in Green’s book. The elements of the matrix are in equation (9.31). These elements contain angles, which depend on time as in equation (9.23) if you are converting from epoch 2000 to some other epoch. Precession isn’t all there is; for precision exceeding $\displaystyle \sim 10”$ you also need to account for nutation of the Earth, which has a random component and is not completely predictable. For the complete story, see Green’s chapter 9 and The Astronomical Almanac 1998, pages B39-B43—for interested parties only!

3 Doing all this in Python

Obviously, all this stuff is simple in Python/Numpy, which deals easily with matrices. Before beginning, though, a cautionary note about 2-D arrays on computers:

3.1 Examples in Python

Test ${\displaystyle {R}_{(ha,\delta )\rightarrow (az,alt)}}$, going both forwards and backwards. For an observatory at latitude ${\displaystyle 41.36^{\circ }}$, ${\displaystyle (az,alt)=(137.60^{\circ },32.43^{\circ })}$ transforms to ${\displaystyle (ha,\delta )=(325.05s^{\circ },-6.52^{\circ })}$. Hour angle is usually given in hours using sexagesimal notation: ${\displaystyle ha=21^{h}40^{m}12^{s}}$.) See the PyEphem module for functions that go back and forth between decimal and sexagesimal notations.

Test ${\displaystyle {R}_{(\alpha ,\delta )\rightarrow (ha,\delta )}}$ by making up your own example, using the fact that ${\displaystyle ha=LST-\alpha }$.

Test ${\displaystyle {R}_{(\alpha ,\delta )_{1950}\rightarrow (\ell ,b)}}$ for the Crab Nebula. The Crab has 1950 equatorial coordinates $\displaystyle (\alpha , \delta )=(05^ h31^ m.5, 21^\circ 59’)$ and Galactic coordinates $\displaystyle (\ell ,b)=(184^\circ 33’, -5^\circ 47’)$ .

Finally, put them all together and make sure that works, too. For all of these, make sure you know how to go backwards! For example, suppose you want to convert ${\displaystyle (az,alt)\rightarrow (\alpha ,\delta )}$. You need to first apply ${\displaystyle {R}_{(ha,\delta )\rightarrow (az,alt)}^{-1}}$ and then ${\displaystyle {R}_{(\alpha ,\delta )\rightarrow (ha,\delta )}^{-1}}$. So the full rotation matrix in this case is…

${\displaystyle {R}_{(az,alt)\rightarrow (\alpha ,\delta )}={R}_{(\alpha ,\delta )\rightarrow (ha,\delta )}^{-1}\cdot {R}_{(ha,\delta )\rightarrow (az,alt)}^{-1}\,\!}$

Again, note the order! Applying ${\displaystyle {R}_{(\alpha ,\delta )\rightarrow (ha,\delta )}^{-1}}$ at the beginning in the matrix product means that it operates last on the vector ${\displaystyle {\vec {x}}}$.

[1] If you don’t believe this, check on the ${\displaystyle {R}}$’s below in §2!