Load and Map an APEX observation
You must download the file E-085.B-0964A-2010_merge.apex from the archive first
Requires:
An "install script":
Mac:
curl -O http://09c8d0b2229f813c1b93-c95ac804525aac4b6dba79b00b39d1d3.r79.cf1.rackcdn.com/Anaconda-2.0.1-MacOSX-x86_64.shLinux:
wget http://09c8d0b2229f813c1b93-c95ac804525aac4b6dba79b00b39d1d3.r79.cf1.rackcdn.com/Anaconda-2.0.1-Linux-x86_64.shBoth:
conda install scipy
conda install astropy
pip install https://github.com/keflavich/sdpy/archive/master.zip
pip install https://github.com/keflavich/FITS_tools/archive/master.zip
pip install https://bitbucket.org/pyspeckit/pyspeckit/get/master.tar.gz
pip install https://github.com/keflavich/mpl_plot_templates/archive/master.zip# You can acquire the data with astroquery once a few kinks are worked out (don't try this yet):
from astroquery.eso import Eso
Eso.login('aginsburg') # change this to your username; you may need to enter your password in the terminal!
# this dataset is public, but you still need to log in
filenames = Eso.retrieve_data('E-085.B-0964A')
import os
localpath = os.path.split(filenames[-1])[-1]
if not os.path.exists(localpath):
os.symlink(filenames[-1],localpath)
In [1]:
from shfi_otf_pipeline import make_apex_cubes
import time # for performance analysis
In [2]:
# for bigger figures
import pylab as pl
pl.rcParams['figure.figsize'] = (10.0, 8.0)
In [3]:
filename = 'E-085.B-0964A.2010JUL08/E-085.B-0964A-2010_merge.apex'
datapath = 'E-085.B-0964A.2010JUL08/'
dataset = 'E-085.B-0964A-2010'
In []:
# there is a specific wrapper for this dataset
# the lines are split into low/high frequencies
t0 = time.time()
make_apex_cubes.build_cube_ao(window='high',
datasets=[dataset],
datapath=make_apex_cubes.aorawpath,
outpath='./',
pcakwargs={'smoothing_scale':50},
freq=True,
)
print "Total time: ",time.time()-t0," seconds"
In [5]:
import glob
[x for x in glob.glob(os.path.join(make_apex_cubes.diagplotdir, dataset)+"_AP-H201-F101*diagnostics.png") if '_obs' not in x]
# List the diagnostic plots, excluding the (many!) observation-by-observation diagnostics
Out[5]:
In [6]:
# Same for TSYS
tsys_plots = [x for x in glob.glob(os.path.join(make_apex_cubes.diagplotdir, dataset)+"_AP-H201-F101*tsys.png") if '_obs' not in x]
tsys_plots
Out[6]:
In [7]:
from IPython.display import Image
In [8]:
# The overall diagnostic plot before TSYS flagging:
Image('/Users/adam/work/h2co/apex/diagnostic_plots/E-085.B-0964A-2010_AP-H201-F101_preflag_diagnostics.png',)
Out[8]:
In [9]:
# Before time-domain subtraction, after TSYS subtraction
Image(os.path.join(make_apex_cubes.diagplotdir, 'E-085.B-0964A-2010_AP-H201-F101_presub_diagnostics.png'))
Out[9]:
In [10]:
# After time-domain subtraction: not much difference (which is fine)
Image(os.path.join(make_apex_cubes.diagplotdir, 'E-085.B-0964A-2010_AP-H201-F101_diagnostics.png'))
Out[10]:
In [11]:
# Before flagging out bad TSYS
# The top plot shows TSYS as a function of time (with x & y axes flipped)
print tsys_plots[0]
Image(tsys_plots[0])
Out[11]:
In [12]:
# After flagging out bad TSYS
# There are some extremely low-noise points; these are probably spurious and should be flagged out
print tsys_plots[1]
Image(tsys_plots[1])
Out[12]:
In [13]:
# Regions with lines masked out for baseline fitting
# The mask used is not being stored right now though...
mask_plots = [x for x in glob.glob(os.path.join(make_apex_cubes.diagplotdir, dataset)+"_AP-H201-F101*masked.png") if '_obs' not in x]
mask_plots
Out[13]:
In [14]:
print mask_plots[0]
Image(mask_plots[0])
Out[14]:
In [15]:
# Display some of the results
import aplpy
F = aplpy.FITSFigure('APEX_H2CO_Ao_Freq_high.fits', slices=[3324],)
F.show_grayscale()
F.show_contour('APEX_H2CO_Ao_Freq_high.fits', slices=[3290], levels=[-0.5,0.5,0.75,1,1.25,1.5], colors=['b'])
F.show_contour('APEX_H2CO_Ao_Freq_high.fits', slices=[3245], levels=[-0.5,0.5,0.75,1,1.25,1.5], colors=['r'])
In [16]:
# Do it a second time, but with PCA_cleaning disabled
t0 = time.time()
if not os.path.exists('nopca'): os.mkdir('nopca/')
make_apex_cubes.build_cube_ao(window='high',
datasets=[dataset],
datapath=make_apex_cubes.aorawpath,
outpath='./nopca/',
pca_clean=False,
freq=True,
)
print "Total time without PCA cleaning: ",time.time()-t0," seconds"
In [17]:
F = aplpy.FITSFigure('nopca/APEX_H2CO_Ao_Freq_high.fits', slices=[3324],)
F.show_grayscale()
F.show_contour('nopca/APEX_H2CO_Ao_Freq_high.fits', slices=[3290], levels=[-0.5,0.5,0.75,1,1.25,1.5], colors=['b'])
F.show_contour('nopca/APEX_H2CO_Ao_Freq_high.fits', slices=[3245], levels=[-0.5,0.5,0.75,1,1.25,1.5], colors=['r'])
In [18]:
# Time PCA: This is a technique to mitigate the possibility of subtracting spectral information
t0 = time.time()
if not os.path.exists('timepca'): os.mkdir('timepca/')
make_apex_cubes.build_cube_ao(window='high',
datasets=[dataset],
datapath=make_apex_cubes.aorawpath,
outpath='./timepca/',
timewise_pca=True,
pca_clean=True,
freq=True,
)
print "Total time: ",time.time()-t0," seconds"
In [19]:
F = aplpy.FITSFigure('timepca/APEX_H2CO_Ao_Freq_high.fits', slices=[3324],)
F.show_grayscale()
F.show_contour('timepca/APEX_H2CO_Ao_Freq_high.fits', slices=[3290], levels=[-0.5,0.5,0.75,1,1.25,1.5], colors=['b'])
F.show_contour('timepca/APEX_H2CO_Ao_Freq_high.fits', slices=[3245], levels=[-0.5,0.5,0.75,1,1.25,1.5], colors=['r'])
In [20]:
# PCA -> noPCA comparison
from astropy.io import fits
nopca = fits.open('nopca/APEX_H2CO_Ao_Freq_high_downsampled.fits')
pca = fits.open('APEX_H2CO_Ao_Freq_high_downsampled.fits')
timepca = fits.open('timepca/APEX_H2CO_Ao_Freq_high_downsampled.fits')
pca_diff = fits.PrimaryHDU(pca[0].data-nopca[0].data, header=pca[0].header)
pca_diff.writeto('APEX_H2CO_Ao_Freq_high_downsampled_PCAminusNoPCA.fits', clobber=True)
time_pca_diff = fits.PrimaryHDU(timepca[0].data-nopca[0].data, header=pca[0].header)
time_pca_diff.writeto('APEX_H2CO_Ao_Freq_high_downsampled_TimePCAminusNoPCA.fits', clobber=True)
In [21]:
# The negative feature here shows that the naive PCA-eigenspectrum approach has failed!
F = aplpy.FITSFigure('APEX_H2CO_Ao_Freq_high_downsampled_PCAminusNoPCA.fits', slices=[819],)
F.show_grayscale()
F.show_contour('APEX_H2CO_Ao_Freq_high_downsampled_PCAminusNoPCA.fits', slices=[213])
In [22]:
# By contrast, timewise PCA shows features that look more like artifacts
F = aplpy.FITSFigure('APEX_H2CO_Ao_Freq_high_downsampled_TimePCAminusNoPCA.fits', slices=[819],)
F.show_grayscale()
F.show_contour('APEX_H2CO_Ao_Freq_high_downsampled_TimePCAminusNoPCA.fits', slices=[213], levels=3)
In [23]:
# Cleanup code: this removes excess junk text from the notebook to decrease its size
import json
def strip_stream(filename='ExampleAPEXMappingReduction.ipynb'):
jay = json.load(open(filename))
for cell in jay['worksheets'][0]['cells']:
if 'outputs' in cell:
cell['outputs'] = [x for x in cell['outputs'] if x['output_type'] != 'stream']
with open(filename,'w') as outf:
json.dump(jay, outf)
return jay
jay = strip_stream()
In [23]: