Skip to content

Instantly share code, notes, and snippets.

@demorest
Last active June 29, 2020 16:34
Embed
What would you like to do?
Interferometric UVW calculation using astropy
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]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment