Error attempting to save FITS files of VLT/FLAMES coadded stellar spectra

vforvendetta · July 12, 2023, 6:21pm

Hello folks,

In the files here: WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

example_FLAMES.zip contains the 45 FITS files from VLT/FLAMES with stellar spectra. Note that multiple exposures of the same star are saved in different files.
P106_one_FITS_per_spectrum.ipynb is a code that so far reads these 45 FITS files, creates coadded spectra for each of the nearly ~1000 stars, and saves the spectra for each star (wavelength and flux) into ASCII files
Hur2012_table1.dat is catalog information of the stars and is needed to be read for the code

I would like to save the coadded spectrum for each star into its own FITS file. For context, the reason for creating a FITS file for coadded spectra is to read the FITS files into a Python code (use_ccf_functions.py) which calculates radial velocities for each source through cross-correlation.

Here is roughly how the ccf code works:

in the script use_ccf_functions one can define (commenting and uncommenting) the mask to be used, the data directory for input and the one for output (where the table and the figures will be placed, and some checking is done not to overwrite), and the instrument (in this case, FLAMES).
in MODULES/calc_rv_new_chapman one can set the number of CPUS (for mp, pathos.multiprocessing), the step in wavelength, RV, and vsini. This module also contains the function I am using to read the FLAMES FITS files: ReadSpec_FLAMES_fits(fits_file):

My current python code (P106_one_FITS_per_spectrum.ipynb) attempts to save the coadded spectrum into a single FITS file named after the source ID. However, when I try to have this new FITS file (example_FLAMES_test/new.fits) read by the use_ccf_functions.py, it cannot read the format properly (attempt_output.txt is the terminal output).

use_ccf_functions.py CAN read in the 45 FITS file outputted by the GIRAFFE pipeline properly, but as you may guess, the FLAMES FITS files contain multiple IDs rather than ideally the coadded spectra of source ID.

Could anyone see what I may be missing in how I creating the new FITS files such that it cannot be read by the CCF code at the moment? Any assistance would be greatly appreciated.

saimn · July 14, 2023, 3:35pm

Could you show the code that write the file and the associated error ?

vforvendetta · July 14, 2023, 5:36pm

@saimn Hello, the code that attempts to write the FITS file with the coadded spectra is P106_one_FITS_per_spectrum.ipynb. To see the whole code, please refer to the download link in the original post, but here is the snippet of the code where I attempt to write the FITS file.

# Create new FITS file

from astropy.io import fits

# Find the first FITS file where a fibre containing the spectrum for 2620 is found
for spectra_path in list_spectra_dirs:
    f = fits.open(spectra_path)
    specdata = f[0].data

    for i in range(len(f[1].data['OBJECT'])):
        if f[1].data['OBJECT'][i] == "2620":
            # Use the primary HDU of this FITS file as the template for the new FITS file
            template_hdulist = fits.open(spectra_path)
            break

# Create the new FITS file
hdulist = fits.HDUList()

# Primary HDU
hdulist.append(template_hdulist[0])

# First extension HDU with the coadded spectra
coadded_data = np.array([spectrum.flux.value for spectrum in s.values()])
coadded_spec = Spectrum1D(spectral_axis=new_disp_grid, flux=np.sum(coadded_data, axis=0) * u.adu)
coadded_hdu = fits.ImageHDU(data=coadded_spec.flux.value, header=template_hdulist[1].header)
hdulist.append(coadded_hdu)

hdulist.writeto('new.fits', overwrite=True)

Here is the associated error when I try to read ‘new.fits’ into the ccf python code (use_ccf_functions.py). This may also be found in the file, ‘attempt_output.txt’

['new.fits']


1 / 1 still not in output table



Starting CCF and Vsin(i) procedure...
Traceback (most recent call last):
  File "/home/nick/P106_Research/wd/2020_04_28_version/use_ccf_functions.py", line 89, in <module>
    c_b, ad_sig = master_make_ccf_vsini(files_without_ccf[i], flag='fits', mask_type=mask_name, instrument=instrument)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/nick/P106_Research/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 1266, in master_make_ccf_vsini
    wlen_obs, flux_obs, bary_correc = ReadSpec_FLAMES_fits(file_location)  
                                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/nick/P106_Research/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 441, in ReadSpec_FLAMES_fits
    star_name = data['OBJECT'][ispec]
                ~~~~^^^^^^^^^^
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

dhomeier · July 14, 2023, 11:10pm

vforvendetta:

  File "/home/nick/P106_Research/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 441, in ReadSpec_FLAMES_fits
    star_name = data['OBJECT'][ispec]
                ~~~~^^^^^^^^^^
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

That error cannot originate from the code snippet you posted above, as there is no star_name variable used there anywhere.
But the IndexError points to attempting to index a plain numpy array as if it were a Table with a column ‘OBJECT’; quite likely you have read your data from an ImageHDU rather than a BinaryTableHDU – I’d check if you are reading from the right HDU at that point.

vforvendetta · July 17, 2023, 6:36pm

