Repository: http://github.com/matteobachetti/srt-single-dish-tools
Documentation: http://srt-single-dish-tools.readthedocs.io/en/latest/
You need to install the following libraries as described below:
astropy
pandas and seaborn for the example with a pair plot
The SRT single dish tools, of course, with their dependencies.
Install Anaconda then, in a shell:
$ conda create -n py3-clean python=3.6 astropy seaborn numpy scipy matplotlib h5py statsmodels pyyaml numba pandas
$ source activate py3-clean
$ conda install -c astropy pyregion
$ git clone https://github.com/matteobachetti/srt-single-dish-tools
$ cd srt-single-dish-tools
$ pip install .Open a shell, then
pip install git+https://github.com/matteobachetti/srt-single-dish-tools.git astropy seaborn numpy scipy matplotlib h5py statsmodels pyyaml numba pandas pyregionAnd wait some hours for the compilation to complete :)
In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
from __future__ import print_function, division
from astropy.table import Table
import numpy as np
from matplotlib import pyplot as plt
from seaborn import pairplot, lmplot
import warnings
from pandas import read_csv
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
from astropy.visualization import ImageNormalize, LogStretch, MinMaxInterval
warnings.filterwarnings('ignore')
In [2]:
!SDTfake -h
In [3]:
!SDTfake -s 20 -b 5 --spacing 0.3 --geometry 60 60 60 60 --beam-width 2.5
In [4]:
!ls sim/
Some of these observations might have been done in different bands, or using different receivers, and you might have lost the list of observations (or the user was not the observer). The script SDTinspector is there to help, dividing the observations in groups based on observing time, backend, receiver, etc.
Option -d specifies to create configuration files for future analysis, that contain all pointings of a given source with a given configuration and all calibrator observations in the +- 6 hours around the observations.
In [5]:
!SDTinspect sim/*/ -d
In [6]:
!cat CCB_TP_Dummy_Obs0.ini
We can modify the configuration file, for example changing the pixel size to 1.5, or run SDTinspect with the --options option:
In [9]:
!SDTinspect sim/*/ -d --options '{"pixel_size": 0.7}'
In [10]:
!cat CCB_TP_Dummy_Obs0.ini
The configuration file created by SDTinspect can be analyzed by all other scripts. We are now ready to use it to calculate the calibration parameters.
The calibration is calculated by fitting a Gaussian to cross scans over a "calibrator", a source chosen for its high flux stability over time. The fit can be good or bad, so it is always recommended to use the --show and --check options, that plot some diagnostic pictures like that above and save it in the same directory as the calibrator scans, and check that the fluxes of calibrators calculated by using the calibration on their raw data is consistent with the real flux of the calibrators themselves.
The warnings appearing below are mostly innocuous. They are generated by the fact that the main peak is very strong and the gradient in the scan is interpreted as an outlier.
In [12]:
!SDTcal -c CCB_TP_Dummy_Obs0.ini --show --check
Here is an example of a calibrator scan, quite noisy, fit with a Gaussian plus a slopy baseline (File sim/calibrator1/0_Dec_scanfit/Feed0_chan0.png)
After this step, we have a calibration table saved in readable format (.csv) and in binary format (.hdf5). Let's see what the csv table shows us. For example, we can plot the gain curves.
In [13]:
table = read_csv('CCB_TP_Dummy_Obs0_cal.csv')
table['Flux/T'] = table['Flux'] / table['Source_temperature']
In [14]:
plt.figure(figsize=(15, 8))
good = (table['Source'] == 'DummyCal')&(table['Chan'] == 'Feed0_LCP')
plt.scatter(table['Elevation'][good], table['Flux/T'][good],
label='DummyCal', marker='o')
good = (table['Source'] == 'DummyCal2')&(table['Chan'] == 'Feed0_LCP')
plt.scatter(table['Elevation'][good], table['Flux/T'][good],
label='DummyCal2', marker='o')
plt.title('Channel: LCP')
plt.xlabel('Elevation')
plt.ylabel('Temperature / Flux Density (K / Jy)')
plt.legend()
plt.figure(figsize=(15, 8))
good = (table['Source'] == 'DummyCal')&(table['Chan'] == 'Feed0_RCP')
plt.scatter(table['Elevation'][good], table['Flux/T'][good],
label='DummyCal', marker='s')
good = (table['Source'] == 'DummyCal2')&(table['Chan'] == 'Feed0_RCP')
plt.scatter(table['Elevation'][good], table['Flux/T'][good],
label='DummyCal2', marker='s')
plt.title('Channel: RCP')
plt.xlabel('Elevation')
plt.ylabel('Temperature / Flux Density (K / Jy)')
plt.legend()
Out[14]:
In [15]:
fig = plt.figure(figsize=(15., 8))
g = lmplot(x="Elevation", y="Flux/T", row='Source', col='Chan',
truncate=True, size=5, data=table)
In [16]:
fig = plt.figure(figsize=(15., 8))
g = lmplot(x="Elevation", y="Flux/T",
truncate=True, data=table)
But we can do everything we want with that information... For example a pairplot
In [17]:
pairplot(table, vars=['Elevation', 'Azimuth', 'Flux/Counts', 'Flux/T', 'Flux Integral/Counts', 'Calculated Flux', 'Flux'])
Out[17]:
In [18]:
table
Out[18]:
This is done with SDTimage. WCS coordinates are automatically calculated and results are written into a ds9- and casaviewer-compatible FITS file.
To visualize, I'm using the astropy library here.
The -e <regionfile> tells the baseline subtraction algorithm to ignore a given part of the map where strong sources are present. This region file has to be in DS9 format and contain coordinates in the fk5 reference frame. E.g.
In [19]:
!cat center.reg
As usual, please ignore below the warnings about outliers. This meanse that the baseline subtraction algorithm is ignoring some points from the fit.
In [20]:
!SDTimage -c CCB_TP_Dummy_Obs0.ini --sub --refilt -e center.reg
In [21]:
def plot_image(file, ext, logstretch=False):
hdu = fits.open(file)[ext]
wcs = WCS(hdu.header)
plt.figure(figsize=(10, 10))
ax = plt.subplot(projection=wcs)
img = hdu.data
if logstretch:
# Create an ImageNormalize object
norm = ImageNormalize(img, interval=MinMaxInterval(),
stretch=LogStretch())
im = ax.imshow(img, origin='lower', cmap='viridis', norm=norm)
else:
im = ax.imshow(img, origin='lower', cmap='viridis')
print("Image max is", np.max(img))
print("Image sum is", np.sum(img))
print("Image median is", np.median(img))
ax.coords.grid(True, color='white', ls='solid')
ax.coords[0].set_axislabel('Right Ascension (J2000)')
ax.coords[1].set_axislabel('Declination (J2000)')
overlay = ax.get_coords_overlay('galactic')
overlay.grid(color='white', ls='dotted')
overlay[0].set_axislabel('Galactic Longitude')
overlay[1].set_axislabel('Galactic Latitude')
plt.colorbar(im)
plot_image('CCB_TP_Dummy_Obs0.fits', 1)
In [22]:
plot_image('CCB_TP_Dummy_Obs0.fits', 1, logstretch=True)
The log-stretched image shows that the noisy scans can create stripes in the map. There are a couple of different ways to deal with this in SDT. One is interactive and the other is automatic. Let us now use the...
We can cleanup the scans interactively. A window opens and we can point and click to select all scans passing through a pixel, and flag bad scans, zap intervals, redo baseline fitting, etc.
Please follow the instructions that appear in the shell. For example, to select all scans passing through a given pixel and +- 1 pixel around, point the pixel with the mouse and press "a". A new window will open showing all the scans. Click on a scan (it becomes green), and filter out intervals in the scan, fit a new baseline, flag scans, etc. by following the instructions appearing in the shell.
In [19]:
!SDTimage -c CCB_TP_Dummy_Obs0.ini --interactive
Now we can calibrate the image, using the calibration calculated above, using the option --calibrate <CALFILE>.hdf5 and specifying the unit (Jy/beam, Jy/pixel, or Jy/sterad) with -u <UNIT>
Note: instead of the calibration file, we can give the HDF5 file dumped by SDTimage in the previous processing.
Instead of SDTimage -c CCB_TP_Dummy_Obs0.ini, we could use SDTimage CCB_TP_Dummy_Obs0_dump.hdf5
In [23]:
!SDTimage -c CCB_TP_Dummy_Obs0.ini -u 'Jy/pixel' --calibrate CCB_TP_Dummy_Obs0_cal.hdf5
In [24]:
plot_image('CCB_TP_Dummy_Obs0_cal.fits', 1, logstretch=True)
In [36]:
!SDTimage CCB_TP_Dummy_Obs0_dump.hdf5 -u 'Jy/beam' --calibrate CCB_TP_Dummy_Obs0_cal.hdf5 --destripe
In [37]:
plot_image('CCB_TP_Dummy_Obs0_cal_destripe.fits', 1, logstretch=True)
In [ ]: