diff --git a/africanus/rime/parangles.py b/africanus/rime/parangles.py index bcd056706..ec184c409 100644 --- a/africanus/rime/parangles.py +++ b/africanus/rime/parangles.py @@ -4,92 +4,23 @@ 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): @@ -97,17 +28,6 @@ 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` @@ -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. @@ -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) diff --git a/africanus/rime/parangles_astropy.py b/africanus/rime/parangles_astropy.py new file mode 100644 index 000000000..b92f9a3de --- /dev/null +++ b/africanus/rime/parangles_astropy.py @@ -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) diff --git a/africanus/rime/parangles_casa.py b/africanus/rime/parangles_casa.py new file mode 100644 index 000000000..181fd5896 --- /dev/null +++ b/africanus/rime/parangles_casa.py @@ -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]) diff --git a/africanus/rime/tests/conftest.py b/africanus/rime/tests/conftest.py new file mode 100644 index 000000000..576635e7d --- /dev/null +++ b/africanus/rime/tests/conftest.py @@ -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) diff --git a/africanus/rime/tests/test_parangles.py b/africanus/rime/tests/test_parangles.py new file mode 100644 index 000000000..e35e1565c --- /dev/null +++ b/africanus/rime/tests/test_parangles.py @@ -0,0 +1,157 @@ +import numpy as np +import pytest + +from africanus.rime.parangles import _discovered_backends + +no_casa = 'casa' not in _discovered_backends +no_astropy = 'astropy' not in _discovered_backends + + +def _julian_day(year, month, day): + """ + Given a Anno Dominei date, computes the Julian Date in days. + + Parameters + ---------- + year : int + month : int + day : int + + Returns + ------- + float + Julian Date + """ + + # Formula below from + # http://scienceworld.wolfram.com/astronomy/JulianDate.html + # Also agrees with https://gist.github.com/jiffyclub/1294443 + return (367*year - int(7*(year + int((month+9)/12))/4) + - int((3*(int(year + (month - 9)/7)/100)+1)/4) + + int(275*month/9) + day + 1721028.5) + + +def _modified_julian_date(year, month, day): + """ + Given a Anno Dominei date, computes the Modified Julian Date in days. + + Parameters + ---------- + year : int + month : int + day : int + + Returns + ------- + float + Modified Julian Date + """ + + return _julian_day(year, month, day) - 2400000.5 + + +def _observation_endpoints(year, month, date, hour_duration): + """ + Start and end points of an observation starting on + ``year-month-day`` and of duration ``hour_duration`` + in Modified Julian Date seconds + """ + start = _modified_julian_date(year, month, date) + end = start + hour_duration / 24. + + # Convert to seconds + start *= 86400. + end *= 86400. + + return (start, end) + + +@pytest.mark.parametrize('backend', [ + 'test', + pytest.param('casa', marks=pytest.mark.skipif( + no_casa, + reason='python-casascore not installed')), + pytest.param('astropy', marks=pytest.mark.skipif( + no_astropy, + reason="astropy not installed"))]) +@pytest.mark.parametrize('observation', [(2018, 1, 1, 4)]) +def test_parallactic_angles(observation, wsrt_ants, backend): + import numpy as np + from africanus.rime import parallactic_angles + + start, end = _observation_endpoints(*observation) + time = np.linspace(start, end, 5) + ant = wsrt_ants[:4, :] + fc = np.random.random((2,)).astype(np.float64) + + pa = parallactic_angles(time, ant, fc, backend=backend) + assert pa.shape == (5, 4) + + +@pytest.mark.skipif(no_casa or no_astropy, + reason="Neither python-casacore or astropy installed") +# Parametrize on observation length and error tolerance +@pytest.mark.parametrize('obs_and_tol', [ + ((2018, 1, 1, 4), "10s"), + ((2018, 2, 20, 8), "10s"), + ((2018, 11, 2, 4), "10s")]) +def test_compare_astropy_and_casa(obs_and_tol, wsrt_ants): + """ + Compare astropy and python-casacore parallactic angle implementations. + More work needs to be done here to get things lined up closer, + but the tolerances above suggest nothing > 10 arcseconds. + """ + import numpy as np + from africanus.rime.parangles_casa import casa_parallactic_angles + from africanus.rime.parangles_astropy import astropy_parallactic_angles + from astropy import units + from astropy.coordinates import Angle + + obs, rtol = obs_and_tol + start, end = _observation_endpoints(*obs) + + time = np.linspace(start, end, 5) + ant = wsrt_ants[:4, :] + fc = np.array([0., 1.04719755], dtype=np.float64) + + astro_pa = astropy_parallactic_angles(time, ant, fc) + casa_pa = casa_parallactic_angles(time, ant, fc, zenith_frame='AZELGEO') + + # Convert to angle degrees + astro_pa = Angle(astro_pa, unit=units.deg).wrap_at(180*units.deg) + casa_pa = Angle(casa_pa*units.rad, unit=units.deg).wrap_at(180*units.deg) + + # Difference in degrees, wrapped at 180 + diff = np.abs((astro_pa - casa_pa).wrap_at(180*units.deg)) + assert np.all(np.abs(diff) < Angle(rtol)) + + +@pytest.mark.parametrize('backend', [ + 'test', + pytest.param('casa', marks=pytest.mark.skipif( + no_casa, + reason='python-casascore not installed')), + pytest.param('astropy', marks=pytest.mark.skipif( + no_astropy, + reason="astropy not installed"))]) +@pytest.mark.parametrize('observation', [(2018, 1, 1, 4)]) +def test_dask_parallactic_angles(observation, wsrt_ants, backend): + da = pytest.importorskip('dask.array') + from africanus.rime import parallactic_angles as np_parangle + from africanus.rime.dask import parallactic_angles as da_parangle + + start, end = _observation_endpoints(*observation) + np_times = np.linspace(start, end, 5) + np_ants = wsrt_ants[:4, :] + np_fc = np.random.random(size=2) + + np_pa = np_parangle(np_times, np_ants, np_fc, backend=backend) + np_pa = np.asarray(np_pa) + + da_times = da.from_array(np_times, chunks=(2, 3)) + da_ants = da.from_array(np_ants, chunks=((2, 2), 3)) + da_fc = da.from_array(np_fc, chunks=2) + + da_pa = da_parangle(da_times, da_ants, da_fc, backend=backend) + + assert np.all(np_pa == da_pa.compute()) diff --git a/africanus/rime/tests/test_rime.py b/africanus/rime/tests/test_rime.py index 79dfd85b1..141e3fb69 100644 --- a/africanus/rime/tests/test_rime.py +++ b/africanus/rime/tests/test_rime.py @@ -67,149 +67,6 @@ def test_feed_rotation(): assert np.allclose(fr, np_expr.reshape(10, 5, 2, 2)) -def _julian_day(year, month, day): - """ - Given a Anno Dominei date, computes the Julian Date in days. - - Parameters - ---------- - year : int - month : int - day : int - - Returns - ------- - float - Julian Date - """ - - # Formula below from - # http://scienceworld.wolfram.com/astronomy/JulianDate.html - # Also agrees with https://gist.github.com/jiffyclub/1294443 - return (367*year - int(7*(year + int((month+9)/12))/4) - - int((3*(int(year + (month - 9)/7)/100)+1)/4) - + int(275*month/9) + day + 1721028.5) - - -def _modified_julian_date(year, month, day): - """ - Given a Anno Dominei date, computes the Modified Julian Date in days. - - Parameters - ---------- - year : int - month : int - day : int - - Returns - ------- - float - Modified Julian Date - """ - - return _julian_day(year, month, day) - 2400000.5 - - -def _observation_endpoints(year, month, date, hour_duration): - """ - Start and end points of an observation starting on - ``year-month-day`` and of duration ``hour_duration`` - in Modified Julian Date seconds - """ - start = _modified_julian_date(year, month, date) - end = start + hour_duration / 24. - - # Convert to seconds - start *= 86400. - end *= 86400. - - return (start, end) - - -@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) - - -no_casa = 'casa' not in _discovered_backends -no_astropy = 'astropy' not in _discovered_backends - - -@pytest.mark.parametrize('backend', [ - 'test', - pytest.param('casa', marks=pytest.mark.skipif( - no_casa, - reason='python-casascore not installed')), - pytest.param('astropy', marks=pytest.mark.skipif( - no_astropy, - reason="astropy not installed"))]) -@pytest.mark.parametrize('observation', [(2018, 1, 1, 4)]) -def test_parallactic_angles(observation, wsrt_ants, backend): - import numpy as np - from africanus.rime import parallactic_angles - - start, end = _observation_endpoints(*observation) - time = np.linspace(start, end, 5) - ant = wsrt_ants[:4, :] - fc = np.random.random((2,)).astype(np.float64) - - pa = parallactic_angles(time, ant, fc, backend=backend) - assert pa.shape == (5, 4) - - -@pytest.mark.skipif(no_casa or no_astropy, - reason="Neither python-casacore or astropy installed") -# Parametrize on observation length and error tolerance -@pytest.mark.parametrize('obs_and_tol', [ - ((2018, 1, 1, 4), "10s"), - ((2018, 2, 20, 8), "10s"), - ((2018, 11, 2, 4), "10s")]) -def test_compare_astropy_and_casa(obs_and_tol, wsrt_ants): - """ - Compare astropy and python-casacore parallactic angle implementations. - More work needs to be done here to get things lined up closer, - but the tolerances above suggest nothing > 10 arcseconds. - """ - import numpy as np - from africanus.rime import parallactic_angles - from astropy import units - from astropy.coordinates import Angle - - obs, rtol = obs_and_tol - start, end = _observation_endpoints(*obs) - - time = np.linspace(start, end, 5) - ant = wsrt_ants[:4, :] - fc = np.array([0., 1.04719755], dtype=np.float64) - - astro_pa = parallactic_angles(time, ant, fc, backend='astropy') - casa_pa = parallactic_angles(time, ant, fc, backend='casa') - - # Convert to angle degrees - astro_pa = Angle(astro_pa, unit=units.deg).wrap_at(180*units.deg) - casa_pa = Angle(casa_pa*units.rad, unit=units.deg).wrap_at(180*units.deg) - - # Difference in degrees, wrapped at 180 - diff = np.abs((astro_pa - casa_pa).wrap_at(180*units.deg)) - assert np.all(np.abs(diff) < Angle(rtol)) - - def test_jones_2x2_mul(): import numpy as np from africanus.rime.predict import jones_2x2_mul @@ -295,37 +152,6 @@ def test_dask_phase_delay(): assert np.all(np_phase == dask_phase) -@pytest.mark.parametrize('backend', [ - 'test', - pytest.param('casa', marks=pytest.mark.skipif( - no_casa, - reason='python-casascore not installed')), - pytest.param('astropy', marks=pytest.mark.skipif( - no_astropy, - reason="astropy not installed"))]) -@pytest.mark.parametrize('observation', [(2018, 1, 1, 4)]) -def test_dask_parallactic_angles(observation, wsrt_ants, backend): - da = pytest.importorskip('dask.array') - from africanus.rime import parallactic_angles as np_parangle - from africanus.rime.dask import parallactic_angles as da_parangle - - start, end = _observation_endpoints(*observation) - np_times = np.linspace(start, end, 5) - np_ants = wsrt_ants[:4, :] - np_fc = np.random.random(size=2) - - np_pa = np_parangle(np_times, np_ants, np_fc, backend=backend) - np_pa = np.asarray(np_pa) - - da_times = da.from_array(np_times, chunks=(2, 3)) - da_ants = da.from_array(np_ants, chunks=((2, 2), 3)) - da_fc = da.from_array(np_fc, chunks=2) - - da_pa = da_parangle(da_times, da_ants, da_fc, backend=backend) - - assert np.all(np_pa == da_pa.compute()) - - def test_dask_feed_rotation(): da = pytest.importorskip('dask.array') import numpy as np