Querying ALMA archive for OrionKL pointings and plotting them on a 2MASS image

In [91]:
from astroquery.alma import Alma
from astroquery.skyview import SkyView
import string
from astropy import units as u
from astropy import coordinates
import pylab as pl
import aplpy

Retrieve OrionKL 2MASS K-band image:

In [92]:
orionkl_coords = coordinates.SkyCoord.from_name('Orion KL')
In [93]:
orionkl_images = SkyView.get_images(position='Orion KL', survey=['2MASS-K'], pixels=500)
orionkl_images
FITS files must be read as binaries; error is likely.
Out[93]:
[[<astropy.io.fits.hdu.image.PrimaryHDU at 0x11f811a50>]]

Retrieve ALMA archive information including private data and non-science fields:

In [94]:
orionkl = Alma.query_region(coordinate=orionkl_coords, radius=4*u.arcmin, public=False, science=False)
In [95]:
orionkl
Out[95]:
<Table masked=True length=891>
Project codeSource nameRADecBandFrequency resolutionIntegrationRelease dateFrequency supportVelocity resolutionPol productsObservation datePI namePWVMember ous idAsdm uidProject titleProject typeScan intent
degdegkHzsGHzm / smm
string128string2048float64float64int32float64float64objectstring21600float64string512objectstring2048float32string256string256string2048string128string2048
2011.0.00199.SOrion H2O maser outburst83.808854166666663-5.3768016666666663631250.010.6772013-06-14[240.01..241.99GHz,31250.00kHz, XX YY] U [242.01..243.99GHz,31250.00kHz, XX YY] U [256.01..257.99GHz,31250.00kHz, XX YY] U [258.01..259.99GHz,31250.00kHz, XX YY]37474.05724998501XX YY2012-04-08 20:19:27Hirota, Tomoya2.9074531uid://A002/X303d22/X4cuid://A002/X3ca066/X114Bursting Water Maser Feature in Orion KLSTARGET
2011.0.00199.SOrion H2O maser outburst83.808854166666663-5.37680166666666637244.14062527.9310000000000012013-12-03[320.98..321.45GHz,244.14kHz, XX YY] U [322.12..322.59GHz,244.14kHz, XX YY] U [334.43..334.90GHz,244.14kHz, XX YY] U [335.98..336.45GHz,244.14kHz, XX YY]222.79811202502449XX YY2012-10-21 10:14:20Hirota, Tomoya0.49233711uid://A002/X303d22/X5auid://A002/X518d2b/Xb3aBursting Water Maser Feature in Orion KLSTARGET
2011.0.00199.SOrion H2O maser outburst83.808854166666663-5.37680166666666637244.14062532.1079999999999972013-12-03[320.98..321.45GHz,244.14kHz, XX YY] U [322.12..322.59GHz,244.14kHz, XX YY] U [334.43..334.90GHz,244.14kHz, XX YY] U [335.98..336.45GHz,244.14kHz, XX YY]222.79811202502449XX YY2012-07-16 11:07:23Hirota, Tomoya0.66106993uid://A002/X303d22/X5auid://A002/X460fc9/Xa91Bursting Water Maser Feature in Orion KLSTARGET
2011.0.00199.SOrion H2O maser outburst83.808854166666663-5.37680166666666637244.14062532.7169999999999992013-12-03[320.98..321.45GHz,244.14kHz, XX YY] U [322.12..322.59GHz,244.14kHz, XX YY] U [334.43..334.90GHz,244.14kHz, XX YY] U [335.98..336.45GHz,244.14kHz, XX YY]222.79811202502449XX YY2012-08-25 10:30:47Hirota, Tomoya0.73374134uid://A002/X303d22/X5auid://A002/X4afc39/Xa3aBursting Water Maser Feature in Orion KLSTARGET
2011.0.00028.SOrionField1-183.817916666666648-5.38955555555555597976.562532.6540000000000032014-01-03[342.36..344.24GHz,976.56kHz, XX YY] U [344.24..346.11GHz,976.56kHz, XX YY] U [354.28..356.16GHz,976.56kHz, XX YY] U [355.79..357.67GHz,976.56kHz, XX YY]835.99889987775487XX YY2012-10-24 05:23:07Mann, Rita0.72024173uid://A002/X327408/X205uid://A002/X51ee7e/X1c2The Effect of Extreme Environment on Protoplanetary Disks in OrionSTARGET
2011.0.00028.SOrionField1-183.817916666666648-5.38955555555555597976.562532.6640000000000012014-01-03[342.36..344.24GHz,976.56kHz, XX YY] U [344.24..346.11GHz,976.56kHz, XX YY] U [354.28..356.16GHz,976.56kHz, XX YY] U [355.79..357.67GHz,976.56kHz, XX YY]835.99889987775487XX YY2012-10-24 08:48:42Mann, Rita0.70387971uid://A002/X327408/X205uid://A002/X51ee7e/X4ddThe Effect of Extreme Environment on Protoplanetary Disks in OrionSTARGET
2011.0.00028.SOrionField1-183.817916666666648-5.38955555555555597976.562532.6940000000000032014-01-03[342.36..344.24GHz,976.56kHz, XX YY] U [344.24..346.11GHz,976.56kHz, XX YY] U [354.28..356.16GHz,976.56kHz, XX YY] U [355.79..357.67GHz,976.56kHz, XX YY]835.99889987775487XX YY2012-10-24 07:14:43Mann, Rita0.73464209uid://A002/X327408/X205uid://A002/X51ee7e/X35fThe Effect of Extreme Environment on Protoplanetary Disks in OrionSTARGET
2011.0.00028.SOrionField1-183.817916666666648-5.38955555555555597976.562532.7070000000000012014-01-03[342.36..344.24GHz,976.56kHz, XX YY] U [344.24..346.11GHz,976.56kHz, XX YY] U [354.28..356.16GHz,976.56kHz, XX YY] U [355.79..357.67GHz,976.56kHz, XX YY]835.99889987775487XX YY2012-10-24 03:47:01Mann, Rita0.70661682uid://A002/X327408/X205uid://A002/X51ee7e/X44The Effect of Extreme Environment on Protoplanetary Disks in OrionSTARGET
2011.0.00028.SOrionField1-183.817916666666648-5.38955555555555597976.562532.7449999999999972014-01-03[342.36..344.24GHz,976.56kHz, XX YY] U [344.24..346.11GHz,976.56kHz, XX YY] U [354.28..356.16GHz,976.56kHz, XX YY] U [355.79..357.67GHz,976.56kHz, XX YY]835.99889987775487XX YY2012-10-24 10:30:52Mann, Rita0.71500158uid://A002/X327408/X205uid://A002/X51ee7e/X67aThe Effect of Extreme Environment on Protoplanetary Disks in OrionSTARGET
.........................................................
2012.1.00352.SOrion_Bar83.835833333333326-5.42222222222222167976.56252226.3829999999998[345.40..346.39GHz,488.28kHz, XX YY] U [356.22..357.22GHz,488.28kHz, XX YY] U [357.29..358.29GHz,488.28kHz, XX YY]410.93402462306449XX YY2014-10-07 09:28:44Goicoechea, Javier1.3864298uid://A001/X144/X21uid://A002/X8f1c7c/X1528The fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET WVR
2012.1.00352.SOrion_Bar83.835833333333326-5.42222222222222167976.56252226.652[345.40..346.40GHz,488.28kHz, XX YY] U [356.22..357.22GHz,488.28kHz, XX YY] U [357.29..358.29GHz,488.28kHz, XX YY]410.93362488666679XX YY2014-10-05 09:41:58Goicoechea, Javier0.28458548uid://A001/X144/X21uid://A002/X8ec83d/XbaaThe fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET WVR
2012.1.00352.SOrion_Bar83.835833333333326-5.42222222222222167976.56252227.9760000000001[345.40..346.40GHz,488.28kHz, XX YY] U [356.22..357.22GHz,488.28kHz, XX YY] U [357.29..358.29GHz,488.28kHz, XX YY]410.93170820116194XX YY2014-09-10 10:11:13Goicoechea, Javier0.44827378uid://A001/X144/X21uid://A002/X8b7acb/X355The fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET WVR
2012.1.00352.SOrion_Bar83.835833333333326-5.42222222222222167976.56252228.877[345.40..346.40GHz,488.28kHz, XX YY] U [356.22..357.22GHz,488.28kHz, XX YY] U [357.29..358.29GHz,488.28kHz, XX YY]410.93386189627125XX YY2014-10-07 08:27:54Goicoechea, Javier1.3269638uid://A001/X144/X21uid://A002/X8f1c7c/X1269The fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET WVR
2012.1.00352.SOrion_Bar83.835833333333326-5.42222222222222167976.56252230.4259999999999[345.40..346.39GHz,488.28kHz, XX YY] U [356.22..357.22GHz,488.28kHz, XX YY] U [357.29..358.29GHz,488.28kHz, XX YY]410.93415100721489XX YY2014-10-08 08:55:11Goicoechea, Javier0.62932658uid://A001/X144/X21uid://A002/X8f3a3d/X16f0The fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET WVR
2012.1.00352.SOrion_Bar83.83667697797128-5.41576521119719437976.562537.659999999999997[345.34..346.34GHz,488.28kHz, XX YY] U [356.21..357.21GHz,488.28kHz, XX YY] U [357.24..358.24GHz,488.28kHz, XX YY]410.9394076007394XX YY2014-04-08 23:24:29Goicoechea, Javier--uid://A002/X6dddc4/X6cuid://A002/X7ea111/X22dThe fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET
2012.1.00352.SOrion_Bar83.837432914740404-5.42624693968655517976.562537.659999999999997[345.34..346.34GHz,488.28kHz, XX YY] U [356.21..357.21GHz,488.28kHz, XX YY] U [357.24..358.24GHz,488.28kHz, XX YY]410.9394076007394XX YY2014-04-08 23:24:29Goicoechea, Javier--uid://A002/X6dddc4/X6cuid://A002/X7ea111/X22dThe fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET
2012.1.00352.SOrion_Bar83.838465543570635-5.42241036078388837976.562537.659999999999997[345.34..346.34GHz,488.28kHz, XX YY] U [356.21..357.21GHz,488.28kHz, XX YY] U [357.24..358.24GHz,488.28kHz, XX YY]410.9394076007394XX YY2014-04-08 23:24:29Goicoechea, Javier--uid://A002/X6dddc4/X6cuid://A002/X7ea111/X22dThe fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET
2012.1.00352.SOrion_Bar83.83949817240088-5.41857378188122327976.562537.659999999999997[345.34..346.34GHz,488.28kHz, XX YY] U [356.21..357.21GHz,488.28kHz, XX YY] U [357.24..358.24GHz,488.28kHz, XX YY]410.9394076007394XX YY2014-04-08 23:24:29Goicoechea, Javier--uid://A002/X6dddc4/X6cuid://A002/X7ea111/X22dThe fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET
2012.1.00352.SOrion_Bar83.842319366830452-5.42138235256522017976.562537.659999999999997[345.34..346.34GHz,488.28kHz, XX YY] U [356.21..357.21GHz,488.28kHz, XX YY] U [357.24..358.24GHz,488.28kHz, XX YY]410.9394076007394XX YY2014-04-08 23:24:29Goicoechea, Javier--uid://A002/X6dddc4/X6cuid://A002/X7ea111/X22dThe fundamental structure of molecular cloud edges: from clumps to photoevaporationSTARGET

Parse components of the ALMA data. Specifically, find the frequency support - the frequency range covered - and convert that into a central frequency for beam radius estimation.

In [96]:
def parse_frequency_support(frequency_support_str):
    supports = frequency_support_str.split("U")
    freq_ranges = [(float(sup.strip('[] ').split("..")[0]),
                    float(sup.strip('[] ').split("..")[1].split(',')[0].strip(string.letters)))
                   *u.Unit(sup.strip('[] ').split("..")[1].split(',')[0].strip(string.punctuation+string.digits))
                   for sup in supports]
    return u.Quantity(freq_ranges)

def approximate_primary_beam_sizes(frequency_support_str):
    freq_ranges = parse_frequency_support(frequency_support_str)
    beam_sizes = [(1.22*fr.mean().to(u.m, u.spectral())/(12*u.m)).to(u.arcsec,
                                                                     u.dimensionless_angles())
                  for fr in freq_ranges]
    return u.Quantity(beam_sizes)
In [97]:
primary_beam_radii = [approximate_primary_beam_sizes(row['Frequency support']) for row in orionkl]

Compute primary beam parameters for the public and private components of the data for plotting below.

In [109]:
bands = np.unique(orionkl['Band'])
print "The bands used include: "
print bands
band_colors_priv = dict(zip(bands, ('red','darkred','orange','brown','maroon')))
band_colors_pub = dict(zip(bands, ('blue','cyan','green','turquoise','teal')))
The bands used include: 
Band
----
   3
   4
   6
   7
   9
