[AstroPy] [pywcs] Unit question
astropy at liska.ath.cx
Fri Apr 1 09:03:11 EDT 2011
I am a bit struggling with the PyWCS coordinate transformation. I want
to extract a subset of a data cube and keep the WCS information valid.
In my case, the third coordinate is the wavelength, in "nm" units:
(error-checking etc. is removed here)
import pyfits, pywcs
l0 = 6e-7
l1 = 7e-7
wcs = pywcs.WCS(hdu.header)
edges = [ [ wcs.wcs.crval, wcs.wcs.crval, l0],
[ wcs.wcs.crval, wcs.wcs.crval, l1] ]
pix_coords = numpy.copy(wcs.wcs_sky2pix(edges, 0))
h = pyfits.ImageHDU(
data = hdu.data[int(pix_coords):int(pix_coords)+1,:, :],
header = hdu.header)
# Here comes the funny stuff
start = [0, 0, int(pix_coords)]
h.header['CRVAL3'] = wcs.wcs_pix2sky([start], 1)
The problem here is that "pix2sky" returns the wavelength in meters
while I want to put it in the unit that was used in the original file.
Is there an easy way to convert this?
More information about the AstroPy