Interferometric UVW calculation using astropy

uvw_compare.py

import numpy as np
import astropy.coordinates as coord
import astropy.time
import astropy.units as u

# Can do this to get updated IERS B values into astropy
from astropy.utils import iers
from astropy.utils import data
iers_b = iers.IERS_B.open(data.download_file(iers.IERS_B_URL, cache=True))
iers_auto = iers.IERS_Auto.open()

# VLA B-config positions, ITRF, m
antpos = np.array([[-1604008.7431 , -5042135.8194 ,  3553403.7084 ],
       [-1601315.9011 , -5041985.30447,  3554808.3081 ],
       [-1604865.6579 , -5042190.0302 ,  3552962.3597 ],
       [-1601173.9801 , -5041902.6559 ,  3554987.5267 ],
       [-1596127.7286 , -5045193.7343 ,  3552652.4174 ],
       [-1600863.6816 , -5039885.303  ,  3557965.315  ],
       [-1601068.7948 , -5042051.9181 ,  3554824.8427 ],
       [-1602592.8713 , -5042054.9913 ,  3554140.7111 ],
       [-1601061.9634 , -5041175.8811 ,  3556058.0394 ],
       [-1600801.9175 , -5042219.3706 ,  3554706.4492 ],
       [-1600781.0408 , -5039347.416  ,  3558761.5242 ],
       [-1599926.1065 , -5042772.9632 ,  3554319.8098 ],
       [-1603249.6777 , -5042091.4126 ,  3553797.8106 ],
       [-1597899.905  , -5044068.6843 ,  3553432.4423 ],
       [-1601110.0378 , -5041488.073  ,  3555597.4397 ],
       [-1599340.8237 , -5043150.9642 ,  3554065.1933 ],
       [-1600690.6036 , -5038758.7058 ,  3559632.0653 ],
       [-1601034.4102 , -5040996.5142 ,  3556322.916  ],
       [-1600930.0767 , -5040316.4046 ,  3557330.4031 ],
       [-1600416.512  , -5042462.443  ,  3554536.0387 ],
       [-1602044.8966 , -5042025.8044 ,  3554427.8217 ],
       [-1598663.0898 , -5043581.3793 ,  3553767.0201 ],
       [-1605808.6423 , -5042230.088  ,  3552459.2035 ],
       [-1601147.9459 , -5041733.8322 ,  3555235.9444 ],
       [-1606841.9757 , -5042279.6672 ,  3551913.0239 ],
       [-1597053.13   , -5044604.6754 ,  3553058.9804 ],
       [-1601614.0825 , -5042001.6537 ,  3554652.508  ]]) * u.m

# Time of observation:
t = astropy.time.Time(57000.12345, format='mjd', scale='utc')

# Format antenna positions and VLA center as EarthLocation.
antpos_ap = coord.EarthLocation(x=antpos[:,0], y=antpos[:,1], z=antpos[:,2])
vla = coord.EarthLocation.of_site('vla')

# Convert antenna pos terrestrial to celestial.  For astropy use 
# get_gcrs_posvel(t)[0] rather than get_gcrs(t) because if a velocity 
# is attached to the coordinate astropy will not allow us to do additional 
# transformations with it (https://github.com/astropy/astropy/issues/6280)
vla_p, vla_v = vla.get_gcrs_posvel(t)
antpos_c_ap = coord.GCRS(antpos_ap.get_gcrs_posvel(t)[0], 
        obstime=t, obsgeoloc=vla_p, obsgeovel=vla_v)

# Define the UVW frame relative to a certain point on the sky.  There are
# two versions, depending on whether the sky offset is done in ICRS 
# or GCRS:
pnt = coord.SkyCoord('12h34m56.7s', '65d43m21.0s', frame='icrs')
#frame_uvw = pnt.skyoffset_frame() # ICRS
frame_uvw = pnt.transform_to(antpos_c_ap).skyoffset_frame() # GCRS

# Rotate antenna positions into UVW frame.
antpos_uvw_ap = antpos_c_ap.transform_to(frame_uvw)

# Full set of baselines would be differences between all pairs of 
# antenna positions, but we'll just do relative to the first antenna
# for simplicity.
bl_uvw_ap = antpos_uvw_ap.cartesian - antpos_uvw_ap.cartesian[0]

# SkyOffsetFrame coords seem to come out as WUV, so shuffle into UVW order:
bl_uvw_ap = coord.CartesianRepresentation(
        x = bl_uvw_ap.y,
        y = bl_uvw_ap.z,
        z = bl_uvw_ap.x)

# Rest of script does the same(?) thing in CASA for comparison.
# Requires casatools from CASA 6 which can be installed via:
# pip install --index-url https://casa-pip.nrao.edu/repository/pypi-casa-release/simple casatools==6.0.0.27
try:
    import casatools
except ImportError:
    casatools = None

if casatools is not None:

    def casa_to_astropy(c):
        """Convert CASA spherical coords to astropy CartesianRepresentation"""
        sph = coord.SphericalRepresentation(
                lon=c['m0']['value']*u.Unit(c['m0']['unit']), 
                lat=c['m1']['value']*u.Unit(c['m1']['unit']), 
                distance=c['m2']['value']*u.Unit(c['m2']['unit']))
        return sph.represent_as(coord.CartesianRepresentation)

    # The necessary interfaces:
    me = casatools.measures()
    qa = casatools.quanta()
    qq = qa.quantity

    # Init CASA frame info:
    me.doframe(me.observatory('VLA'))
    me.doframe(me.epoch('UTC',qq(t.mjd,'d')))
    me.doframe(me.direction('J2000',
        qq(pnt.ra.to(u.rad).value, 'rad'),
        qq(pnt.dec.to(u.rad).value, 'rad')))

    # Format antenna positions for CASA:
    antpos_casa = me.position('ITRF',
            qq(antpos[:,0].to(u.m).value,'m'),
            qq(antpos[:,1].to(u.m).value,'m'),
            qq(antpos[:,2].to(u.m).value,'m'))

    # Converts from ITRF to "J2000":
    antpos_c_casa = me.asbaseline(antpos_casa)

    # Rotate into UVW frame
    antpos_uvw_casa = me.touvw(antpos_c_casa)[0]

    # me.expand would compute all pairs of baselines but here we convert
    # to astropy CartesianRepresentation, and only do baselines to the
    # first antenna for easier comparison
    bl_uvw_casa = casa_to_astropy(antpos_uvw_casa)
    bl_uvw_casa -= bl_uvw_casa[0]