Moving polygons around the sky

My friend Sam is part of a team building a Gravitational Wave Treasure Map. Due to the way gravitational wave telescopes (e.g., LIGO) work, gravitational wave events can only be roughly localized on the sky. For example, gravitational wave event S190923y could have come from anywhere in the dark red regions here:

As soon as a detection is made, alerts are sent out to the astronomical community. Teams of astronomers around the world then search the indicated region for new points of light, in hopes of finding an optical (electromagnetic) counterpart. They report back which regions of the sky they imaged, and the treasure map shows where to focus further searches.

In the course of developing this map, he ran in to a classic coordinate singularity problem, and I helped him work out the math behind it.

The singularity¶

Astronomers locate objects in angular coordinates, right ascension (RA) and declination (Dec)—analogous to longitude and latitude respectively. The problem is that, when your declination goes to $\pm 90^\circ$, your right ascension is undefined. (Just like when you're standing at the North Pole, you don't have a particular longitude.)

So, say you have a telescope that images a $10^\circ \times 10^\circ$ field of view. That's 10 degrees from the telescope pointing location, not 10 degrees RA by 10 degrees Dec. It's easy to see why: define a square centered at $(0^\circ, 0^\circ)$:

In [1]:
%matplotlib inline
%config InlineBackend.figure_formats = ['retina']
%config InlineBackend.print_figure_kwargs = {'facecolor': (1.0, 1.0, 1.0, 0.0)}
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

matplotlib.rcParams.update({
'image.origin': 'lower',
'image.interpolation': 'nearest',
'image.cmap': 'magma',
'font.family': 'serif',
'axes.grid': True,
'figure.figsize': (7, 7)
})

In [2]:
footprint_zero_center = [
[-5, 5],
[5, 5],
[5, -5],
[-5, -5],
]
footprint_zero_center_ra = np.asarray([
pt[0] for pt in footprint_zero_center
])
footprint_zero_center_dec = np.asarray([
pt[1] for pt in footprint_zero_center
])
ax = plt.subplot(111, projection='aitoff')
ax.scatter(
)
plt.tight_layout()


And try to move it up to the pole:

In [3]:
ax = plt.subplot(111, projection='aitoff')
ax.scatter(
)
ax.scatter(
)
ax.scatter(
)
ax.scatter(
)
plt.tight_layout()


If you have centered on $(0^\circ, 85^\circ)$, you'll have a corner at $(-5^\circ, 90^\circ)$ and $(5^\circ, 90^\circ)$. But, of course, RA is undefined when Dec is $90^\circ$! In other words, the two top corners end up occupying the same point, and your footprint is clearly not $10^\circ$ by $10^\circ$ any more.

My friend encountered this when he tried to place detector footprint outlines at the pole for a visualization. There were division by zero errors, wild distortions... basically, it didn't work.

It turns out there's no simple formula to move your polygon around in angular coordinates while preserving surface area (solid angle).

Transforming into a higher dimensional space¶

Actually, that's not entirely true. There's no simple 2-D formula. If you convert into 3D Cartesian coordinates, you can apply standard rotation matrices to move your footprint where-ever you like.

If we have $\text{RA} \in [0^\circ,\ 360^\circ]$ and $\text{Dec} \in [-90^\circ,\ 90^\circ]$, we can convert to Cartesian coordinates by first defining

\begin{align*} \phi &= (90^\circ - \text{Dec}) \times \frac{\pi}{180^\circ}\\ \theta &= \text{RA} \times \frac{\pi}{180^\circ}\\ \end{align*}

Then using the usual conventions, our coordinates are

\begin{align*} x &= r \cos \theta \sin \phi \\ y &= r \sin \theta \sin \phi \\ z &= r \cos \phi \end{align*}

but we only care about the surface of the unit sphere, so we always have $r = 1$. Let's define the functions to convert back and forth:

In [4]:
def ra_dec_to_uvec(ra, dec):
x = np.cos(theta) * np.sin(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(phi)
return x, y, z


With the conventions we've chosen, $(\text{RA},\ \text{Dec}) = (0^\circ,\ 0^\circ)$ should be along the X axis. Check and see:

In [5]:
ra_dec_to_uvec(0, 0)

Out[5]:
(1.0, 0.0, 6.123233995736766e-17)

Slight imprecisions in the conversion of the Dec coordinate led to a negligible nonzero $z$ component, but otherwise everything seems correct.

Verify that $\text{Dec}\,=\,90^\circ$ is along the Z axis:

In [6]:
ra_dec_to_uvec(0, 90)

Out[6]:
(0.0, 0.0, 1.0)

$\text{RA}\,=\,90^\circ$ is along Y.

In [7]:
ra_dec_to_uvec(90, 0)

Out[7]:
(6.123233995736766e-17, 1.0, 6.123233995736766e-17)

Let's go the other way:

In [8]:
def uvec_to_ra_dec(x, y, z):
r = np.sqrt(x**2 + y ** 2 + z ** 2)
theta = np.arctan2(y, x)
phi = np.arccos(z)
if theta < 0:
else:
return ra, dec


Two things to note: we're using the arctan2 function, which produces $\theta$ with the correct sign. We then handle negative angles specially to ensure $\text{RA} \in [0^\circ, 360^\circ]$. Check a few cases:

In [9]:
uvec_to_ra_dec(1.0, 0.0, 0.0)

Out[9]:
(0.0, 0.0)
In [10]:
uvec_to_ra_dec(*ra_dec_to_uvec(90, 0))

Out[10]:
(90.0, 0.0)
In [11]:
uvec_to_ra_dec(*ra_dec_to_uvec(0, 90))

Out[11]:
(0.0, 90.0)

Convert the whole footprint:

In [12]:
footprint_zero_center_uvec = ra_dec_to_uvec(
footprint_zero_center_ra,
footprint_zero_center_dec
)
(
footprint_zero_center_x,
footprint_zero_center_y,
footprint_zero_center_z
) = footprint_zero_center_uvec


And plot it in 3D:

In [13]:
from mpl_toolkits.mplot3d import Axes3D
def make_3d_axes(rows=1, cols=1, loc=1):
ax = plt.gcf().add_subplot(rows, cols, loc, projection='3d')
ax.set_facecolor([0, 0, 0, 0])
minval, maxval = -1.1, 1.1
ax.set_xlim(minval, maxval)
ax.set_ylim(minval, maxval)
ax.set_zlim(minval, maxval)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
return ax

def vec_from_origin(x, y, z):
return [0, x], [0, y], [0, z]

fig = plt.figure()
ax = make_3d_axes()
xs, ys, zs = [], [], []
for idx in range(footprint_zero_center_x.shape[0]):
ax.plot(*vec_from_origin(
footprint_zero_center_x[idx],
footprint_zero_center_y[idx],
footprint_zero_center_z[idx]
))
xs.append(footprint_zero_center_x[idx])
ys.append(footprint_zero_center_y[idx])
zs.append(footprint_zero_center_z[idx])

# close the square
xs.append(xs[0])
ys.append(ys[0])
zs.append(zs[0])
ax.plot(xs, ys, zs)
plt.tight_layout()