如何显示ESO s1d fits文件？(Python)

Question

嗨，我想在python上绘制一个fits文件的光谱。ESO有一个关于如何显示一维谱的指南和一个应该可以工作的代码，如下所示：

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()

但是，当将这段代码放入并使用我有的s1d文件时，我得到IndexError: list index超出范围

HDUlist.info()包含一个索引为0的主HDU

No.    Name      Ver    Type      Cards   Dimensions   Format
0  PRIMARY       1 PrimaryHDU    3051   (313088,)   float32   

打印这个数据显示它是一个一维数组: 631.7653 628.2799 632.91364 ... 45.182148 71.06376 3.1499128

所以我不知道如何绘制这个谱，因为它只是一维的。

Eso还有一个更复杂的代码来绘制光谱，可以在网站https://archive.eso.org/cms/eso-data/help/1dspectra.html#Python上找到，该文件是1dspectrum.py，可以下载。一旦我在那里输入了我的文件，我也会得到一个索引错误，这可以在gyazo https://gyazo.com/e30a1f39e5f33821b5a5c54f3939a028中看到

尝试Iguananaut的解决方案，代码如下：

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   

检查了报头，它们是相同的，所以我删除了CUNIT1并将CTYPE1更改为BUNIT

del header['BUNIT']
header['CTYPE1'] = header['CUNIT1']
del header['CTYPE1']

然后我试着用spectrum1D.read

spec = specutils.Spectrum1D.read(f)
spec

这给出了下面的错误

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

我假设因为我的文件是s1d.fits文件，所以格式应该是wcs1d-fits？我试图在Fits查看器中查看标题，以确定格式，但什么也找不到。当输入wcs1d-fits作为格式时：

spec = specutils.Spectrum1D.read(f, format="wcs1d-fits")
spec

我得到了这个，

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

Answer 1

我看了一下从this page下载的随机频谱。实际上，它只包含一个硬盘单元，所以如果你试图从一个不存在的硬盘单元访问数据，你会得到一个IndexError。

您链接到的代码可能很旧(毕竟是在“存档”页面上)，而且也只是一个示例，不一定适合从任何FITS文件中绘制数据。

specutils软件包具有许多用于分析和绘制一维谱的实用程序。它还可以从许多常见的FITS格式读取光谱，而不必做太多手动操作。

这是我所做的。首先，我打开了文件，特别是这个文件：

>>> 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')

查看标题，尤其是头12个左右的标题关键字，可以为我们提供许多关于如何解释文件中的数据的线索：

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

我们可以从BUNIT = 'Relative Flux'中推断出文件中的单个数组包含一些未指定单位的流量值(该文件在某些方面也有格式错误，但我稍后会介绍)。

它不包含频谱数组，而是频谱从第一个像素开始(我们从CRVAL1和CRPIX1获得；请记住，CRPIXn，就像FITS文件中的其他索引一样，使用从1开始的Fortran风格的索引，而不是从0开始的C/Python，所以这就是说，数组的第0个元素是谱的通量在3781.19ä)。

它还写着CDELT1 = 0.01 / Coordinate increment par pixel，所以每个像素都是+0.01µ。

所以我们可以这样做：

>>> 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原则上可以让这一切变得更容易。要从FITS文件中读取Spectrum1D --如果它是受支持的格式--您可以这样做：

>>> from specutils import Spectrum1D
>>> spec = Spectrum1D.read(f)

不幸的是，在这种情况下，Specutils无法读取文件，我得到：

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

这就是我之前写的文件格式错误的部分意思。对于FITS文件的BUNIT头中的任何单位，“相对通量”都不是一个有效的已知名称。因为我不知道本例中的通量单位是什么，所以我就把它删除了：

>>> del header['BUNIT']

这个文件的另一个问题是它没有正确使用CTYPE1关键字。这应该是CUNIT1。我不知道他们为什么要这么做。我再次更新了头文件：

>>> header['CUNIT1'] = header['CTYPE1]
>>> del header['CTYPE1']

现在它可以工作了：

>>> 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>)>

并且可以像这样绘制(只需遵循specutils文档)：

>>> from astropy.visualization import quantity_support
>>> quantity_support()
>>> fig, ax = plt.subplots()
>>> ax.step(spec.spectral_axis, spec.flux)
>>> fig.show()

​

​