How to display ESO Harps s1d fits file? (Python)
Hi I am looking to plot a spectra on python of a fits file. ESO has a guide on how to display 1D spectra and a code that should work which is the following:
import sys
from astropy.io import fits as pyfits
import numpy as np
#
hdulist = pyfits.open( "your_1d_spectrum_here.fits" )
# print column information
hdulist[1].columns
# get to the data part (in extension 1)
scidata = hdulist[1].data
wave = scidata[0][0]
arr1 = scidata[0][1]
arr2 = scidata[0][2]
# etc.
# where arr1 will contain the data corresponding to the column named: hdulist[1].columns[1]
# where arr2 will contain the data corresponding to the column named: hdulist[1].columns[2]
# etc.
# To plot using maptplotlib:
import matplotlib.pyplot as plt
plt.plot(wave, arr1)
plt.show()
but when putting this code in and using the s1d file I have I get IndexError: list index out of range
The HDUlist.info() consists one one primary HDU with index 0
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 3051 (313088,) float32
and printing this data shows that it is a 1d array: [ 45.182148 71.06376 3.1499128 ... 631.7653 628.2799 632.91364 ]
so I'm clueless on how to plot this spectra as it's only 1d.
Eso also has a more complicated code to plot the spectra as found on the website https://archive.eso.org/cms/eso-data/help/1dspectra.html#Python the file is the 1dspectrum.py which can be downloaded. Once I've inputted my file on there, I get an index error too which can be seen in the gyazo https://gyazo.com/e30a1f39e5f33821b5a5c54f3939a028
Attempted Iguananaut's solution with this as the following code:
import numpy as np
from astropy.io import fits
from astropy.units import u
import matplotlib.pyplot as plt
import specutils
f = fits.open('C:/Users/Rp199/Desktop/J0608_59_harps_2018/HARPS.2018-05-23T23_44_45.005_s1d_A.fits')
f.info()
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 3055 (313093,) float32
Checked the headers and they were the same so I deleted BUNIT and changed CTYPE1 to CUNIT1
del header['BUNIT']
header['CTYPE1'] = header['CUNIT1']
del header['CTYPE1']
Then I tried to use spectrum1D.read
spec = specutils.Spectrum1D.read(f)
spec
which gives me the following error
IORegistryError Traceback (most recent call last)
<ipython-input-28-62aceed26d07> in <module>
----> 1 spec = specutils.Spectrum1D.read(f)
2 spec
~\anaconda3\lib\site-packages\astropy\nddata\mixins\ndio.py in __call__(self, *args, **kwargs)
54
55 def __call__(self, *args, **kwargs):
---> 56 return registry.read(self._cls, *args, **kwargs)
57
58
~\anaconda3\lib\site-packages\astropy\io\registry.py in read(cls, format, *args, **kwargs)
517 fileobj = args[0]
518
--> 519 format = _get_valid_format(
520 'read', cls, path, fileobj, args, kwargs)
521
~\anaconda3\lib\site-packages\astropy\io\registry.py in _get_valid_format(mode, cls, path, fileobj, args, kwargs)
597 if len(valid_formats) == 0:
598 format_table_str = _get_format_table_str(cls, mode.capitalize())
--> 599 raise IORegistryError("Format could not be identified based on the"
600 " file name or contents, please provide a"
601 " 'format' argument.\n"
IORegistryError: Format could not be identified based on the file name or contents, please provide a 'format' argument.
The available formats are:
Format Read Write Auto-identify
----------------- ---- ----- -------------
APOGEE apStar Yes No Yes
APOGEE apVisit Yes No Yes
APOGEE aspcapStar Yes No Yes
ASCII Yes No Yes
ECSV Yes No Yes
HST/COS Yes No Yes
HST/STIS Yes No Yes
IPAC Yes No Yes
JWST s2d Yes No Yes
JWST s3d Yes No Yes
JWST x1d Yes No Yes
SDSS-I/II spSpec Yes No Yes
SDSS-III/IV spec Yes No Yes
Subaru-pfsObject Yes No Yes
iraf Yes No Yes
muscles-sed Yes No Yes
tabular-fits Yes Yes Yes
wcs1d-fits Yes No Yes
I made the assumption that because my file is a s1d.fits file the format would be wcs1d-fits? I tried to look through the headers in Fits viewer to identify a format but couldn't find anything. When putting in wcs1d-fits as the format:
spec = specutils.Spectrum1D.read(f, format="wcs1d-fits")
spec
I am returned this,
TypeError Traceback (most recent call last)
<ipython-input-29-c038104effad> in <module>
----> 1 spec = specutils.Spectrum1D.read(f, format="wcs1d-fits")
2 spec
~\anaconda3\lib\site-packages\astropy\nddata\mixins\ndio.py in __call__(self, *args, **kwargs)
54
55 def __call__(self, *args, **kwargs):
---> 56 return registry.read(self._cls, *args, **kwargs)
57
58
~\anaconda3\lib\site-packages\astropy\io\registry.py in read(cls, format, *args, **kwargs)
521
522 reader = get_reader(format, cls)
--> 523 data = reader(*args, **kwargs)
524
525 if not isinstance(data, cls):
~\anaconda3\lib\site-packages\specutils\io\default_loaders\wcs_fits.py in wcs1d_fits_loader(file_name, spectral_axis_unit, flux_unit, hdu_idx, **kwargs)
63 logging.info("Spectrum file looks like wcs1d-fits")
64
---> 65 with fits.open(file_name, **kwargs) as hdulist:
66 header = hdulist[hdu_idx].header
67 wcs = WCS(header)
~\anaconda3\lib\site-packages\astropy\io\fits\hdu\hdulist.py in fitsopen(name, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
162 raise ValueError(f'Empty filename: {name!r}')
163
--> 164 return HDUList.fromfile(name, mode, memmap, save_backup, cache,
165 lazy_load_hdus, **kwargs)
166
~\anaconda3\lib\site-packages\astropy\io\fits\hdu\hdulist.py in fromfile(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
401 """
402
--> 403 return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,
404 save_backup=save_backup, cache=cache,
405 lazy_load_hdus=lazy_load_hdus, **kwargs)
~\anaconda3\lib\site-packages\astropy\io\fits\hdu\hdulist.py in _readfrom(cls, fileobj, data, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
1052 if not isinstance(fileobj, _File):
1053 # instantiate a FITS file object (ffo)
-> 1054 fileobj = _File(fileobj, mode=mode, memmap=memmap, cache=cache)
1055 # The Astropy mode is determined by the _File initializer if the
1056 # supplied mode was None
~\anaconda3\lib\site-packages\astropy\utils\decorators.py in wrapper(*args, **kwargs)
533 warnings.warn(message, warning_type, stacklevel=2)
534
--> 535 return function(*args, **kwargs)
536
537 return wrapper
~\anaconda3\lib\site-packages\astropy\io\fits\file.py in __init__(self, fileobj, mode, memmap, overwrite, cache)
193 self._open_filename(fileobj, mode, overwrite)
194 else:
--> 195 self._open_filelike(fileobj, mode, overwrite)
196
197 self.fileobj_mode = fileobj_mode(self._file)
~\anaconda3\lib\site-packages\astropy\io\fits\file.py in _open_filelike(self, fileobj, mode, overwrite)
544
545 if mode == 'ostream':
--> 546 self._overwrite_existing(overwrite, fileobj, False)
547
548 # Any "writeable" mode requires a write() method on the file object
~\anaconda3\lib\site-packages\astropy\io\fits\file.py in _overwrite_existing(self, overwrite, fileobj, closed)
442 # The file will be overwritten...
443 if ((self.file_like and hasattr(fileobj, 'len') and fileobj.len > 0) or
--> 444 (os.path.exists(self.name) and os.path.getsize(self.name) != 0)):
445 if overwrite:
446 if self.file_like and hasattr(fileobj, 'truncate'):
~\anaconda3\lib\genericpath.py in exists(path)
17 """Test whether a path exists. Returns False for broken symbolic links"""
18 try:
---> 19 os.stat(path)
20 except (OSError, ValueError):
21 return False
TypeError: stat: path should be string, bytes, os.PathLike or integer, not method
I had a look at a random spectrum downloaded from this page. Indeed it contains only one HDU so you will get an IndexError if you try to access data from an HDU that does not exist.
The code you linked to might be old (it is on an "archive" page after all) and is also merely an example, not necessarily appropriate for plotting data from any FITS file.
The specutils package has many utilities for analyzing and plotting 1D spectra. It can also read spectra from many common FITS formats without having to do too much manually.
Here's what I did. First I opened the file, this one in particular:
>>> f = fits.open('https://www.eso.org/sci/facilities/lasilla/instruments/harps/inst/monitoring/sundata/ceres_2006-05-22_s1d.fits')
>>> f.info()
f = fits.open('https://www.eso.org/sci/facilities/lasilla/instruments/harps/inst/monitoring/sundata/ceres_2006-05-22_s1d.fits')
Looking at the header, especially the first dozen or so header keywords, can give us a lot of clues as to how the data in the file is to be interpreted:
>>> header = f[0].header
>>> header[:12]
SIMPLE = T / file does conform to FITS standard
BITPIX = -32 / number of bits per data pixel
NAXIS = 1 / number of data axes
NAXIS1 = 313237 / length of data axis 1
EXTEND = T / FITS dataset may contain extensions
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
CRPIX1 = 1. / Reference pixel
CRVAL1 = 3781.19 / Coordinate at reference pixel
CDELT1 = 0.01 / Coordinate increment par pixel
CTYPE1 = 'Angstrom' / Units of coordinate
BUNIT = 'Relative Flux' / Units of data values
We can infer from BUNIT = 'Relative Flux' that the single array in the file contains flux values in some unspecified unit (this file is also malformatted in a few ways but I'll get to that later).
It doesn't contain a spectral array, instead the spectrum just starts from 3781.19 Å, at the first pixel (we get this from CRVAL1 and CRPIX1; keep in mind CRPIXn, like other indices in FITS files, uses Fortran-style indexing which starts with 1, as opposed to C/Python which starts at 0, so this is saying the 0-th element of the array is the flux of the spectrum at 3781.19 Å).
It also says CDELT1 = 0.01 / Coordinate increment par pixel so every pixel is +0.01 Å.
So we could do:
>>> import numpy as np
>>> from astropy.units import u
>>> flux = f[0].data
>>> ref = header['CRVAL1']
>>> step = header['CDELT1']
>>> wave = np.arange(ref, ref + (len(flux) * step), step) * u.AA
>>> plt.plot(wave, flux)
Specutils can in principle make this easier though. To read a Spectrum1D from a FITS file--if it's in a supported format--you can do:
>>> from specutils import Spectrum1D
>>> spec = Spectrum1D.read(f)
Unfortunately in this case Specutils failed to read the file and I got:
ValueError: 'Relative Flux' did not parse as unit: At col 0, Relative is not a valid unit. If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html
This is part of what I meant when I wrote earlier that the file is malformatted. "Relative Flux" is not a valid, known name for any unit that would go in the BUNIT header of a FITS file. Since I don't know what the flux units are supposed to be in this case, I just deleted it:
>>> del header['BUNIT']
A further problem with this file is that it isn't using the CTYPE1 keyword correctly. This should be CUNIT1. I don't know why they did this. I again updated the header:
>>> header['CUNIT1'] = header['CTYPE1]
>>> del header['CTYPE1']
Now it works:
>>> spec = Spectrum1D.read(f)
>>> spec
<Spectrum1D(flux=<Quantity [0.0043928 , 0.00500706, 0.00511265, ..., 0.00435138, 0.00433926,
0.00433523]>, spectral_axis=<SpectralAxis [3781.19, 3781.2 , 3781.21, ..., 6913.53, 6913.54, 6913.55] Angstrom>)>
and can be plotted like (just following the specutils documentation):
>>> from astropy.visualization import quantity_support
>>> quantity_support()
>>> fig, ax = plt.subplots()
>>> ax.step(spec.spectral_axis, spec.flux)
>>> fig.show()
