[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?

Regards

Ole



More information about the AstroPy mailing list