In [1]:

    

url = 'http://files.figshare.com/1213488/iras05358_12co21.fits'


    

In [2]:

    

from astropy import log
log.setLevel('ERROR')


    

In [3]:

    

from astropy.utils.data import download_file
from astropy.io import fits
from astropy import wcs
from astropy import coordinates
from astropy import units as u
import os


    

In [4]:

    

f = download_file(url, cache=True)
fitsfile = fits.open(f)
hdu = fitsfile[0]


    

In [5]:

    

w = wcs.WCS(hdu.header)


    

In [6]:

    

hdu.data.shape


    

    
Out[6]:





(1009, 21, 21)

In [7]:

    

import aplpy


    

In [8]:

    

F = aplpy.FITSFigure(hdu,slices=[500])
F.show_grayscale()


    

    




/Users/adam/repos/aplpy/aplpy/labels.py:432: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if self.coord == x or self.axis.apl_tick_positions_world[ipos] > 0:
/Users/adam/repos/aplpy/aplpy/normalize.py:115: RuntimeWarning: invalid value encountered in less
  negative = result < 0.






    

In [9]:

    

from pvextractor import extract_pv_slice
from pvextractor.geometry import Path


    

In [10]:

    

endpoints = [(0,0),(20,20)]
xy = Path(endpoints)
pv = extract_pv_slice(hdu, xy)


    

In [11]:

    

pv.shape


    

    
Out[11]:





(1009, 28)

In [12]:

    

F2 = aplpy.FITSFigure(pv)
F2.show_grayscale(aspect='auto', stretch='arcsinh')


    

    




/Users/adam/anaconda/envs/astropy27/lib/python2.7/site-packages/numpy/ma/core.py:785: RuntimeWarning: invalid value encountered in greater_equal
  return umath.absolute(a) * self.tolerance >= umath.absolute(b)






    

In [13]:

    

# use the figure WCS to overplot a line on the "flat" image
endpoints_wcs = w.sub([wcs.WCSSUB_CELESTIAL]).wcs_pix2world([[x,y] for x,y in endpoints], 0)
print(endpoints_wcs)


    

    




[[ 84.83347466  35.73221637]
 [ 84.76501436  35.78777385]]

In [14]:

    

F = aplpy.FITSFigure(hdu,slices=[500])
F.show_grayscale()
F.show_lines([endpoints_wcs.T],color='r')


    

In [15]:

    

endpoints = [(8,8),(8,15),(15,15),(15,8),(3,5)]
endpoints_wcs = F._wcs.all_pix2world([[x,y,0] for x,y in endpoints], 0)[:,:2]
F = aplpy.FITSFigure(hdu,slices=[500])
F.show_grayscale()
F.show_lines([endpoints_wcs.T],color='r')


    

In [16]:

    

# First, the same approach as above: use the pixel coordinates
path = Path(endpoints)
pv = extract_pv_slice(hdu, path)


    

In [17]:

    

F2 = aplpy.FITSFigure(pv)
F2.show_grayscale(aspect=pv.shape[1]/float(pv.shape[0]), stretch='arcsinh')


    

In [18]:

    

# Second, we use a pvextractor utility WCSPath, which creates a Path object from the WCS
# Start with an array of coordinates:
ep_wcs = coordinates.ICRS([84.806102,84.8061045,84.782138,84.78214175,84.82321362]*u.deg,
                          [35.75444401,  35.77388845, 35.77388811, 35.75444367, 35.74610801]*u.deg)
path = Path(ep_wcs)
pv = extract_pv_slice(hdu, path)
F3 = aplpy.FITSFigure(pv)
F3.show_grayscale(aspect='auto', stretch='arcsinh')


    

In [19]:

    

# Finally, a demonstration of averaging over finite width
path = Path(endpoints, width=3)
pv = extract_pv_slice(hdu, path)
F3 = aplpy.FITSFigure(pv)
F3.show_grayscale(aspect='auto', stretch='arcsinh')