如何显示 ESO Harps s1d 适合文件? (Python)

How to display ESO Harps s1d fits file? (Python)

您好,我想在拟合文件的 python 上绘制光谱。 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 out of range

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

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

打印此数据显示它是一维数组: [ 45.182148 71.06376 3.1499128 ... 631.7653 628.2799 632.91364 ]

所以我不知道如何绘制这个光谱,因为它只有 1d。

Eso 也有一个更复杂的代码来绘制网站上的光谱 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

尝试了 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   

检查了 headers,它们是一样的,所以我删除了 BUNIT 并将 CTYPE1 更改为 CUNIT1

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 查看器中查看 headers 以识别格式,但找不到任何内容。当输入 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

我查看了从 this page 下载的随机频谱。实际上它只包含一个 HDU,因此如果您尝试从不存在的 HDU 访问数据,您将得到一个 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')

查看 header,尤其是前十几个 header 关键字,可以为我们提供很多关于如何解释文件中数据的线索:

>>> 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' 中推断出文件中的单个数组包含某些未指定单位的通量值(此文件在某些​​方面也存在格式错误,但稍后我会谈到)。

它不包含光谱阵列,而是光谱从第一个像素的 3781.19 Å 开始(我们从 CRVAL1CRPIX1 获得;记住 CRPIXn,与 FITS 文件中的其他索引一样,使用从 1 开始的 Fortran-style 索引,而不是从 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 header 的任何单位,“相对通量”不是有效的已知名称。因为我不知道在这种情况下通量单位应该是什么,所以我只是删除了它:

>>> del header['BUNIT']

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

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