Open in Colab: https://colab.research.google.com/github/casangi/examples/blob/master/community/ALMAephemimagingVenusMfsMosaic.ipynb

---

ALMA ephemeris object imaging: continuum (mosaic gridder - ephemeris data repository)

Use ephemeris data in the data repository

Original Author(s): ttsutsum@nrao.edu

Description

The following demo show an incorrect way to do continuum imaging of ALMA observation of Venus with the attached ephemeris table. For the known major solar system objects such as Venus, the ephemeris data (DE200/DE405) are avaiable in the data repository.

Data

The data, venus_ephem_test.ms is available in the casatestdata repository, under https://open-bitbucket.nrao.edu/projects/CASA/repos/casatestdata/browse/measurementset/alma/.

Or at NRAO, casatestdata is accessible at /home/casa/data/casatestdata

Installation

Install CASA

Skip this step if CASA is already installed

import os
print("Installing CASA ...")
os.system('pip install casaconfig')
os.system('pip install casatools==6.7.0.31')
os.system('pip install casatasks==6.7.0.31')

# for Colab
mydatapath = '/content/.casa/data'
#
#mydatapath='data'

import pathlib
from casaconfig import config
if not pathlib.Path(mydatapath).exists():
  pathlib.Path(mydatapath).mkdir(parents=True)
config.measurespath=mydatapath

Install additional modules used in this notebook

import os
os.system("pip install aplpy")
os.system("pip install ipympl")
print("complete")

This notebook uses functions from EphemerisObjectImagingDemoFunctions.py so that file must be present in the same directory running this notebook and need to import them

# if you haven't downloaded EphemerisObjectImagingDemoFuntions.py
!wget https://raw.githubusercontent.com/casangi/examples/refs/heads/CAS-13662/community/EphemerisObjectImagingDemoFunctions.py

from EphemerisObjectImagingDemoFunctions import *

Import required tools and tasks

from casatasks import getephemtable, tclean, exportfits, listobs, imstat, imhead
from casatools import measures, quanta, table, image
import os
from astropy.io import fits
from astropy.wcs import WCS
from matplotlib import pylab as pl
import numpy as np
import aplpy
me = measures()
qa = quanta()
tb = table()
ia = image()

Imaging using the default ephemeris data in the data repository

In order to track the ephemeris source during the observations, the proper ephemeris data is used. Ususally such ephemeris data table(s) are attached to the Measurement Set automatically during importadm when they are available. The example using the attached epehemeris table is described in another notebook (ALMAephemimagingVenusMfsMosaicTrackfield.ipynb).

Obtain data

msfile = 'venus_ephem_test.ms'

# Use the ALMA Venus data in casatestdata
# At NRAO, point to the casatestdata directory
#rootdatapath = '/home/casa/data/casatestdata'
# or set a proper path to the local copy of the casatestdata repository
# e.g.
# rootdatapath = '/Volumes/ssd1/casatestdata'
#inputdatapath=f'{rootdatapath}/measurementset/alma/'

For Colab or for downloading the input data locally execute next cell, otherwise skip it

# If running on Colab

# sparse-checkout of the data repository

# install git lfs
os.system('pip install git-lfs -q')
os.system('git-lfs install')

os.system('git clone --filter=blob:none --sparse https://open-bitbucket.nrao.edu/scm/casa/casatestdata.git mytestdata')
%cd mytestdata
os.system('git sparse-checkout set measurementset/alma/venus_ephem_test.ms')
%cd ..

Updated Git hooks.
Git LFS initialized.

Cloning into 'mytestdata'...
/Users/ttsutsum/SWDevel/casa_mod/mytestenv/lib/python3.10/site-packages/IPython/core/magics/osm.py:417: UserWarning: This is now an optional IPython functionality, setting dhist requires you to install the `pickleshare` library.
  self.shell.db['dhist'] = compress_dhist(dhist)[-100:]

/Users/ttsutsum/SWDevel/casa_mod/examples/community/mytestdata

0

# If running on Colab
inputdatapath = 'mytestdata/measurementset/alma/'

# Check the data exists
from casatasks import listobs
listobs(inputdatapath+msfile)

Imaging

# set msfile to include the full path
msfile = inputdatapath+msfile

imgname='Venus-mosaic-mfs'

delete_tcleanimages(imgname)

Run mosaic with usepointing = True to correctly shift the beam

ret=tclean(vis=msfile, field='Venus', imagename=imgname, imsize=[480, 432], cell=['0.14arcsec'],
           phasecenter='VENUS', specmode='mfs', gridder='mosaic', usepointing=True, niter=100)
ret['iterdone']

100

Examination of the resultant image

Convert a CASA image to FITS image for displaying the image via aplpy

