Open in Colab: https://colab.research.google.com/github/casangi/cngi_prototype/blob/master/docs/coordinates.ipynb


Open In Colab

Coordinate Conversions

The goal of this notebook is to demonstrate conversion of visibility and image data converted between spectral, spatial, and temporal reference frames using the next generation CASA infrastructure.

Ideally, we would treat our units as a property of the data by wrapping dask arrays with pint quantities, optionally using a custom-defined context (for domain-specific units e.g., Jy/Beam). For details, see relevant discussion and new standalone package for ongoing public concept demonstration efforts.

That is not sufficiently mature or supported yet, so the alternative is to use existing conversion routines that operate on numpy ndarrays and parallelize across dask chunks using mapping functions or custom delayed calls.

Note: work in progress

this section of the documentation is still under development and the code cells in this notebook will probably all throw some form of Exception if executed.

The following representative function definition should produce parallelized calls to a topo_to_lsrk (or analagous conversion) routine with array size corresponding to the chunk size of (time,baseline,pol), for all channels.

def topo_to_lsrk(input_array)
    # call SkyCoord and other methods
    output_array = function(input_array)
    return output_array(axis=-1, keepdims=True)

Then we would call it like

test = xarray.apply_ufunc(topo_to_lsrk, xds.DATA.chunk({'chan':-1}), input_core_dims=[['chan']], dask='parallelized', output_dtypes=[xds.DATA.dtype])
output = test.compute()

This approach requires finding or creating a function that operates on numpy arrays instead of the astropy ndarray subclass or custom objects.

Frequency Reference Frames

The previous approach, following official demos and other example notebooks from astropy docs using SkyCoord objects, yields

TypeError: Position information stored on coordinate frame is insufficient to do a full-space position transformation (representation class: <class 'astropy.coordinates.representation.UnitSphericalRepresentation'>)
[ ]:
!pip install --quiet astropy --upgrade
import astropy
astropy.__version__
'4.0.1.post1'

Unless this is >= v4.1 we won’t be able to use SpectralCoord class, but here is how it could look

[ ]:
import numpy as np
from astropy import units as u
from astropy.coordinates import Angle
from astropy.coordinates import EarthLocation
from astropy.coordinates import SkyCoord
from astropy.coordinates import SpectralCoord
from astropy.time import Time

We’ll need to specify the absolute position of observer and target to make use of the frame transformation methods. Luckily this information is usually carried with the ASDM/MS/image and if it isn’t, there are ways to access reasonable defaults.

[ ]:
try:
  vis_xds = cngi.dio.read_vis('some_vis.zarr', ddi=0)
  global_xds = cngi.dio.read_vis('some_vis.zarr', ddi='global')
  # these are just hacky ways to get a position for now and might now even work
  location = global_xds['ASDM_POSITION'].mean().values
  location = global_xds['ANT_POSITION'].mean().values
except:
  location = EarthLocation.of_site(global_xds['OBS_TELESCOPE_NAME'])

A spectral reference requires time information in addition to spatial coordinates. Some string of the form '2019-04-24T02:32:10' is required, and to be in the TOPO frame (e.g., ALMA), this should be the start of the observation.

[ ]:
try:
  observatory = location.get_itrs(obstime=vis_xds.time[0])
except:
  # if that fails maybe can convert from global properties
  import datetime
  time = global_xds['ASDM_startValidTime'][0]
  time = global_xds.OBS_TIME_RANGE.values[0]
  observatory = location.get_itrs(obstime=datetime.fromtimestamp(time))

Not only is the reference location of the “observer” required, the source properties must also be defined in the same frame of a form like SkyCoord('04h21m59.43s +19d32m06.4', frame='icrs', radial_velocity=23.9 * u.km / u.s, distance=144.321 * u.pc).

TODO: assign coordinate variables to d1, d2, d3 in vis.zarr

[ ]:
# assuming these direction axes are sensibly and consistently defined in radians
# either way, this will be far easier for images
radec_string = Angle(global_xds.SRC_DIRECTION * u.rad,
                     global_xds.SRC_DIRECTION * u.rad)

# telescope dependent... is this is kept anywhere in the ASDM/MS/image formats?
target_frame = 'FK5'

# another assumption lacking coordinate associated with the d2 dimension
recession = global_xds.SRC_PROPER_MOTION * u.Hz

source = SkyCoord(radec_string, frame=target_frame, radial_velocity=recession)
[ ]:
# initialize SpectralCoord instance from input xarray.Dataset

# compute for return
new_sc = sc.with_observer_stationary_relative_to('lsrk')

astropy adopts a standard definition of the LSR following Schönrich et al. 2010.

CASA supports a number of spectral reference frames and the definition of LSR is encoded in casacore and seems to follow the references published by Frank Ghigo at GBO.

Spatial Coordinates

Currently it is assumed that all image coordinates have native units of radians, which should allow for consistent selection and conversion between different spatial quantities at the ngCASA and/or application layer.

[ ]:
from astropy import units as u
from astropy.coordinates import EarthLocation
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.coordinates import SpectralCoord
from astropy.coordinates import AltAz

These astropy modules allow us to manipulate quantities for unit-aware computation, including transformation between the reference frames in which data are represented.

[ ]:
VLA_lat = 34.1*u.degree
VLA_lon = -107.6*u.degree
VLA_alt = 2114.89*u.m
print((VLA_lat, VLA_lon, VLA_alt))

observing_location = EarthLocation(lat=VLA_lat,
                                   lon=VLA_lon,
                                   height=VLA_alt,
                                   ellipsoid='WGS84')
print(observing_location, type(observing_location))
(<Quantity 34.1 deg>, <Quantity -107.6 deg>, <Quantity 2114.89 m>)
(-1599173.52082635, -5041233.67723585, 3556822.79344969) m <class 'astropy.coordinates.earth.EarthLocation'>

Space and time are notoriously difficult to untangle, but since we can define locations, coordinates on the celestial sphere, and times, we can transform between these quantities.

[ ]:
reference_time = Time(['2019-10-4T00:00:00'],
                        format='isot',
                        scale='utc',
                        location=observing_location)
print(reference_time, type(reference_time))

phase_center = SkyCoord(ra='19h59m28.5s', dec='+40d44m01.5s', frame='fk5')
print(phase_center)
['2019-10-04T00:00:00.000'] <class 'astropy.time.core.Time'>
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (299.86875, 40.73375)>

Converting between units is then as straightforward as:

[ ]:
alt_az_frame = AltAz(location=observing_location, obstime=reference_time)
print(alt_az_frame)
new_frame = phase_center.transform_to(alt_az_frame)
print()
print("Compare",
      (new_frame.az.radian,
       new_frame.alt.radian),
      "and",
      (new_frame.az.degree,
       new_frame.alt.degree))
<AltAz Frame (obstime=['2019-10-04T00:00:00.000'], location=(-1599173.52082635, -5041233.67723585, 3556822.79344969) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron)>

Compare (array([1.1522305]), array([1.07222892])) and (array([66.01794453]), array([61.43419156]))

We should ensure we are making use of the higher-precision mode when we implement radial velocity (and possibly all) corrections.

It should also be possible to rely on 'axisunits' or similar keys in the image Dataset attributes to set up a spectral frame for eager computation.

[ ]:
sc = SpectralCoord(img_xds.DATA.chan,
                   unit=astropy.units(img_xds.attrs['rest_frequency'].split('')[:-1]),
                   observer=observatory,
                   target= source,
                   doppler_reference=img_xds.attrs['spectral__reference'],
                   doppler_convention=img_xds.attrs['velocity__type'])