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
(1009, 21, 21)
import aplpy
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.
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
(1009, 28)
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)
# 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]]
F = aplpy.FITSFigure(hdu,slices=[500])
F.show_grayscale()
F.show_lines([endpoints_wcs.T],color='r')
Now for a weirder case, with more endpoints
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')