Skip to content

Commit

Permalink
Make CASA and 'AZEL' frame the parangle standard.
Browse files Browse the repository at this point in the history
astropy parangles are out, see here:
#21

* Split parallactic angle test cases into separate files.
  • Loading branch information
sjperkins committed Oct 22, 2018
1 parent b38465a commit e71e90b
Show file tree
Hide file tree
Showing 6 changed files with 303 additions and 267 deletions.
112 changes: 19 additions & 93 deletions africanus/rime/parangles.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,110 +4,30 @@
from __future__ import division
from __future__ import print_function

import numpy as np
import warnings

from ..util.requirements import requires_optional
from .parangles_astropy import (have_astropy_parangles,
astropy_parallactic_angles)
from .parangles_casa import (have_casa_parangles,
casa_parallactic_angles)

_discovered_backends = []
_standard_backends = ['casa', 'astropy']

try:
import pyrap.measures
import pyrap.quanta as pq
except ImportError:
pass
else:
_discovered_backends.append('casa')
# Create a measures server
meas_serv = pyrap.measures.measures()

try:
from astropy.coordinates import (EarthLocation, SkyCoord,
AltAz, CIRS)
from astropy.time import Time
from astropy import units
except ImportError:
pass
else:
_discovered_backends.append('astropy')


@requires_optional('pyrap.measures', 'pyrap.quanta')
def casa_parallactic_angles(times, antenna_positions, field_centre):
"""
Computes parallactic angles per timestep for the given
reference antenna position and field centre.
"""

# Create direction measure for the zenith
zenith = meas_serv.direction('AZELGEO', '0deg', '90deg')

# Create position measures for each antenna
reference_positions = [meas_serv.position(
'itrf',
*(pq.quantity(x, 'm') for x in pos))
for pos in antenna_positions]

# Compute field centre in radians
fc_rad = meas_serv.direction('J2000', *(pq.quantity(f, 'rad')
for f in field_centre))

return np.asarray([
# Set current time as the reference frame
meas_serv.do_frame(meas_serv.epoch("UTC", pq.quantity(t, "s")))
and
[ # Set antenna position as the reference frame
meas_serv.do_frame(rp)
and
meas_serv.posangle(fc_rad, zenith).get_value("rad")
for rp in reference_positions
]
for t in times])


@requires_optional('astropy')
def astropy_parallactic_angles(times, antenna_positions, field_centre):
"""
Computes parallactic angles per timestep for the given
reference antenna position and field centre.
"""
ap = antenna_positions
fc = field_centre

# Convert from MJD second to MJD
times = Time(times / 86400.00, format='mjd', scale='utc')
if have_astropy_parangles:
_discovered_backends.append('astropy')

ap = EarthLocation.from_geocentric(
ap[:, 0], ap[:, 1], ap[:, 2], unit='m')
fc = SkyCoord(ra=fc[0], dec=fc[1], unit=units.rad, frame='fk5')
pole = SkyCoord(ra=0, dec=90, unit=units.deg, frame='fk5')
if have_casa_parangles:
_discovered_backends.append('casa')

cirs_frame = CIRS(obstime=times)
pole_cirs = pole.transform_to(cirs_frame)
fc_cirs = fc.transform_to(cirs_frame)

altaz_frame = AltAz(location=ap[None, :], obstime=times[:, None])
pole_altaz = pole_cirs[:, None].transform_to(altaz_frame)
fc_altaz = fc_cirs[:, None].transform_to(altaz_frame)
return fc_altaz.position_angle(pole_altaz)
_standard_backends = set(['casa', 'astropy', 'test'])


