Skip to content

Instantly share code, notes, and snippets.

@demorest
Last active Jun 29, 2020
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