Skip to content

Commit

Permalink
Merge pull request #17 from lsst/tickets/DM-33326
Browse files Browse the repository at this point in the history 
DM-33326: Fix round-tripping of spectractor spectra
  • Loading branch information
@mfisherlevine
mfisherlevine committed Feb 17, 2022
2 parents dd7be1d + a95638f commit 2c6a49a
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 40 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,3 +26,5 @@ spectractor/version.py
.sconsign.dblite
config.log
tests/.tests/
*__pycache__/

121 changes: 82 additions & 39 deletions spectractor/extractor/spectrum.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
import matplotlib.pyplot as plt
import numpy as np
import os
import astropy

from spectractor import parameters
from spectractor.config import set_logger, load_config, apply_rebinning_to_parameters
Expand All @@ -20,6 +21,22 @@
from spectractor.simulation.throughput import TelescopeTransmission


fits_mappings = {'date_obs': 'DATE-OBS',
'expo': 'EXPTIME',
'airmass': 'AIRMASS',
'disperser_label': 'GRATING',
'units': 'UNIT2',
'rotation_angle': 'ROTANGLE',
'dec': 'DEC',
'hour_angle': 'HA',
'temperature': 'OUTTEMP',
'pressure': 'OUTPRESS',
'humidity': 'OUTHUM',
'lambda_ref': 'LBDA_REF',
'parallactic_angle': 'PARANGLE',
}


class Spectrum:
""" Class used to store information and methods relative to spectra and their extraction.
Expand Down Expand Up @@ -123,7 +140,9 @@ class Spectrum:
Size of the spectrogram along the y axis.
"""

def __init__(self, file_name="", image=None, order=1, target=None, config="", fast_load=False):
def __init__(self, file_name="", image=None, order=1, target=None, config="", fast_load=False,
spectrogram_file_name_override=None,
psf_file_name_override=None,):
""" Class used to store information and methods relative to spectra and their extraction.
Parameters
Expand Down Expand Up @@ -203,9 +222,17 @@ def __init__(self, file_name="", image=None, order=1, target=None, config="", fa
self.lambdas_order2 = None
self.data_order2 = None
self.err_order2 = None
self.dec = None
self.hour_angle = None
self.temperature = None
self.pressure = None
self.humidity = None
self.parallactic_angle = None
self.filename = file_name
if file_name != "":
self.load_spectrum(file_name)
self.load_spectrum(file_name,
spectrogram_file_name_override=spectrogram_file_name_override,
psf_file_name_override=psf_file_name_override)
if image is not None:
self.header = image.header
self.date_obs = image.date_obs
Expand Down Expand Up @@ -484,6 +511,18 @@ def save_spectrum(self, output_file_name, overwrite=False):
'third column the corresponding errors.'
hdu1 = fits.PrimaryHDU()
hdu1.header = self.header

for attribute, header_key in fits_mappings.items():
try:
value = getattr(self, attribute)
except AttributeError:
print(f"Failed to get {attribute}")
continue
if isinstance(value, astropy.coordinates.angles.Angle):
value = value.degree
hdu1.header[header_key] = value
print(f"Set header key {header_key} to {value} from attr {attribute}")

hdu1.header["EXTNAME"] = "SPECTRUM"
hdu2 = fits.ImageHDU()
hdu2.header["EXTNAME"] = "SPEC_COV"
Expand Down Expand Up @@ -553,14 +592,21 @@ def save_spectrogram(self, output_file_name, overwrite=False):
hdu.writeto(output_file_name, overwrite=overwrite)
self.my_logger.info('\n\tSpectrogram saved in %s' % output_file_name)

def load_spectrum(self, input_file_name):
def load_spectrum(self, input_file_name, spectrogram_file_name_override=None,
psf_file_name_override=None):
"""Load the spectrum from a fits file (data, error and wavelengths).
Parameters
----------
input_file_name: str
Path to the input fits file
spectrogram_file_name_override : str
Manually specify a path to the spectrogram file.
psf_file_name_override : str
Manually specify a path to the psf file.
Examples
--------
>>> s = Spectrum()
Expand All @@ -575,48 +621,33 @@ def load_spectrum(self, input_file_name):
self.data = raw_data[1]
if len(raw_data) > 2:
self.err = raw_data[2]
if self.header['DATE-OBS'] != "":
self.date_obs = self.header['DATE-OBS']
if self.header['EXPTIME'] != "":
self.expo = self.header['EXPTIME']
if self.header['AIRMASS'] != "":
self.airmass = self.header['AIRMASS']
if self.header['GRATING'] != "":
self.disperser_label = self.header['GRATING']
if self.header['TARGET'] != "":
self.target = load_target(self.header['TARGET'], verbose=parameters.VERBOSE)

# set the simple items from the mappings. More complex items, i.e.
# those needing function calls, follow
for attribute, header_key in fits_mappings.items():
if (item := self.header.get(header_key)) is not None:
setattr(self, attribute, item)
print(f'set {attribute} to {item}')
else:
print(f'Failed to set spectrum attribute {attribute} using header {header_key}')

# set the more complex items by hand here
if target := self.header.get('TARGET'):
self.target = load_target(target, verbose=parameters.VERBOSE)
self.lines = self.target.lines
if self.header['UNIT2'] != "":
self.units = self.header['UNIT2']
if self.header['ROTANGLE'] != "":
self.rotation_angle = self.header['ROTANGLE']
if self.header['TARGETX'] != "" and self.header['TARGETY'] != "":
self.x0 = [self.header['TARGETX'], self.header['TARGETY']]
if self.header['D2CCD'] != "":
parameters.DISTANCE2CCD = float(self.header["D2CCD"])
if 'DEC' in self.header and self.header['DEC'] != "":
self.dec = self.header['DEC']
if 'RA' in self.header and self.header['HA'] != "":
self.hour_angle = self.header['HA']
if 'OUTTEMP' in self.header and self.header['OUTTEMP'] != "":
self.temperature = self.header['OUTTEMP']
if 'OUTPRESS' in self.header and self.header['OUTPRESS'] != "":
self.pressure = self.header['OUTPRESS']
if 'OUTHUM' in self.header and self.header['OUTHUM'] != "":
self.humidity = self.header['OUTHUM']
if self.header['LBDA_REF'] != "":
self.lambda_ref = self.header['LBDA_REF']
if 'PARANGLE' in self.header and self.header['PARANGLE'] != "":
self.parallactic_angle = self.header['PARANGLE']
if 'CCDREBIN' in self.header and self.header['CCDREBIN'] != "":
if parameters.CCD_REBIN != self.header['CCDREBIN']:
if (targetx := self.header.get('TARGETX')) and (targety := self.header.get('TARGETY')):
self.x0 = [targetx, targety] # should be a tuple not a list
if rebin := self.header.get('CCDREBIN'):
if parameters.CCD_REBIN != rebin:
raise ValueError("Different values of rebinning parameters between config file and header. Choose.")
parameters.CCD_REBIN = rebin
if dist := self.header.get('D2CCD'):
parameters.DISTANCE2CCD = float(dist)

self.my_logger.info('\n\tLoading disperser %s...' % self.disperser_label)
self.disperser = Hologram(self.disperser_label, D=parameters.DISTANCE2CCD,
data_dir=parameters.DISPERSER_DIR, verbose=parameters.VERBOSE)
self.my_logger.info('\n\tSpectrum loaded from %s' % input_file_name)
spectrogram_file_name = input_file_name.replace('spectrum', 'spectrogram')
if parameters.OBS_OBJECT_TYPE == "STAR":
self.adr_params = [self.dec, self.hour_angle, self.temperature,
self.pressure, self.humidity, self.airmass]
Expand All @@ -632,13 +663,25 @@ def load_spectrum(self, input_file_name):
self.target.image_y0 = float(hdu_list["ORDER0"].header["IM_Y0"])
else:
self.cov_matrix = np.diag(self.err ** 2)

# original, hard-coded spectrogram/table relative paths
spectrogram_file_name = input_file_name.replace('spectrum', 'spectrogram')
psf_file_name = input_file_name.replace('spectrum.fits', 'table.csv')

# for LSST-DM supplied filenames
if spectrogram_file_name_override and psf_file_name_override:
self.fast_load = False
spectrogram_file_name = spectrogram_file_name_override
psf_file_name = psf_file_name_override
self.my_logger.info(f'Applying spectrogram filename override {spectrogram_file_name}')
self.my_logger.info(f'Applying psf filename override {psf_file_name}')

if not self.fast_load:
self.my_logger.info(f'\n\tLoading spectrogram from {spectrogram_file_name}...')
if os.path.isfile(spectrogram_file_name):
self.load_spectrogram(spectrogram_file_name)
else:
raise FileNotFoundError(f"Spectrogram file {spectrogram_file_name} does not exist.")
psf_file_name = input_file_name.replace('spectrum.fits', 'table.csv')
self.my_logger.info(f'\n\tLoading PSF from {psf_file_name}...')
if os.path.isfile(psf_file_name):
self.load_chromatic_psf(psf_file_name)
Expand Down
4 changes: 3 additions & 1 deletion spectractor/extractor/targets.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,9 @@ def load_spectra(self):
else:
self.wavelengths.append(data.wave)
self.spectra.append(data.flux)
elif 'HD' in self.label or self.label in parameters.STAR_NAMES: # it is a star
# TODO DM-33731: the use of self.label in parameters.STAR_NAMES:
# below works for running but breaks a test so needs fixing for DM
elif 'HD' in self.label: # or self.label in parameters.STAR_NAMES: # it is a star
self.emission_spectrum = False
self.hydrogen_only = False
self.lines = Lines(ATMOSPHERIC_LINES + HYDROGEN_LINES + STELLAR_LINES,
Expand Down

0 comments on commit 2c6a49a

Please sign in to comment.