# Define functions to use in displaying the image
def displayImage(imgname, markers={}, title=''):
    ''' convert casa image to fits and display the fits image'''
    # export to a FITS image (overwrite if the output exists)
    from casatools import quanta
    _qa = quanta()
    fitsimagename = imgname+'fits'
    exportfits(imagename=imgname, fitsimage=fitsimagename, overwrite=True)
    fits.setval(fitsimagename, 'TIMESYS', value='utc')
    fig = pl.figure(figsize=(15,15))
    img = aplpy.FITSFigure(fitsimagename, subplot=[0.1,0.1,0.5,0.5])
    img.show_colorscale()
    if markers!={}:
       for ky, pos in markers.items():
          if 'ra' in pos and 'dec' in pos:
             print(f'ky={ky}')
             if 'ext' in ky:
                 color = 'red'
             else:
                 color = 'yellow'
             img.show_markers(pos['ra'], pos['dec'], edgecolor=color, marker='o',s=10, alpha=0.5)
             marker_ra = _qa.time(_qa.quantity(pos['ra'],'deg'),prec=9)[0]
             marker_dec = _qa.angle(_qa.quantity(pos['dec'],'deg'),prec=9)[0]
             print(f'{ky} at: {marker_ra}, {marker_dec}')
    pl.title(title)
    fig.canvas.draw()
    print('image center:',printImageCenter(imgname))


def printImageCenter(imgname):
    from casatools import image, quanta
    _ia = image()
    _qa = quanta()
    _ia.open(imgname)
    shape = _ia.shape()
    csys = _ia.coordsys()
    center = csys.toworld([shape[0]/2.,shape[1]/2.])
    _ia.done()
    return _qa.time(_qa.quantity(center['numeric'][0], 'rad'),prec=9)[0], _qa.angle(_qa.quantity(center['numeric'][1],'rad'),prec=9)[0]

# Returns direction of the ephemeris object at a given time (parallax collected)
def ephem_dir(ephemtab, refep, observatory):
    from casatools import measures
    _me = measures()
    _me.framecomet(ephemtab)
    _me.doframe(_me.observatory(observatory))
    _me.doframe(_me.epoch('utc',refep))
    return _me.measure(_me.measure(_me.direction('COMET'), 'AZELGEO'),'ICRS')

Get the direction of Venus at the first timestamp of the select data from the ephemeris table

# Find the path of the attached ephem. table for fieldid=0
ephemtab = get_attachedEphemtablepath(msfile,'Venus')
print(ephemtab)

mytestdata/measurementset/alma/venus_ephem_test.ms/FIELD/EPHEM0_Venus_58491.4.tab/

from casatasks import imhead
iminfo=imhead(imgname+'.image', mode='list')
iminfo['date-obs']

'2019/01/08/11:21:50.208000'

ineph_dir = ephem_dir(ephemtab, iminfo['date-obs'],'ALMA')
print(ineph_dir)
ineph_rad = qa.convert(ineph_dir['m0'],'deg')
ineph_decd = qa.convert(ineph_dir['m1'],'deg')
inephmarker = {'ineph':{'ra':ineph_rad['value'], 'dec':ineph_decd['value']}}

{'m0': {'unit': 'rad', 'value': -2.1053946095161686}, 'm1': {'unit': 'rad', 'value': -0.29607600119188454}, 'refer': 'ICRS', 'type': 'direction'}

displayImage(imgname+'.image', markers=inephmarker, title='Venus mosaic image using the default ephem table in the data repository')

INFO: Setting slices=[0, 0] [aplpy.core]
INFO: Auto-setting vmin to -3.726e+00 [aplpy.core]
INFO: Auto-setting vmax to  3.681e+00 [aplpy.core]
ky=ineph
ineph at: 15:57:28.746, -016.57.50.059
image center: ('15:57:28.543', '-016.57.49.411')

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 58491.473498 from DATE-OBS'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to   -67.754929 from OBSGEO-[XYZ].
Set OBSGEO-B to   -23.022886 from OBSGEO-[XYZ].
Set OBSGEO-H to     5053.796 from OBSGEO-[XYZ]'. [astropy.wcs.wcs]

<Figure size 1500x1500 with 0 Axes>

Note: the marker in yellow is set to the direction based on the default ephemeris table at the first time of the selected data and its coordintes are shifted with the coordinates of the center of the image. This is because in the observation, the updated ephemeris data was used.

print("internal ephem dir = ",qa.time(ineph_dir['m0'],prec=9),qa.angle(ineph_dir['m1'],prec=9))
printImageCenter(imgname+'.image')

internal ephem dir =  ['15:57:28.746'] ['-016.57.50.059']

('15:57:28.543', '-016.57.49.411')

stats = imstat(imgname+'.image')
print(f"{stats['max']}@{stats['maxpos']}")

[3.33142209]@[155 170   0   0]