@dhomeier You are correct - I was reading the data from an ImageHDU. Here is the result when I run hdulist.info().

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     532   (133, 6181)   float32   
  1  FIBER_SETUP    1 ImageHDU         9   (4286,)   float64

The star_name variable is from the module file (calc_rv_new_chapman.py) of the ccf code (use_ccf_functions.py). Here the relevant snippet/function in the code that is meant to read in the FITS file I am trying to create:

def ReadSpec_FLAMES_fits(fits_file):

        """Function to read in fits file of high
        resolution spectra and then return wlen 
        and flux arrays"""

        # open fits file
        image=pf.open(fits_file)
        
        # retrieve flux and header
        data = image[0].header
        header = image[1].data 
        
        # get object data
        ispec = 1
        #wlen, flux, err = data[0]
        star_name = data['OBJECT'][ispec]
        flux =  image['FLUX_SKY_SUB'].data[:,ispec]
        err = np.sqrt(flux)*0.01      #dummy errors, need to be loaded from science_rbnerrors_0000.fits
        
        # get coords
        ra_val = data['RA'][ispec]
        dec_val = data['DEC'][ispec]
        
        wlen = (header['ESO PRO REBIN WLEN MIN'] + np.arange(header['ESO PRO REBIN NX'])*header['ESO PRO REBIN LAMBDA STEP']) * 10.
        
        print ("Min. - Max. wlen: %i - %i \n" % (wlen[0], wlen[-1]))
        
        # file type to distinguish between        
        file_type = header['HIERARCH ESO PRO CATG']
        
        # date of observations
        mjd_obs = header['MJD-OBS']
                  
        # barycentric correction already applied for FEROS
        bary_correc = 0
                  
        image.close()
        
        # define star_name again with file type and obs date
        star_name = "%s_%s" % (star_name.replace(" ", ""), \
        file_type)
        
        print (star_name, ra_val, dec_val, mjd_obs, wlen, flux, bary_correc     )
        
        return star_name, ra_val, dec_val, mjd_obs, wlen, flux, bary_correc

Since the primary HDU is also not a BinaryTableHDU, I believe the manner in which I am attempting to create FITS files that contain the coadded spectrum for each star (new.fits) is not correct. Is it possible for a second set of eyes to read over the files where I create newfits ( P106_one_FITS_per_spectrum.ipynb) so I could reconcile having it read correctly by the ReadSpec_FLAMES_fits in my ccf code?

dhomeier · August 3, 2023, 3:13pm

The FITS standard does not allow anything but an ImageHDU as primary HDU, so the standard procedure for storing tabular data is to have a basically empty PrimaryHDU and write the tables to the first or further extensions.
It would probably help if you could directly attach your scripts and perhaps error logs, as extracting everything from a huge zip file may be too time consuming for most contributors.

vforvendetta · August 5, 2023, 10:26pm

Very well. If it helps with troubleshooting, please let me know if you would like the full file/scripts anyway.

This is a code snippet from the file where I am creating a new fits file of the coadded spectra of one star.

# Create new FITS file

# convert to numpy arrays for saving to FITS file
wavelength = coadded_spec.spectral_axis.value
flux = coadded_spec.flux.value

# Create Primary HDU
primary_hdu = fits.PrimaryHDU()
primary_hdu.header = f[0].header  # f from the original file reading, assuming it's the last one read

# Create Binary Table HDU
# Using header data of original Binary Table HDU
table_hdu = fits.BinTableHDU.from_columns(columns=f[1].columns, header=f[1].header)

# Create Image HDUs
# Assuming they're zero-filled like in the original FITS files.
# If you have real data to put here, replace np.zeros_like(flux) with your data.
skyback_hdu = fits.ImageHDU(data=np.zeros_like(flux), header=f[2].header)
flux_sky_sub_hdu = fits.ImageHDU(data=flux, header=f[3].header)  # replace with co-added flux

# Construct HDU list and write to file
hdulist = fits.HDUList([primary_hdu, table_hdu, skyback_hdu, flux_sky_sub_hdu])
hdulist.writeto('new.fits', overwrite=True)

I am attempting to read this file (news.fits) into the module file (calc_rv_new_chapman.py) of the ccf code (use_ccf_functions.py). Here the relevant snippet/function in the code that is meant to read in the FITS file I am trying to create:

def ReadSpec_FLAMES_fits(fits_file):

        """Function to read in fits file of high
        resolution spectra and then return wlen 
        and flux arrays"""

        # open fits file
        image=pf.open(fits_file)
        
        # retrieve flux and header
        data = image[0].data
        header=image[1].header
        

        # get object data
        ispec = 1
        #wlen, flux, err = data[0]

        star_name = data['OBJECT'][ispec]
        flux =  image['FLUX_SKY_SUB'].data[:,ispec]
        err = np.sqrt(flux)*0.01                                #dummy errors, need to be loaded from science_rbnerrors_0000.fits
        
        # get coords
        ra_val = data['RA'][ispec]
        dec_val = data['DEC'][ispec]

        wlen = (header['ESO PRO REBIN WLEN MIN'] + np.arange(header['ESO PRO REBIN NX'])*header['ESO PRO REBIN LAMBDA STEP']) * 10.

        
        print ("Min. - Max. wlen: %i - %i \n" % (wlen[0], wlen[-1]))
        
        # file type to distinguish between        
        file_type = header['HIERARCH ESO PRO CATG']
        
        # date of observations
        mjd_obs = header['MJD-OBS']
                  
        # barycentric correction already applied for FEROS
        bary_correc = 0
                  
        image.close()
        
        # define star_name again with file type and obs date
        star_name = "%s_%s" % (star_name.replace(" ", ""), \
        file_type)
        
        print (star_name, ra_val, dec_val, mjd_obs, wlen, flux, bary_correc     )
        
        return star_name, ra_val, dec_val, mjd_obs, wlen, flux, bary_correc