In [99]:
private_circle_parameters = {band: [(row['RA'],row['Dec'],np.mean(rad).to(u.deg).value)
                             for row,rad in zip(orionkl, primary_beam_radii)
                             if row['Release date']!='' and row['Band']==band]
                             for band in bands}
public_circle_parameters = {band: [(row['RA'],row['Dec'],np.mean(rad).to(u.deg).value)
                             for row,rad in zip(orionkl, primary_beam_radii)
                             if row['Release date']=='' and row['Band']==band]
                             for band in bands}

unique_private_circle_parameters = {band: np.array(list(set(private_circle_parameters[band])))
                                    for band in bands}
unique_public_circle_parameters = {band: np.array(list(set(public_circle_parameters[band])))
                                   for band in bands}

for band in bands:
    print "BAND {0}".format(band)
    privrows = sum((orionkl['Band']==band) & (orionkl['Release date'] != ''))
    pubrows  = sum((orionkl['Band']==band) & (orionkl['Release date'] == ''))
    print "PUBLIC:  Number of rows: {0}.  Unique pointings: {1}".format(pubrows, len(unique_public_circle_parameters[band]))
    print "PRIVATE: Number of rows: {0}.  Unique pointings: {1}".format(privrows, len(unique_private_circle_parameters[band]))
