I am tying to get the physical sky coordinates of a given pixel from within a python script. I would like to use astropy's WCS, but I'll do anything from within python.
I have tried these two snippets of code.
from astropy.io import fits
from astropy.wcs import WCS
def astropymethod1(img):
# from http://astropy.readthedocs.org/en/latest/wcs/
w = WCS(img)
lon, lat = w.all_pix2world( 100., 100., 1)
print lon, lat
def astropymethod2(img):
# from http://astropy.readthedocs.org/en/latest/wcs/
hdu = fits.open(img)
w = WCS(hdu[0].header)
lon, lat = w.wcs_pix2world(100., 100., 1)
print lon, lat
The issues are I get an error the first time I try to use WCS and the result is only ever the pixel values I put in.
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
The problem is that you have a multi-extension FITS file. Here's an example session showing how you can get access to the appropriate WCS:
In [1]: from astropy.io import fits
In [2]: h = fits.getheader('SN1415_F625W_1_drz.fits')
In [3]: f = fits.open('SN1415_F625W_1_drz.fits')
In [4]: f
Out[4]:
[<astropy.io.fits.hdu.image.PrimaryHDU at 0x106735490>,
<astropy.io.fits.hdu.image.ImageHDU at 0x106749750>,
<astropy.io.fits.hdu.image.ImageHDU at 0x106751310>,
<astropy.io.fits.hdu.image.ImageHDU at 0x106751d10>,
<astropy.io.fits.hdu.table.BinTableHDU at 0x1067dfdd0>]
In [5]: from astropy import wcs
In [6]: w = wcs.WCS(f[0].header)
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
In [7]: w.wcs.naxis
Out[7]: 2
In [8]: f[0].data
In [9]: w = wcs.WCS(f[1].header)
In [10]: w.wcs.naxis
Out[10]: 2
In [11]: f[1].data
Out[11]:
array([[ 0.01986978, -0.04018363, 0.03330525, ..., 0. ,
0. , 0. ],
[ 0.0695872 , -0.00979143, 0.00147662, ..., 0. ,
0. , 0. ],
[-0.09292094, 0.02481506, -0.01057338, ..., 0. ,
0. , 0. ],
...,
[ 0. , 0. , 0. , ..., 0.02375774,
0.0389731 , 0.03825707],
[ 0. , 0. , 0. , ..., -0.01570918,
-0.01053802, 0.00461219],
[ 0. , 0. , 0. , ..., -0.0638448 ,
-0.0240754 , 0.02679451]], dtype=float32)
In [12]: w.wcs_pix2world(100., 100., 1)
Out[12]: [array(6.113076380801787), array(0.616758775753701)]
So you probably want to redefine your method:
def astropymethod2(img, hduid=1):
# from http://astropy.readthedocs.org/en/latest/wcs/
hdu = fits.open(img)
w = WCS(hdu[hduid].header)
lon, lat = w.wcs_pix2world(100., 100., 1)
print lon, lat