Here is the error output:

python use_ccf_functions.py
['F1-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC1.fits', 'F1-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC2.fits', 'F1-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC3.fits', 'F1-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC4.fits', 'F1-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC5.fits', 'F1-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC6.fits', 'F2-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC7.fits', 'F2-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC8.fits', 'F2-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC9.fits', 'F2-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC10.fits', 'F2-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC11.fits', 'F2-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC12.fits', 'F3-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC13.fits', 'F3-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC14.fits', 'F3-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC15.fits', 'F3-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC16.fits', 'F3-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC17.fits', 'F3-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC18.fits', 'F4-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC19.fits', 'F4-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC20.fits', 'F4-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC21.fits', 'F4-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC22.fits', 'F4-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC23.fits', 'F4-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC24.fits', 'F5-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC25.fits', 'F5-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC26.fits', 'F5-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC27.fits', 'F5-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC28.fits', 'F5-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC29.fits', 'F5-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC30.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC31.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC32.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC33.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC34.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC35.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC36.fits', 'F6-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC37.fits', 'F7-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC38.fits', 'F7-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC39.fits', 'F7-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC40.fits', 'F7-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC41.fits', 'F7-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC42.fits', 'F7-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC43.fits', 'F8-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC44.fits', 'F8-Carina-FLAMES_SCIENCE_RBNSPECTRA_CRCM=PYCOSMIC45.fits']


45 / 45 still not in output table



Starting CCF and Vsin(i) procedure...
Traceback (most recent call last):
  File "/home/nick/wd/2020_04_28_version/use_ccf_functions.py", line 88, in <module>
    c_b, ad_sig = master_make_ccf_vsini(files_without_ccf[i], flag='fits', mask_type=mask_name, instrument=instrument)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/nick/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 1269, in master_make_ccf_vsini
    wlen_obs, flux_obs, bary_correc = ReadSpec_FLAMES_fits(file_location)  
                                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/nick/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 443, in ReadSpec_FLAMES_fits
    star_name = data['OBJECT'][ispec]
                ~~~~^^^^^^^^^^
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

I believe the issue is how I am creating the new fits file in the first place (hence why I intend to include every relevant code/script for context to the issue). What is my mistake with creating the new FITS file?

dhomeier · August 9, 2023, 9:55pm

As I guessed above, you are trying to reference data['OBJECT'] as if data was a table or structured array, when it is the data section from the (Image) HDU 0, a plain ndarray.
Is the object information in image[1].data perhaps?

vforvendetta · August 11, 2023, 8:36pm

Even when I run with image[1].data, it seems neither extension of this FITS file is structured as a table or structured array.

  File "/home/nick/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 1269, in master_make_ccf_vsini
    wlen_obs, flux_obs, bary_correc = ReadSpec_FLAMES_fits(file_location)  
                                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/nick/wd/2020_04_28_version/./MODULES/calc_rv_new_chapman.py", line 443, in ReadSpec_FLAMES_fits
    star_name = data['OBJECT'][ispec]
                ~~~~^^^^^^^^^^
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

dhomeier · August 11, 2023, 9:39pm

Indeed, could have seen this already from the file info you posted earlier:

vforvendetta:

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     532   (133, 6181)   float32   
  1  FIBER_SETUP    1 ImageHDU         9   (4286,)   float64

This may need someone familiar with the FLAMES file structure, or some background documentation why you expect the files to have a data section with that information. OBJECT names or IDs may commonly be stored just in the header (although in your case it seems you need arrays of names and coordinates etc., which points to a table again).

With the file above probably the next step

would fail just as well, as it tries to access a HDU of Name (EXTNAME) FLUX_SKY_SUB, which does not exist in the listing.

Topic		Replies	Views
Best practices for writing spectra files Astropy specutils , fits	1	324	November 10, 2021
Problem with annotation script for FIT files Astropy astropy , question	0	314	June 2, 2023
Fits to Jpeg Conversion Using Funpack (cfitsio) package OpenAstronomy fits , question	0	353	January 13, 2023
Can't open fits files using astropy, "fits files are empty or corrupted" Astropy question	3	173	January 24, 2024
Fits file and WCS issue Astropy	1	419	June 10, 2023

Error attempting to save FITS files of VLT/FLAMES coadded stellar spectra

Related Topics