BAND 3
PUBLIC:  Number of rows: 2.  Unique pointings: 1
PRIVATE: Number of rows: 1.  Unique pointings: 1
BAND 4
PUBLIC:  Number of rows: 0.  Unique pointings: 0
PRIVATE: Number of rows: 2.  Unique pointings: 2
BAND 6
PUBLIC:  Number of rows: 388.  Unique pointings: 388
PRIVATE: Number of rows: 262.  Unique pointings: 155
BAND 7
PUBLIC:  Number of rows: 25.  Unique pointings: 14
PRIVATE: Number of rows: 190.  Unique pointings: 151
BAND 9
PUBLIC:  Number of rows: 0.  Unique pointings: 0
PRIVATE: Number of rows: 21.  Unique pointings: 17

Show all of the private observation pointings that have been acquired

In [110]:
fig = aplpy.FITSFigure(orionkl_images[0])
fig.show_grayscale(stretch='arcsinh')
for band in bands:
    if unique_private_circle_parameters[band].any():
        fig.show_circles(unique_private_circle_parameters[band][:,0],
                         unique_private_circle_parameters[band][:,1],
                         unique_private_circle_parameters[band][:,2],
                         color=band_colors_priv[band], alpha=0.2)
INFO:astropy:Auto-setting vmin to  2.393e+02
INFO:astropy:Auto-setting vmax to  3.155e+03
INFO: Auto-setting vmin to  2.393e+02 [aplpy.core]
INFO: Auto-setting vmax to  3.155e+03 [aplpy.core]

In principle, all of the pointings shown below should be downloadable from the archive:

In [111]:
fig = aplpy.FITSFigure(orionkl_images[0])
fig.show_grayscale(stretch='arcsinh')
for band in bands:
    if unique_public_circle_parameters[band].any():
        fig.show_circles(unique_public_circle_parameters[band][:,0],
                         unique_public_circle_parameters[band][:,1],
                         unique_public_circle_parameters[band][:,2],
                         color=band_colors_pub[band], alpha=0.2)
INFO:astropy:Auto-setting vmin to  2.393e+02
INFO:astropy:Auto-setting vmax to  3.155e+03
INFO: Auto-setting vmin to  2.393e+02 [aplpy.core]
INFO: Auto-setting vmax to  3.155e+03 [aplpy.core]

Use pyregion to write the observed regions to disk. Pyregion has a very awkward API; there is (in principle) work in progress to improve that situation but for now one must do all this extra work.

In [112]:
import pyregion
from pyregion.parser_helper import Shape
prv_regions = {band: pyregion.ShapeList([Shape('circle',[x,y,r]) for x,y,r in private_circle_parameters[band]])
               for band in bands}
pub_regions = {band: pyregion.ShapeList([Shape('circle',[x,y,r]) for x,y,r in public_circle_parameters[band]])
               for band in bands}
for band in bands:
    circle_pars = np.vstack([x for x in (private_circle_parameters[band],
                                    public_circle_parameters[band]) if any(x)])
    for r,(x,y,c) in zip(prv_regions[band]+pub_regions[band],
                         circle_pars):
        r.coord_format = 'fk5'
        r.coord_list = [x,y,c]
        r.attr = ([], {'color': 'green',  'dash': '0 ',  'dashlist': '8 3 ',  'delete': '1 ',  'edit': '1 ',
                       'fixed': '0 ',  'font': '"helvetica 10 normal roman"',  'highlite': '1 ',
                       'include': '1 ',  'move': '1 ',  'select': '1 ',  'source': '1',  'text': '',
                       'width': '1 '})
        
    if prv_regions[band]:
        prv_regions[band].write('OrionKL_observed_regions_band{0}_private_March2015.reg'.format(band))
    if pub_regions[band]:
        pub_regions[band].write('OrionKL_observed_regions_band{0}_public_March2015.reg'.format(band))
