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

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

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

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

w = wcs.WCS(hdu.header)

hdu.data.shape

import aplpy

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

from pvextractor import extract_pv_slice
from pvextractor.geometry import Path

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

pv.shape

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

# 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)

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

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')

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

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

# 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')

# 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')