[AstroPy] [pywcs] Unit question

Ole Streicher astropy at liska.ath.cx
Fri Apr 1 09:03:11 EDT 2011

Dear list,

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[0], wcs.wcs.crval[1], l0],
          [ wcs.wcs.crval[0], wcs.wcs.crval[1], l1] ]

pix_coords = numpy.copy(wcs.wcs_sky2pix(edges, 0))
h = pyfits.ImageHDU(
      data = hdu.data[int(pix_coords[0][2]):int(pix_coords[1][2])+1,:, :],
      header = hdu.header)

# Here comes the funny stuff
start = [0, 0, int(pix_coords[0][2])]
h.header['CRVAL3'] = wcs.wcs_pix2sky([start], 1)[0][2]

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 mailing list