In [113]:
from astropy.io import fits
In [114]:
prv_mask = {band: fits.PrimaryHDU(prv_regions[band].get_mask(orionkl_images[0][0]).astype('int'),
                           header=orionkl_images[0][0].header) for band in bands
            if prv_regions[band]}
pub_mask = {band: fits.PrimaryHDU(pub_regions[band].get_mask(orionkl_images[0][0]).astype('int'),
                           header=orionkl_images[0][0].header) for band in bands
            if pub_regions[band]}
In [115]:
for band in pub_mask:
    pub_mask[band].writeto('public_orionkl_band{0}_almaobs_mask.fits'.format(band), clobber=True)
for band in prv_mask:
    prv_mask[band].writeto('private_orionkl_band{0}_almaobs_mask.fits'.format(band), clobber=True)    
In [117]:
fig = aplpy.FITSFigure(orionkl_images[0])
fig.show_grayscale(stretch='arcsinh')
for band in bands:
    if band in prv_mask:
        fig.show_contour(prv_mask[band], levels=[0.5,1], colors=[band_colors_priv[band]]*2)
    if band in pub_mask:
        fig.show_contour(pub_mask[band], levels=[0.5,1], colors=[band_colors_pub[band]]*2)
INFO:astropy:Auto-setting vmin to  2.393e+02
INFO:astropy:Auto-setting vmax to  3.155e+03
INFO: Auto-setting vmin to  2.393e+02 [aplpy.core]
INFO: Auto-setting vmax to  3.155e+03 [aplpy.core]

More advanced

Now we create a 'hit mask' showing the relative depth of each observed field in each band

In [118]:
hit_mask_public = {band: np.zeros_like(orionkl_images[0][0].data) for band in pub_mask}
hit_mask_private = {band: np.zeros_like(orionkl_images[0][0].data) for band in prv_mask}
from astropy import wcs
mywcs = wcs.WCS(orionkl_images[0][0].header)
In [119]:
for band in bands:
    for row,rad in zip(orionkl, primary_beam_radii):
        shape = Shape('circle', (row['RA'], row['Dec'],np.mean(rad).to(u.deg).value))
        shape.coord_format = 'fk5'
        shape.coord_list = (row['RA'], row['Dec'],np.mean(rad).to(u.deg).value)
        shape.attr = ([], {'color': 'green',  'dash': '0 ',  'dashlist': '8 3 ',  'delete': '1 ',  'edit': '1 ',
                       'fixed': '0 ',  'font': '"helvetica 10 normal roman"',  'highlite': '1 ',
                       'include': '1 ',  'move': '1 ',  'select': '1 ',  'source': '1',  'text': '',
                       'width': '1 '})
        if row['Release date']=='' and row['Band']==band and band in prv_mask:
            (xlo,xhi,ylo,yhi),mask = pyregion_subset(shape, hit_mask_private[band], mywcs) 
            hit_mask_private[band][ylo:yhi,xlo:xhi] += row['Integration']*mask
        if row['Release date']!='' and row['Band']==band and band in pub_mask:
            (xlo,xhi,ylo,yhi),mask = pyregion_subset(shape, hit_mask_public[band], mywcs) 
            hit_mask_public[band][ylo:yhi,xlo:xhi] += row['Integration']*mask
In [120]:
fig = aplpy.FITSFigure(orionkl_images[0])
fig.show_grayscale(stretch='arcsinh')
for band in bands:
    if band in pub_mask:
        fig.show_contour(fits.PrimaryHDU(data=hit_mask_public[band], header=orionkl_images[0][0].header),
                         levels=np.logspace(0,5,base=2, num=6), colors=[band_colors_pub[band]]*6)
    if band in prv_mask:
        fig.show_contour(fits.PrimaryHDU(data=hit_mask_private[band], header=orionkl_images[0][0].header),
                         levels=np.logspace(0,5,base=2, num=6), colors=[band_colors_priv[band]]*6)
INFO:astropy:Auto-setting vmin to  2.393e+02
INFO:astropy:Auto-setting vmax to  3.155e+03
INFO: Auto-setting vmin to  2.393e+02 [aplpy.core]
INFO: Auto-setting vmax to  3.155e+03 [aplpy.core]