def parallactic_angles(times, antenna_positions, field_centre, **kwargs):
"""
Computes parallactic angles per timestep for the given
reference antenna position and field centre.
Notes
-----
* The python-casacore backend uses an ``AZELGEO`` frame in order
to more closely agree with the astropy backend,
but slightly differs from ``AZEL`` frame using in MeqTrees.
* The astropy backend is slightly more than 2x faster than
the casa backend
* The astropy and python-casacore differ by at most
10 arcseconds
Parameters
----------
times : :class:`numpy.ndarray`
Expand All @@ -118,12 +38,12 @@ def parallactic_angles(times, antenna_positions, field_centre, **kwargs):
in *metres* in the *ITRF* frame.
field_centre : :class:`numpy.ndarray`
Field centre of shape :code:`(2,)` in *radians*
backend : {'casa', 'astropy', 'test'}, optional
backend : {'casa', 'test'}, optional
Backend to use for calculating the parallactic angles.
* ``casa`` defers to an implementation
depending on ``python-casacore``
* ``astropy`` defers to astropy.
depending on ``python-casacore``.
This backend should be used by default.
* ``test`` creates parallactic angles
by multiplying the ``times`` and ``antenna_position``
arrays. It exist solely for testing.
Expand All @@ -143,11 +63,17 @@ def parallactic_angles(times, antenna_positions, field_centre, **kwargs):
raise ValueError("None of the standard backends "
"%s are installed" % _standard_backends)

if backend not in _standard_backends:
raise ValueError("'%s' is not one of the "
"standard backends '%s'"
% (backend, _standard_backends))

if not field_centre.shape == (2,):
raise ValueError("Invalid field_centre shape %s" %
(field_centre.shape,))

if backend == 'astropy':
warnings.warn('astropy backend currently returns the incorrect values')
return astropy_parallactic_angles(times,
antenna_positions,
field_centre)
Expand Down
44 changes: 44 additions & 0 deletions africanus/rime/parangles_astropy.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
# -*- coding: utf-8 -*-

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from ..util.requirements import requires_optional

try:
from astropy.coordinates import (EarthLocation, SkyCoord,
AltAz, CIRS)
from astropy.time import Time
from astropy import units
except ImportError:
have_astropy_parangles = False
else:
have_astropy_parangles = True


@requires_optional('astropy')
def astropy_parallactic_angles(times, antenna_positions, field_centre):
"""
Computes parallactic angles per timestep for the given
reference antenna position and field centre.
"""
ap = antenna_positions
fc = field_centre

# Convert from MJD second to MJD
times = Time(times / 86400.00, format='mjd', scale='utc')

ap = EarthLocation.from_geocentric(
ap[:, 0], ap[:, 1], ap[:, 2], unit='m')
fc = SkyCoord(ra=fc[0], dec=fc[1], unit=units.rad, frame='fk5')
pole = SkyCoord(ra=0, dec=90, unit=units.deg, frame='fk5')

cirs_frame = CIRS(obstime=times)
pole_cirs = pole.transform_to(cirs_frame)
fc_cirs = fc.transform_to(cirs_frame)

altaz_frame = AltAz(location=ap[None, :], obstime=times[:, None])
pole_altaz = pole_cirs[:, None].transform_to(altaz_frame)
fc_altaz = fc_cirs[:, None].transform_to(altaz_frame)
return fc_altaz.position_angle(pole_altaz)
54 changes: 54 additions & 0 deletions africanus/rime/parangles_casa.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
# -*- coding: utf-8 -*-

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np

from ..util.requirements import requires_optional

try:
import pyrap.measures
import pyrap.quanta as pq
except ImportError:
have_casa_parangles = False
else:
have_casa_parangles = True

# Create a measures server
meas_serv = pyrap.measures.measures()


@requires_optional('pyrap.measures', 'pyrap.quanta')
def casa_parallactic_angles(times, antenna_positions, field_centre,
zenith_frame='AZEL'):
"""
Computes parallactic angles per timestep for the given
reference antenna position and field centre.
"""

# Create direction measure for the zenith
zenith = meas_serv.direction(zenith_frame, '0deg', '90deg')

# Create position measures for each antenna
reference_positions = [meas_serv.position(
'itrf',
*(pq.quantity(x, 'm') for x in pos))
for pos in antenna_positions]

# Compute field centre in radians
fc_rad = meas_serv.direction('J2000', *(pq.quantity(f, 'rad')
for f in field_centre))

return np.asarray([
# Set current time as the reference frame
meas_serv.do_frame(meas_serv.epoch("UTC", pq.quantity(t, "s")))
and
[ # Set antenna position as the reference frame
meas_serv.do_frame(rp)
and
meas_serv.posangle(fc_rad, zenith).get_value("rad")
for rp in reference_positions
]
for t in times])
29 changes: 29 additions & 0 deletions africanus/rime/tests/conftest.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Tests for `codex-africanus` package."""


import numpy as np
import pytest


@pytest.fixture
def wsrt_ants():
""" Westerbork antenna positions """
return np.array([
[3828763.10544699, 442449.10566454, 5064923.00777],
[3828746.54957258, 442592.13950824, 5064923.00792],
[3828729.99081359, 442735.17696417, 5064923.00829],
[3828713.43109885, 442878.2118934, 5064923.00436],
[3828696.86994428, 443021.24917264, 5064923.00397],
[3828680.31391933, 443164.28596862, 5064923.00035],
[3828663.75159173, 443307.32138056, 5064923.00204],
[3828647.19342757, 443450.35604638, 5064923.0023],
[3828630.63486201, 443593.39226634, 5064922.99755],
[3828614.07606798, 443736.42941621, 5064923.],
[3828609.94224429, 443772.19450029, 5064922.99868],
[3828601.66208572, 443843.71178407, 5064922.99963],
[3828460.92418735, 445059.52053929, 5064922.99071],
[3828452.64716351, 445131.03744105, 5064922.98793]],
dtype=np.float64)
Loading

0 comments on commit e71e90b

Please sign in to comment.