[AstroPy] Loading wavelength from FITS spectrum

David Berry d.berry at jach.hawaii.edu
Thu Jul 4 08:20:09 EDT 2013


Using pyast (http://dsberry.github.io/starlink/pyast.html) you can do
it like this:


import pyfits
import starlink.Atl as Atl

hdulist = pyfits.open( 'test.fit' )
frameset,encoding = Atl.readfitswcs( hdulist )
wcsvalues = frameset.trangrid( 1, hdulist[0].shape[0] )



This assumes your FITS file has only one axis. If you have more axes
(RA, Dec axes for instance), you need additionally to use mapsplit to
split off the mapping for the spectral axis from the complete
n-dimensional FrameSet:

import pyfits
import starlink.Atl as Atl

# The one-based index of the spectral axis
spectralaxis = 1

# Read the FITS headers and get an AST FrameSet describing the
# full WCS.
hdulist = pyfits.open( 'test.fit' )
frameset,encoding = Atl.readfitswcs( hdulist )

# Split off the Mapping that describes just the spectral axis.
wcsaxis,specmap = frameset.mapsplit( spectralaxis )

# Use this Mapping to generate a list of spectral axis values
# corresponding to pixel 1, 2, 3...
wcsvalues = specmap.trangrid( 1, hdulist[0].shape[pixelaxis-1] )



David

On 4 July 2013 10:06, Ian Crossfield <ianc at mpia-hd.mpg.de> wrote:
> The previously suggested solutions work for a linear wavelength
> solution, but things get more complicated for higher-order polynomial
> relations.  In these cases I use the following:
>
> from pyraf import iraf
> import numpy as np
> iraf.wspectext('input.fits', 'output.txt', header=False)
> wavelength_and_flux = np.loadtxt('output.txt')
>
> I agree that it would be convenient to have an AstroPy helper function
> for doing this sort of thing.
>
> -Ian
>
>
>
> On 7/4/13 10:57 AM, Thøger Rivera-Thorsen wrote:
>> I usually go:
>>
>> data, header = pyfits.getdata(fitsfile, header=True)
>> wave = numpy.arange(data.shape[0]) * header['CDELT1'] + header['CRVAL1']
>>
>> That should be about as simple as it gets :-)
>>
>> Cheers,
>> Emil
>>
>>
>>
>> On 2013-07-04 10:43, Éric Depagne wrote:
>>> Le jeudi 4 juillet 2013 à 03:51:20, Gustavo Bragança a écrit :
>>>> Hi,
>>>>
>>> Hi Gustavo,
>>>
>>>> Is there any easy way to load the values of wavelength from a FITS spectrum
>>>> using astropy?
>>>>
>>>> I could obtain the wavelengths by using some values of the header, but I
>>>> was wondering if there is an easy way.
>>>>
>>> I use the following code to get the wavelength and intensity:
>>> def readspec(fitsfile):
>>>           f = pyfits.open(fitsfile)
>>>           cdelt1 = f[0].header['CDELT1']
>>>           crval1 = f[0].header['CRVAL1']
>>>           start = crval1 - f[0].header['CRPIX1'] * cdelt1
>>>           end = start + cdelt1 * f[0].shape[0] - cdelt1/10.
>>>           x = numpy.arange(start, end, cdelt1)
>>>           f.close()
>>>           return x, f[0].data
>>>
>>> I'm not sure it's the best way to do so, but it works.
>>>
>>> HTH.
>>>
>>> Éric.
>>>
>>> Un clavier azerty en vaut deux
>>> ----------------------------------------------------------
>>> Éric Depagne                            eric at depagne.org
>>> _______________________________________________
>>> AstroPy mailing list
>>> AstroPy at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/astropy
>> _______________________________________________
>> AstroPy mailing list
>> AstroPy at scipy.org
>> http://mail.scipy.org/mailman/listinfo/astropy
>
>
> --
> Ian Crossfield
> http://www.mpia-hd.mpg.de/homes/ianc/
> Max-Planck-Institut für Astronomie
> +49-(0)6221 528-405
> Room 308/4
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy



More information about the AstroPy mailing list