From d58f6527fbe1ec6abcbb1b5396dd0cc7233b88e7 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 18 Jan 2023 15:22:46 -0500 Subject: [PATCH 01/31] Working on line pixel refinement --- specreduce/wavelength_calibration.py | 198 +++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100644 specreduce/wavelength_calibration.py diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py new file mode 100644 index 00000000..09de76b7 --- /dev/null +++ b/specreduce/wavelength_calibration.py @@ -0,0 +1,198 @@ +from astropy.modeling.models import Linear1D, Gaussian1D +import astropy.units as u +from specutils import SpectralRegion +from specutils.fitting import fit_lines + + +""" +Starting from Kyle's proposed pseudocode from https://github.com/astropy/specreduce/pull/152 +""" + +class CalibrationLine(): + + def __init__(input_spectrum, wavelength, pixel, refinement_method=None, + refinement_kwargs={}): + """ + input_spectrum: `~specutils.Spectrum1D` + A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. + wavelength: `astropy.units.Quantity` + The rest wavelength of the line. + pixel: int + Initial guess of the integer pixel location of the line center in the spectrum. + refinement_method: str, optional + None: use the actual provided pixel values as anchors for the fit without refining + their location. + 'gradient': Simple gradient descent/ascent to nearest min/max value. Direction + defined by ``direction`` in refinement_kwargs. + 'min'/'max': find min/max within pixel range provided by ``range`` kwarg. + 'gaussian': fit gaussian to spectrum within pixel range provided by ``range`` kwarg. + refinement_kwargs: dict, optional + Keywords to set the parameters for line location refinement. Distance on either side + of line center to include on either side of gaussian fit is set by int ``range``. + Gradient descent/ascent is determined by setting ``direction`` to 'min' or 'max'. + """ + self.input_spectrum = input_spectrum + self.wavelength = wavelength + self.pixel = pixel + self.refinement_method = refinement_method + self.refinement_kwargs = refinement_kwargs + + if self.refinement_method in ("gaussian", "min", "max"): + if 'range' not in self.refinement_kwargs: + raise ValueError(f"You must define 'range' in refinement_kwargs to use " + "{self.refinement_method} refinement.") + elif self.refinement_method == "gradient" and 'direction' not in self.refinement_kwargs: + raise ValueError("You must define 'direction' in refinement_kwargs to use " + "gradient refinement") + + @classmethod + def by_name(cls, input_spectrum, line_name, pixel, + refinement_method=None, refinement_kwargs={}): + # this does wavelength lookup and passes that wavelength to __init__ allowing the user + # to call CalibrationLine.by_name(spectrum, 'H_alpha', pixel=40) + # We could also have a type-check for ``wavelength`` above, this just feels more verbose + return cls(...) + + def _do_refinement(input_spectrum): + if self.refinement_method == 'gaussian': + window_width = self.refinement_kwargs.get('range') + window = SpectralRegion((self.pixel-window_width)*u.pix, + (self.pixel+window_width)*u.pix) + + # Use specutils.fit_lines to do model fitting. Define window in u.pix based on kwargs + input_model = Gaussian1D(mean=self.pixel, amplitude = self.input_spectrum[pixel], + stddev=3) + + fitted_model = fit_lines(self.input_spectrum, input_model, window=window) + new_pixel = fitted_model.mean + + elif self.refinement_method == 'min': + window_width = self.refinement_kwargs.get('range') + new_pixel = np.argmin(self.input_spectrum[self.pixel-window_width:self.pixel+window_width+1]) + + elif self.refinement_method = 'max': + window_width = self.refinement_kwargs.get('range') + new_pixel = np.argmax(self.input_spectrum[self.pixel-window_width:self.pixel+window_width+1]) + + return new_pixel + + def refine(self, input_spectrum=None, return_object=False): + # finds the center of the line according to refinement_method/kwargs and returns + # either the pixel value or a CalibrationLine object with pixel changed to the refined value + input_spectrum = self.input_spectrum if input_spectrum is None else input_spectrum + refined_pixel = _do_refinement(input_spectrum, self.refinement_method, **self.refinement_kwargs) + if return_object: + return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, ...) + + return refined_pixel + + @cached_property + def refined_pixel(self): + # returns the refined pixel for self.input_spectrum + return self.refine() + + @property + def with_refined_pixel(self): + # returns a copy of this object, but with the pixel updated to the refined value + return self.refine(return_object=True) + + +class WavelengthCalibration1D(): + + def __init__(input_spectrum, lines, model=Linear1D, + default_refinement_method=None, default_refinement_kwargs={}): + """ + input_spectrum: `~specutils.Spectrum1D` + A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. + lines: + List of lines to anchor the wavelength solution fit. coerced to CalibrationLine + objects if passed as tuples with default_refinement_method and default_refinement_kwargs as defaults + model: `~astropy.modeling.Model` + The model to fit for the wavelength solution. Defaults to a linear model. + default_refinement_method: str, optional + words + default_refinement_kwargs: dict, optional + words + """ + self.model = model + + + @classmethod + def autoidentify(cls, input_spectrum, line_list, model=Linear1D, ...): + # line_list could be a string ("common stellar") or an object + # this builds a list of CalibrationLine objects and passes to __init__ + return cls(...) + + @property + def refined_lines(self): + return [l.with_refined_pixel for l in self.lines] + + @property + def refined_pixels(self): + # useful for plotting over spectrum to change input lines/refinement options + return [(l.wavelength, l.refined_pixel) for l in self.lines] + + @cached_property + def wcs(self): + # computes and returns WCS after fitting self.model to self.refined_pixels + return WCS(...) + + def __call__(self, apply_to_spectrum=None): + # returns spectrum1d with wavelength calibration applied + # actual line refinement and WCS solution should already be done so that this can be called on multiple science sources + apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum + apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! + return apply_to_spectrum + + +class WavelengthCalibration2D(input_spectrum, trace, lines, model=Linear1D, + default_refinement_method=None, default_refinement_kwargs={}): + # input_spectrum must be 2-dimensional + # lines are coerced to CalibrationLine objects if passed as tuples with default_refinement_method and default_refinement_kwargs as defaults + @classmethod + def autoidentify(cls, input_spectrum, trace, line_list, model=Linear1D, ...): + # does this do 2D identification, or just search a 1d spectrum exactly at the trace? + # or should we require using WavelengthCalibration1D.autoidentify? + return cls(...) + + @property + def refined_lines(self): + return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=False)) for l in self.lines] for row in rows] + + @property + def refined_pixels(self): + return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=True)) for l in self.lines] for row in rows] + + @cached_property + def wcs(self): + # computes and returns WCS after fitting self.model to self.refined_pixels + return WCS(...) + + def __call__(self, apply_to_spectrum=None): + # returns spectrum1d with wavelength calibration applied + # actual line refinement and WCS solution should already be done so that this can be called on multiple science sources + apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum + apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! + return apply_to_spectrum + + + +# various supported syntaxes for creating a calibration object (do we want this much flexibility?): +cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine(5500, 20), CalibrationLine(6200, 32))) +cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine.by_name('H_alpha', 20, 'uphill'), CalibrationLine.by_name('Lyman_alpha', 32, 'max', {'range': 10}))) +cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), ('Lyman_alpha', 32))) +cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), (6000, 30)), default_refinement_method='gaussian') +cal1d = WavelengthCalibration1D.autoidentify(spec1d_lamp, common_stellar_line_list) + +# once we have that object, we can access the fitted wavelength solution via the WCS and apply it to +# a spectrum (either the same or separate from the spectrum used for fitting) +targ_wcs = cal1d.wcs +spec1d_targ1_cal = cal1d(spec1d_targ1) +spec1d_targ2_cal = cal1d(spec1d_targ2) + + +# we sould either start a 2D calibration the same way (except also passing the trace) or by passing +# the refined lines from the 1D calibration as the starting point +cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, ((5500, 20), (6000, 30))) +cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, cal1d.refined_lines) +spec2d_targ_cal = cal2d(spec2d_targ, trace_targ?) \ No newline at end of file From 1f3483aae86af337adae0f4da92ca5539c95ef48 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 19 Jan 2023 15:34:22 -0500 Subject: [PATCH 02/31] Testing/debugging CalibrationLine, starting on WavelengthCalibration --- specreduce/wavelength_calibration.py | 60 +++++++++++++++++----------- 1 file changed, 37 insertions(+), 23 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 09de76b7..55c466e7 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -1,5 +1,7 @@ from astropy.modeling.models import Linear1D, Gaussian1D import astropy.units as u +from functools import cached_property +import numpy as np from specutils import SpectralRegion from specutils.fitting import fit_lines @@ -10,7 +12,7 @@ class CalibrationLine(): - def __init__(input_spectrum, wavelength, pixel, refinement_method=None, + def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, refinement_kwargs={}): """ input_spectrum: `~specutils.Spectrum1D` @@ -53,26 +55,30 @@ def by_name(cls, input_spectrum, line_name, pixel, # We could also have a type-check for ``wavelength`` above, this just feels more verbose return cls(...) - def _do_refinement(input_spectrum): + def _do_refinement(self, input_spectrum): if self.refinement_method == 'gaussian': window_width = self.refinement_kwargs.get('range') window = SpectralRegion((self.pixel-window_width)*u.pix, (self.pixel+window_width)*u.pix) # Use specutils.fit_lines to do model fitting. Define window in u.pix based on kwargs - input_model = Gaussian1D(mean=self.pixel, amplitude = self.input_spectrum[pixel], - stddev=3) + input_model = Gaussian1D(mean=self.pixel, stddev=3, + amplitude = self.input_spectrum.flux[self.pixel]) fitted_model = fit_lines(self.input_spectrum, input_model, window=window) - new_pixel = fitted_model.mean + new_pixel = fitted_model.mean.value elif self.refinement_method == 'min': window_width = self.refinement_kwargs.get('range') - new_pixel = np.argmin(self.input_spectrum[self.pixel-window_width:self.pixel+window_width+1]) + new_pixel = np.argmin(self.input_spectrum.flux[self.pixel-window_width: + self.pixel+window_width+1]) + new_pixel += self.pixel - window_width - elif self.refinement_method = 'max': + elif self.refinement_method == 'max': window_width = self.refinement_kwargs.get('range') - new_pixel = np.argmax(self.input_spectrum[self.pixel-window_width:self.pixel+window_width+1]) + new_pixel = np.argmax(self.input_spectrum.flux[self.pixel-window_width: + self.pixel+window_width+1]) + new_pixel += self.pixel - window_width return new_pixel @@ -80,7 +86,7 @@ def refine(self, input_spectrum=None, return_object=False): # finds the center of the line according to refinement_method/kwargs and returns # either the pixel value or a CalibrationLine object with pixel changed to the refined value input_spectrum = self.input_spectrum if input_spectrum is None else input_spectrum - refined_pixel = _do_refinement(input_spectrum, self.refinement_method, **self.refinement_kwargs) + refined_pixel = self._do_refinement(input_spectrum) if return_object: return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, ...) @@ -118,12 +124,12 @@ def __init__(input_spectrum, lines, model=Linear1D, @classmethod - def autoidentify(cls, input_spectrum, line_list, model=Linear1D, ...): + def autoidentify(cls, input_spectrum, line_list, model=Linear1D): # line_list could be a string ("common stellar") or an object # this builds a list of CalibrationLine objects and passes to __init__ return cls(...) - @property + @property def refined_lines(self): return [l.with_refined_pixel for l in self.lines] @@ -135,22 +141,29 @@ def refined_pixels(self): @cached_property def wcs(self): # computes and returns WCS after fitting self.model to self.refined_pixels + inputs = np.array(self.refined_pixels) + return WCS(...) - def __call__(self, apply_to_spectrum=None): + def apply_to_spectrum(self, spectrum=None): # returns spectrum1d with wavelength calibration applied - # actual line refinement and WCS solution should already be done so that this can be called on multiple science sources - apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum - apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! - return apply_to_spectrum + # actual line refinement and WCS solution should already be done so that this can + # ibe called on multiple science sources + spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum + spectrum.wcs = self.wcs # might need a deepcopy! + return spectrum -class WavelengthCalibration2D(input_spectrum, trace, lines, model=Linear1D, - default_refinement_method=None, default_refinement_kwargs={}): +class WavelengthCalibration2D(): # input_spectrum must be 2-dimensional + # lines are coerced to CalibrationLine objects if passed as tuples with default_refinement_method and default_refinement_kwargs as defaults + def __init__(input_spectrum, trace, lines, model=Linear1D, + default_refinement_method=None, default_refinement_kwargs={}): + pass + @classmethod - def autoidentify(cls, input_spectrum, trace, line_list, model=Linear1D, ...): + def autoidentify(cls, input_spectrum, trace, line_list, model=Linear1D): # does this do 2D identification, or just search a 1d spectrum exactly at the trace? # or should we require using WavelengthCalibration1D.autoidentify? return cls(...) @@ -173,10 +186,10 @@ def __call__(self, apply_to_spectrum=None): # actual line refinement and WCS solution should already be done so that this can be called on multiple science sources apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! - return apply_to_spectrum - + return apply_to_spectrum +''' # various supported syntaxes for creating a calibration object (do we want this much flexibility?): cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine(5500, 20), CalibrationLine(6200, 32))) cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine.by_name('H_alpha', 20, 'uphill'), CalibrationLine.by_name('Lyman_alpha', 32, 'max', {'range': 10}))) @@ -184,7 +197,7 @@ def __call__(self, apply_to_spectrum=None): cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), (6000, 30)), default_refinement_method='gaussian') cal1d = WavelengthCalibration1D.autoidentify(spec1d_lamp, common_stellar_line_list) -# once we have that object, we can access the fitted wavelength solution via the WCS and apply it to +# once we have that object, we can access the fitted wavelength solution via the WCS and apply it to # a spectrum (either the same or separate from the spectrum used for fitting) targ_wcs = cal1d.wcs spec1d_targ1_cal = cal1d(spec1d_targ1) @@ -195,4 +208,5 @@ def __call__(self, apply_to_spectrum=None): # the refined lines from the 1D calibration as the starting point cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, ((5500, 20), (6000, 30))) cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, cal1d.refined_lines) -spec2d_targ_cal = cal2d(spec2d_targ, trace_targ?) \ No newline at end of file +spec2d_targ_cal = cal2d(spec2d_targ, trace_targ?) +''' From c54aed73c1aa79fbdc0e646c0aa9432f43c545c9 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 26 Jan 2023 14:33:29 -0500 Subject: [PATCH 03/31] Working on model fitting and construction of resulting GWCS --- specreduce/wavelength_calibration.py | 41 +++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 55c466e7..9730b4b8 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -1,6 +1,9 @@ from astropy.modeling.models import Linear1D, Gaussian1D +from astropy.modeling.fitting import LMLSQFitter import astropy.units as u from functools import cached_property +from gwcs import wcs +from gwcs import coordinate_frames as cf import numpy as np from specutils import SpectralRegion from specutils.fitting import fit_lines @@ -105,8 +108,9 @@ def with_refined_pixel(self): class WavelengthCalibration1D(): - def __init__(input_spectrum, lines, model=Linear1D, - default_refinement_method=None, default_refinement_kwargs={}): + def __init__(input_spectrum, lines, model=Linear1D, spectral_unit=u.Angstrom, + fitter=LMLSQFitter(calc_uncertainties=True), + default_refinement_method=None, default_refinement_kwargs={}): """ input_spectrum: `~specutils.Spectrum1D` A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. @@ -120,7 +124,24 @@ def __init__(input_spectrum, lines, model=Linear1D, default_refinement_kwargs: dict, optional words """ + self.input_spectrum = input_spectrum self.model = model + self.spectral_unit = spectral_unit + self.default_refinement_method = default_refinement_method + self.default_refinement_kwargs = default_refinement_kwargs + + if self.default_refinement_method in ("gaussian", "min", "max"): + if 'range' not in self.default_refinement_kwargs: + raise ValueError(f"You must define 'range' in default_refinement_kwargs to use " + "{self.refinement_method} refinement.") + elif (self.default_refinement_method == "gradient" and 'direction' not in + self.default_refinement_kwargs): + raise ValueError("You must define 'direction' in default_refinement_kwargs to use " + "gradient refinement") + + if isinstance(lines, str): + if lines in self.available_line_lists: + raise ValueError(f"Line list '{lines}' is not an available line list.") @classmethod @@ -141,9 +162,21 @@ def refined_pixels(self): @cached_property def wcs(self): # computes and returns WCS after fitting self.model to self.refined_pixels - inputs = np.array(self.refined_pixels) + x = np.array(self.refined_pixels)[:,0] + y = np.array(self.refined_pixels)[:,1] - return WCS(...) + # Fit the model + self.model = fitter(self.model, x, y) + + # Build a GWCS pipeline from the fitted model + pixel_frame = cf.SpectralFrame(axes_names=["x",], unit=[u.pix,]) + spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], unit=[self.spectral_unit,]) + + pipeline = [(pixel_frame, model),(spectral_frame, None)] + + wcsobj = wcs.WCS(pipeline) + + return wcsobj def apply_to_spectrum(self, spectrum=None): # returns spectrum1d with wavelength calibration applied From 5eb4aaabf6c689191e57b3c3407f70f266e7f9bb Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 26 Jan 2023 16:24:41 -0500 Subject: [PATCH 04/31] Working through bugs/testing --- specreduce/wavelength_calibration.py | 48 ++++++++++++---------------- 1 file changed, 20 insertions(+), 28 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 9730b4b8..891c625f 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -88,6 +88,8 @@ def _do_refinement(self, input_spectrum): def refine(self, input_spectrum=None, return_object=False): # finds the center of the line according to refinement_method/kwargs and returns # either the pixel value or a CalibrationLine object with pixel changed to the refined value + if self.refinement_method is None: + return self.pixel input_spectrum = self.input_spectrum if input_spectrum is None else input_spectrum refined_pixel = self._do_refinement(input_spectrum) if return_object: @@ -105,18 +107,22 @@ def with_refined_pixel(self): # returns a copy of this object, but with the pixel updated to the refined value return self.refine(return_object=True) + def __str__(self): + return(f"CalibrationLine: ({self.wavelength}, {self.pixel})") + class WavelengthCalibration1D(): - def __init__(input_spectrum, lines, model=Linear1D, spectral_unit=u.Angstrom, + def __init__(self, input_spectrum, lines, model=Linear1D, spectral_unit=u.Angstrom, fitter=LMLSQFitter(calc_uncertainties=True), default_refinement_method=None, default_refinement_kwargs={}): """ input_spectrum: `~specutils.Spectrum1D` A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. - lines: - List of lines to anchor the wavelength solution fit. coerced to CalibrationLine - objects if passed as tuples with default_refinement_method and default_refinement_kwargs as defaults + lines: str, list + List of lines to anchor the wavelength solution fit. List items are coerced to + CalibrationLine objects if passed as tuples of (pixel, wavelength) with + default_refinement_method and default_refinement_kwargs as defaults. model: `~astropy.modeling.Model` The model to fit for the wavelength solution. Defaults to a linear model. default_refinement_method: str, optional @@ -139,9 +145,18 @@ def __init__(input_spectrum, lines, model=Linear1D, spectral_unit=u.Angstrom, raise ValueError("You must define 'direction' in default_refinement_kwargs to use " "gradient refinement") + self.lines = [] if isinstance(lines, str): if lines in self.available_line_lists: raise ValueError(f"Line list '{lines}' is not an available line list.") + else: + for line in lines: + if isinstance(line, CalibrationLine): + self.lines.append(line) + else: + self.lines.append(CalibrationLine(self.input_spectrum, line[1], line[0], + self.default_refinement_method, + self.default_refinement_kwargs)) @classmethod @@ -169,7 +184,7 @@ def wcs(self): self.model = fitter(self.model, x, y) # Build a GWCS pipeline from the fitted model - pixel_frame = cf.SpectralFrame(axes_names=["x",], unit=[u.pix,]) + pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], unit=[self.spectral_unit,]) pipeline = [(pixel_frame, model),(spectral_frame, None)] @@ -220,26 +235,3 @@ def __call__(self, apply_to_spectrum=None): apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! return apply_to_spectrum - - -''' -# various supported syntaxes for creating a calibration object (do we want this much flexibility?): -cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine(5500, 20), CalibrationLine(6200, 32))) -cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine.by_name('H_alpha', 20, 'uphill'), CalibrationLine.by_name('Lyman_alpha', 32, 'max', {'range': 10}))) -cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), ('Lyman_alpha', 32))) -cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), (6000, 30)), default_refinement_method='gaussian') -cal1d = WavelengthCalibration1D.autoidentify(spec1d_lamp, common_stellar_line_list) - -# once we have that object, we can access the fitted wavelength solution via the WCS and apply it to -# a spectrum (either the same or separate from the spectrum used for fitting) -targ_wcs = cal1d.wcs -spec1d_targ1_cal = cal1d(spec1d_targ1) -spec1d_targ2_cal = cal1d(spec1d_targ2) - - -# we sould either start a 2D calibration the same way (except also passing the trace) or by passing -# the refined lines from the 1D calibration as the starting point -cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, ((5500, 20), (6000, 30))) -cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, cal1d.refined_lines) -spec2d_targ_cal = cal2d(spec2d_targ, trace_targ?) -''' From baa706c8c69fd61db26cc338c3fa915d7ffc8507 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 26 Jan 2023 17:22:25 -0500 Subject: [PATCH 05/31] Working initial implementation --- specreduce/wavelength_calibration.py | 29 +++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 891c625f..e447e2aa 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -5,7 +5,7 @@ from gwcs import wcs from gwcs import coordinate_frames as cf import numpy as np -from specutils import SpectralRegion +from specutils import SpectralRegion, Spectrum1D from specutils.fitting import fit_lines @@ -93,7 +93,9 @@ def refine(self, input_spectrum=None, return_object=False): input_spectrum = self.input_spectrum if input_spectrum is None else input_spectrum refined_pixel = self._do_refinement(input_spectrum) if return_object: - return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, ...) + return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, + refinement_method=self.refinement_method, + refinement_kwargs=self.refinement_kwargs) return refined_pixel @@ -110,10 +112,13 @@ def with_refined_pixel(self): def __str__(self): return(f"CalibrationLine: ({self.wavelength}, {self.pixel})") + def __repr__(self): + return(f"CalibrationLine({self.wavelength}, {self.pixel})") + class WavelengthCalibration1D(): - def __init__(self, input_spectrum, lines, model=Linear1D, spectral_unit=u.Angstrom, + def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angstrom, fitter=LMLSQFitter(calc_uncertainties=True), default_refinement_method=None, default_refinement_kwargs={}): """ @@ -133,6 +138,7 @@ def __init__(self, input_spectrum, lines, model=Linear1D, spectral_unit=u.Angstr self.input_spectrum = input_spectrum self.model = model self.spectral_unit = spectral_unit + self.fitter = fitter self.default_refinement_method = default_refinement_method self.default_refinement_kwargs = default_refinement_kwargs @@ -177,17 +183,17 @@ def refined_pixels(self): @cached_property def wcs(self): # computes and returns WCS after fitting self.model to self.refined_pixels - x = np.array(self.refined_pixels)[:,0] - y = np.array(self.refined_pixels)[:,1] + x = [line[1] for line in self.refined_pixels] * u.pix + y = [line[0].value for line in self.refined_pixels] * self.spectral_unit # Fit the model - self.model = fitter(self.model, x, y) + self.model = self.fitter(self.model, x, y) # Build a GWCS pipeline from the fitted model pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], unit=[self.spectral_unit,]) - pipeline = [(pixel_frame, model),(spectral_frame, None)] + pipeline = [(pixel_frame, self.model),(spectral_frame, None)] wcsobj = wcs.WCS(pipeline) @@ -196,10 +202,11 @@ def wcs(self): def apply_to_spectrum(self, spectrum=None): # returns spectrum1d with wavelength calibration applied # actual line refinement and WCS solution should already be done so that this can - # ibe called on multiple science sources - spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum - spectrum.wcs = self.wcs # might need a deepcopy! - return spectrum + # be called on multiple science sources + spectrum = self.input_spectrum if spectrum is None else spectrum + updated_spectrum = Spectrum1D(spectrum.flux, wcs=self.wcs, mask=spectrum.mask, + uncertainty=spectrum.uncertainty) + return updated_spectrum class WavelengthCalibration2D(): From 02ed67e4c4abc8dfda595ffce15f895cd3a3ea6e Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 26 Jan 2023 17:34:35 -0500 Subject: [PATCH 06/31] Codestyle --- specreduce/wavelength_calibration.py | 111 ++++++++++++++------------- 1 file changed, 56 insertions(+), 55 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index e447e2aa..635a244f 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -9,10 +9,6 @@ from specutils.fitting import fit_lines -""" -Starting from Kyle's proposed pseudocode from https://github.com/astropy/specreduce/pull/152 -""" - class CalibrationLine(): def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, @@ -45,7 +41,7 @@ def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, if self.refinement_method in ("gaussian", "min", "max"): if 'range' not in self.refinement_kwargs: raise ValueError(f"You must define 'range' in refinement_kwargs to use " - "{self.refinement_method} refinement.") + f"{self.refinement_method} refinement.") elif self.refinement_method == "gradient" and 'direction' not in self.refinement_kwargs: raise ValueError("You must define 'direction' in refinement_kwargs to use " "gradient refinement") @@ -53,37 +49,36 @@ def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, @classmethod def by_name(cls, input_spectrum, line_name, pixel, refinement_method=None, refinement_kwargs={}): - # this does wavelength lookup and passes that wavelength to __init__ allowing the user - # to call CalibrationLine.by_name(spectrum, 'H_alpha', pixel=40) - # We could also have a type-check for ``wavelength`` above, this just feels more verbose - return cls(...) + # this will do wavelength lookup and passes that wavelength to __init__ allowing the user + # to call CalibrationLine.by_name(spectrum, 'H_alpha', pixel=40) + # We could also have a type-check for ``wavelength`` above, this just feels more verbose + return cls(...) def _do_refinement(self, input_spectrum): - if self.refinement_method == 'gaussian': - window_width = self.refinement_kwargs.get('range') - window = SpectralRegion((self.pixel-window_width)*u.pix, - (self.pixel+window_width)*u.pix) + if self.refinement_method == 'gaussian': + window_width = self.refinement_kwargs.get('range') + window = SpectralRegion((self.pixel-window_width)*u.pix, + (self.pixel+window_width)*u.pix) - # Use specutils.fit_lines to do model fitting. Define window in u.pix based on kwargs - input_model = Gaussian1D(mean=self.pixel, stddev=3, - amplitude = self.input_spectrum.flux[self.pixel]) + input_model = Gaussian1D(mean=self.pixel, stddev=3, + amplitude=self.input_spectrum.flux[self.pixel]) - fitted_model = fit_lines(self.input_spectrum, input_model, window=window) - new_pixel = fitted_model.mean.value + fitted_model = fit_lines(self.input_spectrum, input_model, window=window) + new_pixel = fitted_model.mean.value - elif self.refinement_method == 'min': - window_width = self.refinement_kwargs.get('range') - new_pixel = np.argmin(self.input_spectrum.flux[self.pixel-window_width: - self.pixel+window_width+1]) - new_pixel += self.pixel - window_width + elif self.refinement_method == 'min': + window_width = self.refinement_kwargs.get('range') + new_pixel = np.argmin(self.input_spectrum.flux[self.pixel-window_width: + self.pixel+window_width+1]) + new_pixel += self.pixel - window_width - elif self.refinement_method == 'max': - window_width = self.refinement_kwargs.get('range') - new_pixel = np.argmax(self.input_spectrum.flux[self.pixel-window_width: - self.pixel+window_width+1]) - new_pixel += self.pixel - window_width + elif self.refinement_method == 'max': + window_width = self.refinement_kwargs.get('range') + new_pixel = np.argmax(self.input_spectrum.flux[self.pixel-window_width: + self.pixel+window_width+1]) + new_pixel += self.pixel - window_width - return new_pixel + return new_pixel def refine(self, input_spectrum=None, return_object=False): # finds the center of the line according to refinement_method/kwargs and returns @@ -93,7 +88,7 @@ def refine(self, input_spectrum=None, return_object=False): input_spectrum = self.input_spectrum if input_spectrum is None else input_spectrum refined_pixel = self._do_refinement(input_spectrum) if return_object: - return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, + return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, refinement_method=self.refinement_method, refinement_kwargs=self.refinement_kwargs) @@ -110,10 +105,10 @@ def with_refined_pixel(self): return self.refine(return_object=True) def __str__(self): - return(f"CalibrationLine: ({self.wavelength}, {self.pixel})") + return f"CalibrationLine: ({self.wavelength}, {self.pixel})" def __repr__(self): - return(f"CalibrationLine({self.wavelength}, {self.pixel})") + return f"CalibrationLine({self.wavelength}, {self.pixel})" class WavelengthCalibration1D(): @@ -144,8 +139,8 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs if self.default_refinement_method in ("gaussian", "min", "max"): if 'range' not in self.default_refinement_kwargs: - raise ValueError(f"You must define 'range' in default_refinement_kwargs to use " - "{self.refinement_method} refinement.") + raise ValueError("You must define 'range' in default_refinement_kwargs to use " + f"{self.refinement_method} refinement.") elif (self.default_refinement_method == "gradient" and 'direction' not in self.default_refinement_kwargs): raise ValueError("You must define 'direction' in default_refinement_kwargs to use " @@ -164,7 +159,6 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs self.default_refinement_method, self.default_refinement_kwargs)) - @classmethod def autoidentify(cls, input_spectrum, line_list, model=Linear1D): # line_list could be a string ("common stellar") or an object @@ -173,12 +167,12 @@ def autoidentify(cls, input_spectrum, line_list, model=Linear1D): @property def refined_lines(self): - return [l.with_refined_pixel for l in self.lines] + return [line.with_refined_pixel for line in self.lines] @property def refined_pixels(self): # useful for plotting over spectrum to change input lines/refinement options - return [(l.wavelength, l.refined_pixel) for l in self.lines] + return [(line.wavelength, line.refined_pixel) for line in self.lines] @cached_property def wcs(self): @@ -193,29 +187,32 @@ def wcs(self): pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], unit=[self.spectral_unit,]) - pipeline = [(pixel_frame, self.model),(spectral_frame, None)] + pipeline = [(pixel_frame, self.model), (spectral_frame, None)] wcsobj = wcs.WCS(pipeline) return wcsobj def apply_to_spectrum(self, spectrum=None): - # returns spectrum1d with wavelength calibration applied - # actual line refinement and WCS solution should already be done so that this can - # be called on multiple science sources - spectrum = self.input_spectrum if spectrum is None else spectrum - updated_spectrum = Spectrum1D(spectrum.flux, wcs=self.wcs, mask=spectrum.mask, - uncertainty=spectrum.uncertainty) - return updated_spectrum + # returns spectrum1d with wavelength calibration applied + # actual line refinement and WCS solution should already be done so that this can + # be called on multiple science sources + spectrum = self.input_spectrum if spectrum is None else spectrum + updated_spectrum = Spectrum1D(spectrum.flux, wcs=self.wcs, mask=spectrum.mask, + uncertainty=spectrum.uncertainty) + return updated_spectrum +''' +# WavelengthCalibration2D is a planned future feature class WavelengthCalibration2D(): # input_spectrum must be 2-dimensional - # lines are coerced to CalibrationLine objects if passed as tuples with default_refinement_method and default_refinement_kwargs as defaults + # lines are coerced to CalibrationLine objects if passed as tuples with + # default_refinement_method and default_refinement_kwargs as defaults def __init__(input_spectrum, trace, lines, model=Linear1D, default_refinement_method=None, default_refinement_kwargs={}): - pass + return NotImplementedError("2D wavelength calibration is not yet implemented") @classmethod def autoidentify(cls, input_spectrum, trace, line_list, model=Linear1D): @@ -225,20 +222,24 @@ def autoidentify(cls, input_spectrum, trace, line_list, model=Linear1D): @property def refined_lines(self): - return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=False)) for l in self.lines] for row in rows] + return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=False)) + for l in self.lines] for row in rows] @property def refined_pixels(self): - return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=True)) for l in self.lines] for row in rows] + return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=True)) + for l in self.lines] for row in rows] @cached_property def wcs(self): # computes and returns WCS after fitting self.model to self.refined_pixels - return WCS(...) + pass def __call__(self, apply_to_spectrum=None): - # returns spectrum1d with wavelength calibration applied - # actual line refinement and WCS solution should already be done so that this can be called on multiple science sources - apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum - apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! - return apply_to_spectrum + # returns spectrum1d with wavelength calibration applied + # actual line refinement and WCS solution should already be done so that this + # can be called on multiple science sources + apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum + apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! + return apply_to_spectrum +''' From 1df8621016d0bd0c260b3bea122446405c2cc062 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 27 Jan 2023 13:26:27 -0500 Subject: [PATCH 07/31] Match order of pixel/wavelength between CalibrationLine and WavelengthCalibration1D, add tests --- .../tests/test_wavelength_calibration.py | 36 +++++++++++++++++++ specreduce/wavelength_calibration.py | 4 +-- 2 files changed, 38 insertions(+), 2 deletions(-) create mode 100644 specreduce/tests/test_wavelength_calibration.py diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py new file mode 100644 index 00000000..90082df0 --- /dev/null +++ b/specreduce/tests/test_wavelength_calibration.py @@ -0,0 +1,36 @@ +import numpy as np +import pytest + +import astropy.units as u +from astropy.tests.helper import assert_quantity_allclose +from astropy.utils.exceptions import AstropyUserWarning +from specutils import Spectrum1D + +from specreduce.wavelength_calibration import CalibrationLine, WavelengthCalibration1D + +def test_linear_from_list(): + np.random.seed(7) + flux = np.random.random(50)*u.Jy + sa = np.arange(0,50)*u.pix + spec = Spectrum1D(flux, spectral_axis = sa) + test = WavelengthCalibration1D(spec, [(5000*u.AA, 0),(5100*u.AA, 10),(5198*u.AA, 20),(5305*u.AA, 30)]) + with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): + spec2 = test.apply_to_spectrum(spec) + + assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) + assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) + + +def test_linear_from_calibrationline(): + np.random.seed(7) + flux = np.random.random(50)*u.Jy + sa = np.arange(0,50)*u.pix + spec = Spectrum1D(flux, spectral_axis = sa) + lines = [CalibrationLine(spec, 5000*u.AA, 0), CalibrationLine(spec, 5100*u.AA, 10), + CalibrationLine(spec, 5198*u.AA, 20), CalibrationLine(spec, 5305*u.AA, 30)] + test = WavelengthCalibration1D(spec, lines) + with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): + spec2 = test.apply_to_spectrum(spec) + + assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) + assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) \ No newline at end of file diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 635a244f..a0bc82cb 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -121,7 +121,7 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. lines: str, list List of lines to anchor the wavelength solution fit. List items are coerced to - CalibrationLine objects if passed as tuples of (pixel, wavelength) with + CalibrationLine objects if passed as tuples of (wavelength, pixel) with default_refinement_method and default_refinement_kwargs as defaults. model: `~astropy.modeling.Model` The model to fit for the wavelength solution. Defaults to a linear model. @@ -155,7 +155,7 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs if isinstance(line, CalibrationLine): self.lines.append(line) else: - self.lines.append(CalibrationLine(self.input_spectrum, line[1], line[0], + self.lines.append(CalibrationLine(self.input_spectrum, line[0], line[1], self.default_refinement_method, self.default_refinement_kwargs)) From eb5feeac3072b09ce3daf17129e128a2aa68a243 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 27 Jan 2023 13:29:18 -0500 Subject: [PATCH 08/31] Fix codestyle in test --- .../tests/test_wavelength_calibration.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 90082df0..f6fc1a1e 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -8,15 +8,17 @@ from specreduce.wavelength_calibration import CalibrationLine, WavelengthCalibration1D + def test_linear_from_list(): np.random.seed(7) flux = np.random.random(50)*u.Jy - sa = np.arange(0,50)*u.pix - spec = Spectrum1D(flux, spectral_axis = sa) - test = WavelengthCalibration1D(spec, [(5000*u.AA, 0),(5100*u.AA, 10),(5198*u.AA, 20),(5305*u.AA, 30)]) + sa = np.arange(0, 50)*u.pix + spec = Spectrum1D(flux, spectral_axis=sa) + test = WavelengthCalibration1D(spec, [(5000*u.AA, 0), (5100*u.AA, 10), + (5198*u.AA, 20), (5305*u.AA, 30)]) with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): spec2 = test.apply_to_spectrum(spec) - + assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) @@ -24,13 +26,13 @@ def test_linear_from_list(): def test_linear_from_calibrationline(): np.random.seed(7) flux = np.random.random(50)*u.Jy - sa = np.arange(0,50)*u.pix - spec = Spectrum1D(flux, spectral_axis = sa) + sa = np.arange(0, 50)*u.pix + spec = Spectrum1D(flux, spectral_axis=sa) lines = [CalibrationLine(spec, 5000*u.AA, 0), CalibrationLine(spec, 5100*u.AA, 10), CalibrationLine(spec, 5198*u.AA, 20), CalibrationLine(spec, 5305*u.AA, 30)] test = WavelengthCalibration1D(spec, lines) with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): spec2 = test.apply_to_spectrum(spec) - + assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) - assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) \ No newline at end of file + assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) From 6ea98e491ee0c13be93a59198bb0941b0d9b42f8 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Tue, 31 Jan 2023 11:37:47 -0500 Subject: [PATCH 09/31] Adding tests for CalibrationLine --- specreduce/conftest.py | 31 ++++++++++++ .../tests/test_wavelength_calibration.py | 47 ++++++++++++------- specreduce/wavelength_calibration.py | 3 ++ 3 files changed, 64 insertions(+), 17 deletions(-) diff --git a/specreduce/conftest.py b/specreduce/conftest.py index 672b2733..d5093e1b 100644 --- a/specreduce/conftest.py +++ b/specreduce/conftest.py @@ -6,6 +6,10 @@ import os from astropy.version import version as astropy_version +import astropy.units as u +import numpy as np +import pytest +from specutils import Spectrum1D # For Astropy 3.0 and later, we can use the standalone pytest plugin if astropy_version < '3.0': @@ -19,6 +23,33 @@ except ImportError: ASTROPY_HEADER = False +@pytest.fixture +def spec1d(): + np.random.seed(7) + flux = np.random.random(50)*u.Jy + sa = np.arange(0, 50)*u.pix + spec = Spectrum1D(flux, spectral_axis=sa) + return spec + +@pytest.fixture +def spec1d_with_emission_line(): + np.random.seed(7) + sa = np.arange(0, 200)*u.pix + flux = (np.random.randn(200) + + 10*np.exp(-0.01*((sa.value-130)**2)) + + sa.value/100) * u.Jy + spec = Spectrum1D(flux, spectral_axis=sa) + return spec + +@pytest.fixture +def spec1d_with_absorption_line(): + np.random.seed(7) + sa = np.arange(0, 200)*u.pix + flux = (np.random.randn(200) - + 10*np.exp(-0.01*((sa.value-130)**2)) + + sa.value/100) * u.Jy + spec = Spectrum1D(flux, spectral_axis=sa) + return spec def pytest_configure(config): diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index f6fc1a1e..dfb1d8d2 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -1,4 +1,5 @@ import numpy as np +from numpy.testing import assert_allclose import pytest import astropy.units as u @@ -8,31 +9,43 @@ from specreduce.wavelength_calibration import CalibrationLine, WavelengthCalibration1D - -def test_linear_from_list(): - np.random.seed(7) - flux = np.random.random(50)*u.Jy - sa = np.arange(0, 50)*u.pix - spec = Spectrum1D(flux, spectral_axis=sa) - test = WavelengthCalibration1D(spec, [(5000*u.AA, 0), (5100*u.AA, 10), +def test_linear_from_list(spec1d): + test = WavelengthCalibration1D(spec1d, [(5000*u.AA, 0), (5100*u.AA, 10), (5198*u.AA, 20), (5305*u.AA, 30)]) with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): - spec2 = test.apply_to_spectrum(spec) + spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) -def test_linear_from_calibrationline(): - np.random.seed(7) - flux = np.random.random(50)*u.Jy - sa = np.arange(0, 50)*u.pix - spec = Spectrum1D(flux, spectral_axis=sa) - lines = [CalibrationLine(spec, 5000*u.AA, 0), CalibrationLine(spec, 5100*u.AA, 10), - CalibrationLine(spec, 5198*u.AA, 20), CalibrationLine(spec, 5305*u.AA, 30)] - test = WavelengthCalibration1D(spec, lines) +def test_linear_from_calibrationline(spec1d): + lines = [CalibrationLine(spec1d, 5000*u.AA, 0), CalibrationLine(spec1d, 5100*u.AA, 10), + CalibrationLine(spec1d, 5198*u.AA, 20), CalibrationLine(spec1d, 5305*u.AA, 30)] + test = WavelengthCalibration1D(spec1d, lines) with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): - spec2 = test.apply_to_spectrum(spec) + spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) + +def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line): + with pytest.raises(ValueError, match="You must define 'range' in refinement_kwargs"): + line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, + refinement_method='gaussian') + + with pytest.raises(ValueError, match="You must define 'direction' in refinement_kwargs"): + line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, + refinement_method='gradient') + + line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, refinement_method='gaussian', + refinement_kwargs={'range': 25}) + assert_allclose(line.refine(), 129.44371) + + line2 = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 130, refinement_method='max', + refinement_kwargs={'range': 10}) + assert line2.refine() == 128 + + line3 = CalibrationLine(spec1d_with_absorption_line, 5000*u.AA, 128, refinement_method='min', + refinement_kwargs={'range': 10}) + assert line3.refine() == 130 \ No newline at end of file diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index a0bc82cb..90b8fd46 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -78,6 +78,9 @@ def _do_refinement(self, input_spectrum): self.pixel+window_width+1]) new_pixel += self.pixel - window_width + else: + raise ValueError(f"Refinement method {self.refinement_method} is not implemented") + return new_pixel def refine(self, input_spectrum=None, return_object=False): From af4a5230fbbf13aee7db023897752b422985ea3f Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Tue, 31 Jan 2023 12:00:09 -0500 Subject: [PATCH 10/31] Fix codestlye --- specreduce/conftest.py | 12 ++++++++---- specreduce/tests/test_wavelength_calibration.py | 8 ++++---- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/specreduce/conftest.py b/specreduce/conftest.py index d5093e1b..87ee89f1 100644 --- a/specreduce/conftest.py +++ b/specreduce/conftest.py @@ -23,6 +23,7 @@ except ImportError: ASTROPY_HEADER = False + @pytest.fixture def spec1d(): np.random.seed(7) @@ -31,26 +32,29 @@ def spec1d(): spec = Spectrum1D(flux, spectral_axis=sa) return spec + @pytest.fixture def spec1d_with_emission_line(): np.random.seed(7) sa = np.arange(0, 200)*u.pix flux = (np.random.randn(200) + - 10*np.exp(-0.01*((sa.value-130)**2)) + - sa.value/100) * u.Jy + 10*np.exp(-0.01*((sa.value-130)**2)) + + sa.value/100) * u.Jy spec = Spectrum1D(flux, spectral_axis=sa) return spec + @pytest.fixture def spec1d_with_absorption_line(): np.random.seed(7) sa = np.arange(0, 200)*u.pix flux = (np.random.randn(200) - - 10*np.exp(-0.01*((sa.value-130)**2)) + - sa.value/100) * u.Jy + 10*np.exp(-0.01*((sa.value-130)**2)) + + sa.value/100) * u.Jy spec = Spectrum1D(flux, spectral_axis=sa) return spec + def pytest_configure(config): if ASTROPY_HEADER: diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index dfb1d8d2..e0420e6b 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -1,17 +1,16 @@ -import numpy as np from numpy.testing import assert_allclose import pytest import astropy.units as u from astropy.tests.helper import assert_quantity_allclose from astropy.utils.exceptions import AstropyUserWarning -from specutils import Spectrum1D from specreduce.wavelength_calibration import CalibrationLine, WavelengthCalibration1D + def test_linear_from_list(spec1d): test = WavelengthCalibration1D(spec1d, [(5000*u.AA, 0), (5100*u.AA, 10), - (5198*u.AA, 20), (5305*u.AA, 30)]) + (5198*u.AA, 20), (5305*u.AA, 30)]) with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): spec2 = test.apply_to_spectrum(spec1d) @@ -29,6 +28,7 @@ def test_linear_from_calibrationline(spec1d): assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) + def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line): with pytest.raises(ValueError, match="You must define 'range' in refinement_kwargs"): line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, @@ -48,4 +48,4 @@ def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line) line3 = CalibrationLine(spec1d_with_absorption_line, 5000*u.AA, 128, refinement_method='min', refinement_kwargs={'range': 10}) - assert line3.refine() == 130 \ No newline at end of file + assert line3.refine() == 130 From 568fd0911d955d4292f69c12eaa7f2d659d18243 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Tue, 31 Jan 2023 12:02:54 -0500 Subject: [PATCH 11/31] Increase atol for gaussian center --- specreduce/tests/test_wavelength_calibration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index e0420e6b..d7ccb6f3 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -40,7 +40,7 @@ def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line) line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, refinement_method='gaussian', refinement_kwargs={'range': 25}) - assert_allclose(line.refine(), 129.44371) + assert_allclose(line.refine(), 129.44371, atol=0.01) line2 = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 130, refinement_method='max', refinement_kwargs={'range': 10}) From be35c0919b3d8f31ab209f773a0146fd01354307 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 1 Feb 2023 11:21:57 -0500 Subject: [PATCH 12/31] Remove autoidentify for now, add Polynomial1D test --- specreduce/tests/test_wavelength_calibration.py | 12 ++++++++++++ specreduce/wavelength_calibration.py | 8 +------- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index d7ccb6f3..9723d134 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -2,6 +2,7 @@ import pytest import astropy.units as u +from astropy.modeling.models import Polynomial1D from astropy.tests.helper import assert_quantity_allclose from astropy.utils.exceptions import AstropyUserWarning @@ -28,6 +29,17 @@ def test_linear_from_calibrationline(spec1d): assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) +def test_poly_from_calibrationline(spec1d): + # This test is mostly to prove that you can use other models + lines = [CalibrationLine(spec1d, 5005*u.AA, 0), CalibrationLine(spec1d, 5110*u.AA, 10), + CalibrationLine(spec1d, 5214*u.AA, 20), CalibrationLine(spec1d, 5330*u.AA, 30), + CalibrationLine(spec1d, 5438*u.AA, 40)] + test = WavelengthCalibration1D(spec1d, lines, model=Polynomial1D(2)) + with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): + spec2 = test.apply_to_spectrum(spec1d) + + assert_allclose(test.model.parameters, [5.00477143e+03, 1.03457143e+01, 1.28571429e-02]) + def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line): with pytest.raises(ValueError, match="You must define 'range' in refinement_kwargs"): diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 90b8fd46..df0ba534 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -102,7 +102,7 @@ def refined_pixel(self): # returns the refined pixel for self.input_spectrum return self.refine() - @property + @cached_property def with_refined_pixel(self): # returns a copy of this object, but with the pixel updated to the refined value return self.refine(return_object=True) @@ -162,12 +162,6 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs self.default_refinement_method, self.default_refinement_kwargs)) - @classmethod - def autoidentify(cls, input_spectrum, line_list, model=Linear1D): - # line_list could be a string ("common stellar") or an object - # this builds a list of CalibrationLine objects and passes to __init__ - return cls(...) - @property def refined_lines(self): return [line.with_refined_pixel for line in self.lines] From f3b1b9fda616a505cef8af20cedc27dcf368c7ee Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 1 Feb 2023 12:07:29 -0500 Subject: [PATCH 13/31] Add cache clearing, avoid AstropyUserWarning --- .../tests/test_wavelength_calibration.py | 13 +++---- specreduce/wavelength_calibration.py | 37 +++++++++++++++++-- 2 files changed, 38 insertions(+), 12 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 9723d134..ab209ab4 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -4,7 +4,6 @@ import astropy.units as u from astropy.modeling.models import Polynomial1D from astropy.tests.helper import assert_quantity_allclose -from astropy.utils.exceptions import AstropyUserWarning from specreduce.wavelength_calibration import CalibrationLine, WavelengthCalibration1D @@ -12,8 +11,7 @@ def test_linear_from_list(spec1d): test = WavelengthCalibration1D(spec1d, [(5000*u.AA, 0), (5100*u.AA, 10), (5198*u.AA, 20), (5305*u.AA, 30)]) - with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): - spec2 = test.apply_to_spectrum(spec1d) + spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) @@ -23,21 +21,20 @@ def test_linear_from_calibrationline(spec1d): lines = [CalibrationLine(spec1d, 5000*u.AA, 0), CalibrationLine(spec1d, 5100*u.AA, 10), CalibrationLine(spec1d, 5198*u.AA, 20), CalibrationLine(spec1d, 5305*u.AA, 30)] test = WavelengthCalibration1D(spec1d, lines) - with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): - spec2 = test.apply_to_spectrum(spec1d) + spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) + def test_poly_from_calibrationline(spec1d): # This test is mostly to prove that you can use other models lines = [CalibrationLine(spec1d, 5005*u.AA, 0), CalibrationLine(spec1d, 5110*u.AA, 10), CalibrationLine(spec1d, 5214*u.AA, 20), CalibrationLine(spec1d, 5330*u.AA, 30), CalibrationLine(spec1d, 5438*u.AA, 40)] test = WavelengthCalibration1D(spec1d, lines, model=Polynomial1D(2)) - with pytest.warns(AstropyUserWarning, match="Model is linear in parameters"): - spec2 = test.apply_to_spectrum(spec1d) - + test.apply_to_spectrum(spec1d) + assert_allclose(test.model.parameters, [5.00477143e+03, 1.03457143e+01, 1.28571429e-02]) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index df0ba534..4f370743 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -1,5 +1,5 @@ from astropy.modeling.models import Linear1D, Gaussian1D -from astropy.modeling.fitting import LMLSQFitter +from astropy.modeling.fitting import LMLSQFitter, LinearLSQFitter import astropy.units as u from functools import cached_property from gwcs import wcs @@ -37,6 +37,7 @@ def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, self.pixel = pixel self.refinement_method = refinement_method self.refinement_kwargs = refinement_kwargs + self._cached_properties = ['refined_pixel', 'with_refined_pixel'] if self.refinement_method in ("gaussian", "min", "max"): if 'range' not in self.refinement_kwargs: @@ -46,6 +47,16 @@ def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, raise ValueError("You must define 'direction' in refinement_kwargs to use " "gradient refinement") + def _clear_cache(self, *attrs): + """ + provide convenience function to clearing the cache for cached_properties + """ + if not len(attrs): + attrs = self._cached_properties + for attr in attrs: + if attr in self.__dict__: + del self.__dict__[attr] + @classmethod def by_name(cls, input_spectrum, line_name, pixel, refinement_method=None, refinement_kwargs={}): @@ -117,8 +128,7 @@ def __repr__(self): class WavelengthCalibration1D(): def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angstrom, - fitter=LMLSQFitter(calc_uncertainties=True), - default_refinement_method=None, default_refinement_kwargs={}): + fitter=None, default_refinement_method=None, default_refinement_kwargs={}): """ input_spectrum: `~specutils.Spectrum1D` A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. @@ -136,9 +146,18 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs self.input_spectrum = input_spectrum self.model = model self.spectral_unit = spectral_unit - self.fitter = fitter self.default_refinement_method = default_refinement_method self.default_refinement_kwargs = default_refinement_kwargs + self._cached_properties = ['wcs',] + + if fitter is None: + print("Got here") + if self.model.linear: + self.fitter = LinearLSQFitter(calc_uncertainties=True) + else: + self.fitter = LMLSQFitter(calc_uncertainties=True) + else: + self.fitter = fitter if self.default_refinement_method in ("gaussian", "min", "max"): if 'range' not in self.default_refinement_kwargs: @@ -162,6 +181,16 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs self.default_refinement_method, self.default_refinement_kwargs)) + def _clear_cache(self, *attrs): + """ + provide convenience function to clearing the cache for cached_properties + """ + if not len(attrs): + attrs = self._cached_properties + for attr in attrs: + if attr in self.__dict__: + del self.__dict__[attr] + @property def refined_lines(self): return [line.with_refined_pixel for line in self.lines] From 8117864fc38c3271363da34cded013ac188638f6 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 1 Feb 2023 13:32:47 -0500 Subject: [PATCH 14/31] Add input_spectrum property and clear cache on setting it --- .../tests/test_wavelength_calibration.py | 17 ++++++++++++ specreduce/wavelength_calibration.py | 27 ++++++++++++++++--- 2 files changed, 41 insertions(+), 3 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index ab209ab4..75e99c2a 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -58,3 +58,20 @@ def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line) line3 = CalibrationLine(spec1d_with_absorption_line, 5000*u.AA, 128, refinement_method='min', refinement_kwargs={'range': 10}) assert line3.refine() == 130 + + +def test_replace_spectrum(spec1d, spec1d_with_emission_line): + lines = [CalibrationLine(spec1d, 5000*u.AA, 0), CalibrationLine(spec1d, 5100*u.AA, 10), + CalibrationLine(spec1d, 5198*u.AA, 20), CalibrationLine(spec1d, 5305*u.AA, 30)] + test = WavelengthCalibration1D(spec1d, lines) + # Accessing this property causes fits the model and caches the resulting WCS + test.wcs + assert "wcs" in test.__dict__ + for line in test.lines: + assert "refined_pixel" in line.__dict__ + + # Replace the input spectrum, which should clear the cached properties + test.input_spectrum = spec1d_with_emission_line + assert "wcs" not in test.__dict__ + for line in test.lines: + assert "refined_pixel" not in line.__dict__ diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 4f370743..83f69367 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -32,7 +32,7 @@ def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, of line center to include on either side of gaussian fit is set by int ``range``. Gradient descent/ascent is determined by setting ``direction`` to 'min' or 'max'. """ - self.input_spectrum = input_spectrum + self._input_spectrum = input_spectrum self.wavelength = wavelength self.pixel = pixel self.refinement_method = refinement_method @@ -118,6 +118,16 @@ def with_refined_pixel(self): # returns a copy of this object, but with the pixel updated to the refined value return self.refine(return_object=True) + @property + def input_spectrum(self): + return self._input_spectrum + + @input_spectrum.setter + def input_spectrum(self, new_spectrum): + # We want to clear the refined locations if a new calibration spectrum is provided + self._clear_cache() + self._input_spectrum = new_spectrum + def __str__(self): return f"CalibrationLine: ({self.wavelength}, {self.pixel})" @@ -143,7 +153,7 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs default_refinement_kwargs: dict, optional words """ - self.input_spectrum = input_spectrum + self._input_spectrum = input_spectrum self.model = model self.spectral_unit = spectral_unit self.default_refinement_method = default_refinement_method @@ -151,7 +161,6 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs self._cached_properties = ['wcs',] if fitter is None: - print("Got here") if self.model.linear: self.fitter = LinearLSQFitter(calc_uncertainties=True) else: @@ -191,6 +200,18 @@ def _clear_cache(self, *attrs): if attr in self.__dict__: del self.__dict__[attr] + @property + def input_spectrum(self): + return self._input_spectrum + + @input_spectrum.setter + def input_spectrum(self, new_spectrum): + # We want to clear the refined locations if a new calibration spectrum is provided + self._clear_cache() + for line in self.lines: + line.input_spectrum = new_spectrum + self._input_spectrum = new_spectrum + @property def refined_lines(self): return [line.with_refined_pixel for line in self.lines] From 8f1535b1d1693e4fc6f60104b27f0261a213338c Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 2 Feb 2023 11:21:24 -0500 Subject: [PATCH 15/31] Add default range if method is set but not range, remove check for unimplemented method --- specreduce/tests/test_wavelength_calibration.py | 7 ------- specreduce/wavelength_calibration.py | 15 ++++----------- 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 75e99c2a..8e26fe1d 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -39,13 +39,6 @@ def test_poly_from_calibrationline(spec1d): def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line): - with pytest.raises(ValueError, match="You must define 'range' in refinement_kwargs"): - line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, - refinement_method='gaussian') - - with pytest.raises(ValueError, match="You must define 'direction' in refinement_kwargs"): - line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, - refinement_method='gradient') line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, refinement_method='gaussian', refinement_kwargs={'range': 25}) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 83f69367..961a4c2e 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -41,11 +41,8 @@ def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, if self.refinement_method in ("gaussian", "min", "max"): if 'range' not in self.refinement_kwargs: - raise ValueError(f"You must define 'range' in refinement_kwargs to use " - f"{self.refinement_method} refinement.") - elif self.refinement_method == "gradient" and 'direction' not in self.refinement_kwargs: - raise ValueError("You must define 'direction' in refinement_kwargs to use " - "gradient refinement") + # We may want to adjust this default based on real calibration spectra + self.refinement_kwargs['range'] = 10 def _clear_cache(self, *attrs): """ @@ -170,12 +167,8 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs if self.default_refinement_method in ("gaussian", "min", "max"): if 'range' not in self.default_refinement_kwargs: - raise ValueError("You must define 'range' in default_refinement_kwargs to use " - f"{self.refinement_method} refinement.") - elif (self.default_refinement_method == "gradient" and 'direction' not in - self.default_refinement_kwargs): - raise ValueError("You must define 'direction' in default_refinement_kwargs to use " - "gradient refinement") + # We may want to adjust this default based on real calibration spectra + self.default_refinement_method['range'] = 10 self.lines = [] if isinstance(lines, str): From 22e2d68d8e4f3b19b18c2560c0978d5d1e689349 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 3 Feb 2023 14:46:11 -0500 Subject: [PATCH 16/31] Also clear cache when updating lines or model --- .../tests/test_wavelength_calibration.py | 1 - specreduce/wavelength_calibration.py | 32 +++++++++++++++---- 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 8e26fe1d..86005b82 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -1,5 +1,4 @@ from numpy.testing import assert_allclose -import pytest import astropy.units as u from astropy.modeling.models import Polynomial1D diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 961a4c2e..43daaea6 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -151,7 +151,7 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs words """ self._input_spectrum = input_spectrum - self.model = model + self._model = model self.spectral_unit = spectral_unit self.default_refinement_method = default_refinement_method self.default_refinement_kwargs = default_refinement_kwargs @@ -170,18 +170,18 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs # We may want to adjust this default based on real calibration spectra self.default_refinement_method['range'] = 10 - self.lines = [] + self._lines = [] if isinstance(lines, str): if lines in self.available_line_lists: raise ValueError(f"Line list '{lines}' is not an available line list.") else: for line in lines: if isinstance(line, CalibrationLine): - self.lines.append(line) + self._lines.append(line) else: - self.lines.append(CalibrationLine(self.input_spectrum, line[0], line[1], - self.default_refinement_method, - self.default_refinement_kwargs)) + self._lines.append(CalibrationLine(self.input_spectrum, line[0], line[1], + self.default_refinement_method, + self.default_refinement_kwargs)) def _clear_cache(self, *attrs): """ @@ -205,6 +205,24 @@ def input_spectrum(self, new_spectrum): line.input_spectrum = new_spectrum self._input_spectrum = new_spectrum + @property + def lines(self): + return self._lines + + @lines.setter + def lines(self, new_lines): + self._clear_cache() + self._lines = new_lines + + @property + def model(self): + return self._model + + @model.setter + def model(self, new_model): + self._clear_cache() + self._model = new_model + @property def refined_lines(self): return [line.with_refined_pixel for line in self.lines] @@ -221,7 +239,7 @@ def wcs(self): y = [line[0].value for line in self.refined_pixels] * self.spectral_unit # Fit the model - self.model = self.fitter(self.model, x, y) + self._model = self.fitter(self._model, x, y) # Build a GWCS pipeline from the fitted model pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) From 4a508e31e0ddb70cd257308fa22fadb5ca0d4a75 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 3 Feb 2023 15:56:18 -0500 Subject: [PATCH 17/31] Move fitter default to wcs call --- specreduce/wavelength_calibration.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 43daaea6..4c614a60 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -156,14 +156,7 @@ def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angs self.default_refinement_method = default_refinement_method self.default_refinement_kwargs = default_refinement_kwargs self._cached_properties = ['wcs',] - - if fitter is None: - if self.model.linear: - self.fitter = LinearLSQFitter(calc_uncertainties=True) - else: - self.fitter = LMLSQFitter(calc_uncertainties=True) - else: - self.fitter = fitter + self.fitter = fitter if self.default_refinement_method in ("gaussian", "min", "max"): if 'range' not in self.default_refinement_kwargs: @@ -238,6 +231,12 @@ def wcs(self): x = [line[1] for line in self.refined_pixels] * u.pix y = [line[0].value for line in self.refined_pixels] * self.spectral_unit + if self.fitter is None: + if self.model.linear: + self.fitter = LinearLSQFitter(calc_uncertainties=True) + else: + self.fitter = LMLSQFitter(calc_uncertainties=True) + # Fit the model self._model = self.fitter(self._model, x, y) From a8ca23fed824690d52c2ea0145ce408a57f3d6e6 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 3 Feb 2023 15:57:55 -0500 Subject: [PATCH 18/31] Keep storing fitter as None --- specreduce/wavelength_calibration.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 4c614a60..93ea1c78 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -232,13 +232,16 @@ def wcs(self): y = [line[0].value for line in self.refined_pixels] * self.spectral_unit if self.fitter is None: + # Flexible defaulting if self.fitter is None if self.model.linear: - self.fitter = LinearLSQFitter(calc_uncertainties=True) + fitter = LinearLSQFitter(calc_uncertainties=True) else: - self.fitter = LMLSQFitter(calc_uncertainties=True) + fitter = LMLSQFitter(calc_uncertainties=True) + else: + fitter = self.fitter # Fit the model - self._model = self.fitter(self._model, x, y) + self._model = fitter(self._model, x, y) # Build a GWCS pipeline from the fitted model pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) From 1941a2650618f106265a2cd24dc1466f5b0a138b Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Tue, 7 Feb 2023 12:11:52 -0500 Subject: [PATCH 19/31] Add docs, improve importing of new classes --- docs/conf.py | 3 +- docs/index.rst | 1 + docs/wavelength_calibration.rst | 71 +++++++++++++++++++ specreduce/__init__.py | 1 + .../tests/test_wavelength_calibration.py | 2 +- specreduce/wavelength_calibration.py | 3 + 6 files changed, 79 insertions(+), 2 deletions(-) create mode 100644 docs/wavelength_calibration.rst diff --git a/docs/conf.py b/docs/conf.py index 7a668af0..ea4f7f0b 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -174,7 +174,8 @@ { 'astropy': ('https://docs.astropy.org/en/stable/', None), 'ccdproc': ('https://ccdproc.readthedocs.io/en/stable/', None), - 'specutils': ('https://specutils.readthedocs.io/en/stable/', None) + 'specutils': ('https://specutils.readthedocs.io/en/stable/', None), + 'gwcs': ('https://gwcs.readthedocs.io/en/stable/', None) } ) # diff --git a/docs/index.rst b/docs/index.rst index e05ac72b..785b4397 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -50,6 +50,7 @@ Calibration .. toctree:: :maxdepth: 1 + wavelength_calibration.rst extinction.rst specphot_standards.rst diff --git a/docs/wavelength_calibration.rst b/docs/wavelength_calibration.rst new file mode 100644 index 00000000..5b05e4ed --- /dev/null +++ b/docs/wavelength_calibration.rst @@ -0,0 +1,71 @@ +.. _wavelength_calibration: + +Wavelength Calibration +====================== + +Wavelength calibration is currently supported for 1D spectra. Given a list of spectral +lines with known wavelengths and estimated pixel positions on an input calibration +spectrum, you can currently use ``specreduce`` to: + +#. Refine the pixel position estimates of the lines based on features in the input + calibration spectrum. +#. Fit an ``astropy`` model to the wavelength/pixel pairs to generate a spectral WCS + solution for the dispersion. +#. Apply the generated spectral WCS to other `~specutils.Spectrum1D` objects. + +Calibration Lines +----------------- + +``specreduce`` provides a `~specreduce.wavelength_calibration.CalibrationLine` class for +encoding the information about each spectral line to be used in the dispersion model +solution fitting. The minimum information required is the wavelength and estimated pixel +location on the calibration spectrum of each line. By default, no refinement of the pixel +position is done; to enable this feature you must provide a refinement method and any +associated keywords (currently ``range`` defines the number of pixels on either side +of the estimated location to use for all available refinement options). For example, +to define a line and use ``specreduce`` to update the pixel location the the centroid +of a gaussian fit to the line, you would do the following:: + + import astropy.units as u + from specreduce import CalibrationLine + test_line = CalibrationLine(calibration_spectrum, 6562.81*u.AA, 500, + refinement_method="gaussian", + refinement_kwargs={"range": 10}) + refined_line = test_line.with_refined_pixel() + +Note that in this example and in the example below, ``calibration_spectrum`` is a +`~specutils.Spectrum1D` spectrum to be used to refine the line locations (for example +a lamp spectrum). + +1D Wavelength Calibration +------------------------- + +Of course, refining the location of a line is only really useful when done for multiple +lines to be fit as part of the dispersion model. This can be accomplished with the +`~specreduce.wavelength_calibration.WavelengthCalibration1D` class as follows below. +Note that the lines can be input as simple (wavelength, pixel) pairs rather than as +`~specreduce.wavelength_calibration.CalibrationLine` objects - in this case the +``default_refinement_method`` and ``default_refinement_kwargs`` will be applied to +each input line:: + + import astropy.units as u + from specreduce import WavelengthCalibration1D + test_cal = WavelengthCalibration1D(calibration_spectrum, + [(5000*u.AA, 0), (5100*u.AA, 10), + (5198*u.AA, 20), (5305*u.AA, 30)], + default_refinement_method="gaussian") + calibrated_spectrum = test_cal.apply_to_spectrum(science_spectrum) + +The example above uses the default model (`~astropy.modeling.functional_models.Linear1D`) +to fit the input spectral lines, and then applies the calculated WCS solution to a second +spectrum (``science_spectrum``). Any other ``astropy`` model can be provided as the +input ``model`` parameter to the `~specreduce.wavelength_calibration.WavelengthCalibration1D`. +In the above example, the model fit and WCS construction is all done as part of the +``apply_to_spectrum()`` call, but you could also access the `~gwcs.wcs.WCS` object itself +by calling:: + + test_cal.wcs + +The calculated WCS is a cached property that will be cleared if the ``lines``, ``model``, +or ``input_spectrum`` properties are updated, since these will alter the calculated dispersion +fit. \ No newline at end of file diff --git a/specreduce/__init__.py b/specreduce/__init__.py index 6bff38bb..5bf2362e 100644 --- a/specreduce/__init__.py +++ b/specreduce/__init__.py @@ -7,3 +7,4 @@ # ---------------------------------------------------------------------------- from specreduce.core import * # noqa +from specreduce.wavelength_calibration import * # noqa diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 86005b82..f186c1da 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -4,7 +4,7 @@ from astropy.modeling.models import Polynomial1D from astropy.tests.helper import assert_quantity_allclose -from specreduce.wavelength_calibration import CalibrationLine, WavelengthCalibration1D +from specreduce import CalibrationLine, WavelengthCalibration1D def test_linear_from_list(spec1d): diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 93ea1c78..01eaeddc 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -9,6 +9,9 @@ from specutils.fitting import fit_lines +__all__ = ['CalibrationLine', 'WavelengthCalibration1D'] + + class CalibrationLine(): def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, From 829dc5cba1b2a9e1fddad3483b2ca06a6f571788 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 6 Apr 2023 17:08:26 -0400 Subject: [PATCH 20/31] Remove CalibrationLine, move to astropy tables internallyt --- .../tests/test_wavelength_calibration.py | 55 ++-- specreduce/wavelength_calibration.py | 266 ++++++------------ 2 files changed, 112 insertions(+), 209 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index f186c1da..521c373f 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -1,69 +1,52 @@ from numpy.testing import assert_allclose +from astropy.table import QTable import astropy.units as u from astropy.modeling.models import Polynomial1D from astropy.tests.helper import assert_quantity_allclose -from specreduce import CalibrationLine, WavelengthCalibration1D +from specreduce import WavelengthCalibration1D def test_linear_from_list(spec1d): - test = WavelengthCalibration1D(spec1d, [(5000*u.AA, 0), (5100*u.AA, 10), - (5198*u.AA, 20), (5305*u.AA, 30)]) + centers = [0, 10, 20, 30] + w = [5000, 5100, 5198, 5305]*u.AA + test = WavelengthCalibration1D(spec1d, centers, line_wavelengths=w) spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) - -def test_linear_from_calibrationline(spec1d): - lines = [CalibrationLine(spec1d, 5000*u.AA, 0), CalibrationLine(spec1d, 5100*u.AA, 10), - CalibrationLine(spec1d, 5198*u.AA, 20), CalibrationLine(spec1d, 5305*u.AA, 30)] - test = WavelengthCalibration1D(spec1d, lines) +def test_linear_from_table(spec1d): + centers = [0, 10, 20, 30] + w = [5000, 5100, 5198, 5305]*u.AA + table = QTable([centers, w], names=["pixel_center", "wavelength"]) + test = WavelengthCalibration1D(spec1d, table) spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) - -def test_poly_from_calibrationline(spec1d): +def test_poly_from_table(spec1d): # This test is mostly to prove that you can use other models - lines = [CalibrationLine(spec1d, 5005*u.AA, 0), CalibrationLine(spec1d, 5110*u.AA, 10), - CalibrationLine(spec1d, 5214*u.AA, 20), CalibrationLine(spec1d, 5330*u.AA, 30), - CalibrationLine(spec1d, 5438*u.AA, 40)] - test = WavelengthCalibration1D(spec1d, lines, model=Polynomial1D(2)) + centers = [0, 10, 20, 30, 40] + w = [5005, 5110, 5214, 5330, 5438]*u.AA + table = QTable([centers, w], names=["pixel_center", "wavelength"]) + + test = WavelengthCalibration1D(spec1d, table, model=Polynomial1D(2)) test.apply_to_spectrum(spec1d) assert_allclose(test.model.parameters, [5.00477143e+03, 1.03457143e+01, 1.28571429e-02]) -def test_calibrationline(spec1d_with_emission_line, spec1d_with_absorption_line): - - line = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 128, refinement_method='gaussian', - refinement_kwargs={'range': 25}) - assert_allclose(line.refine(), 129.44371, atol=0.01) - - line2 = CalibrationLine(spec1d_with_emission_line, 5000*u.AA, 130, refinement_method='max', - refinement_kwargs={'range': 10}) - assert line2.refine() == 128 - - line3 = CalibrationLine(spec1d_with_absorption_line, 5000*u.AA, 128, refinement_method='min', - refinement_kwargs={'range': 10}) - assert line3.refine() == 130 - - def test_replace_spectrum(spec1d, spec1d_with_emission_line): - lines = [CalibrationLine(spec1d, 5000*u.AA, 0), CalibrationLine(spec1d, 5100*u.AA, 10), - CalibrationLine(spec1d, 5198*u.AA, 20), CalibrationLine(spec1d, 5305*u.AA, 30)] - test = WavelengthCalibration1D(spec1d, lines) + centers = [0, 10, 20, 30]*u.pix + w = [5000, 5100, 5198, 5305]*u.AA + test = WavelengthCalibration1D(spec1d, centers, line_wavelengths=w) # Accessing this property causes fits the model and caches the resulting WCS test.wcs assert "wcs" in test.__dict__ - for line in test.lines: - assert "refined_pixel" in line.__dict__ # Replace the input spectrum, which should clear the cached properties test.input_spectrum = spec1d_with_emission_line assert "wcs" not in test.__dict__ - for line in test.lines: - assert "refined_pixel" not in line.__dict__ diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 01eaeddc..d0a9ab22 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -1,183 +1,119 @@ -from astropy.modeling.models import Linear1D, Gaussian1D +from astropy.modeling.models import Linear1D from astropy.modeling.fitting import LMLSQFitter, LinearLSQFitter +from astropy.table import QTable, hstack import astropy.units as u from functools import cached_property from gwcs import wcs from gwcs import coordinate_frames as cf import numpy as np -from specutils import SpectralRegion, Spectrum1D -from specutils.fitting import fit_lines +from specutils import Spectrum1D -__all__ = ['CalibrationLine', 'WavelengthCalibration1D'] +__all__ = ['WavelengthCalibration1D'] -class CalibrationLine(): +def get_available_catalogs(): + """ + ToDo: Decide in what format to store calibration line catalogs (e.g., for lamps) + and write this function to determine the list of available catalog names. + """ + return [] - def __init__(self, input_spectrum, wavelength, pixel, refinement_method=None, - refinement_kwargs={}): - """ - input_spectrum: `~specutils.Spectrum1D` - A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. - wavelength: `astropy.units.Quantity` - The rest wavelength of the line. - pixel: int - Initial guess of the integer pixel location of the line center in the spectrum. - refinement_method: str, optional - None: use the actual provided pixel values as anchors for the fit without refining - their location. - 'gradient': Simple gradient descent/ascent to nearest min/max value. Direction - defined by ``direction`` in refinement_kwargs. - 'min'/'max': find min/max within pixel range provided by ``range`` kwarg. - 'gaussian': fit gaussian to spectrum within pixel range provided by ``range`` kwarg. - refinement_kwargs: dict, optional - Keywords to set the parameters for line location refinement. Distance on either side - of line center to include on either side of gaussian fit is set by int ``range``. - Gradient descent/ascent is determined by setting ``direction`` to 'min' or 'max'. - """ - self._input_spectrum = input_spectrum - self.wavelength = wavelength - self.pixel = pixel - self.refinement_method = refinement_method - self.refinement_kwargs = refinement_kwargs - self._cached_properties = ['refined_pixel', 'with_refined_pixel'] - - if self.refinement_method in ("gaussian", "min", "max"): - if 'range' not in self.refinement_kwargs: - # We may want to adjust this default based on real calibration spectra - self.refinement_kwargs['range'] = 10 - def _clear_cache(self, *attrs): - """ - provide convenience function to clearing the cache for cached_properties - """ - if not len(attrs): - attrs = self._cached_properties - for attr in attrs: - if attr in self.__dict__: - del self.__dict__[attr] - - @classmethod - def by_name(cls, input_spectrum, line_name, pixel, - refinement_method=None, refinement_kwargs={}): - # this will do wavelength lookup and passes that wavelength to __init__ allowing the user - # to call CalibrationLine.by_name(spectrum, 'H_alpha', pixel=40) - # We could also have a type-check for ``wavelength`` above, this just feels more verbose - return cls(...) - - def _do_refinement(self, input_spectrum): - if self.refinement_method == 'gaussian': - window_width = self.refinement_kwargs.get('range') - window = SpectralRegion((self.pixel-window_width)*u.pix, - (self.pixel+window_width)*u.pix) - - input_model = Gaussian1D(mean=self.pixel, stddev=3, - amplitude=self.input_spectrum.flux[self.pixel]) - - fitted_model = fit_lines(self.input_spectrum, input_model, window=window) - new_pixel = fitted_model.mean.value - - elif self.refinement_method == 'min': - window_width = self.refinement_kwargs.get('range') - new_pixel = np.argmin(self.input_spectrum.flux[self.pixel-window_width: - self.pixel+window_width+1]) - new_pixel += self.pixel - window_width - - elif self.refinement_method == 'max': - window_width = self.refinement_kwargs.get('range') - new_pixel = np.argmax(self.input_spectrum.flux[self.pixel-window_width: - self.pixel+window_width+1]) - new_pixel += self.pixel - window_width - - else: - raise ValueError(f"Refinement method {self.refinement_method} is not implemented") - - return new_pixel - - def refine(self, input_spectrum=None, return_object=False): - # finds the center of the line according to refinement_method/kwargs and returns - # either the pixel value or a CalibrationLine object with pixel changed to the refined value - if self.refinement_method is None: - return self.pixel - input_spectrum = self.input_spectrum if input_spectrum is None else input_spectrum - refined_pixel = self._do_refinement(input_spectrum) - if return_object: - return CalibrationLine(input_spectrum, self.wavelength, refined_pixel, - refinement_method=self.refinement_method, - refinement_kwargs=self.refinement_kwargs) - - return refined_pixel - - @cached_property - def refined_pixel(self): - # returns the refined pixel for self.input_spectrum - return self.refine() - - @cached_property - def with_refined_pixel(self): - # returns a copy of this object, but with the pixel updated to the refined value - return self.refine(return_object=True) - - @property - def input_spectrum(self): - return self._input_spectrum - - @input_spectrum.setter - def input_spectrum(self, new_spectrum): - # We want to clear the refined locations if a new calibration spectrum is provided - self._clear_cache() - self._input_spectrum = new_spectrum - - def __str__(self): - return f"CalibrationLine: ({self.wavelength}, {self.pixel})" - - def __repr__(self): - return f"CalibrationLine({self.wavelength}, {self.pixel})" +def concatenate_catalogs(): + """ + ToDo: Code logic to combine the lines from multiple catalogs if needed + """ + pass class WavelengthCalibration1D(): - def __init__(self, input_spectrum, lines, model=Linear1D(), spectral_unit=u.Angstrom, - fitter=None, default_refinement_method=None, default_refinement_kwargs={}): + def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=None, + model=Linear1D(), spectral_unit=u.Angstrom, fitter=None): """ input_spectrum: `~specutils.Spectrum1D` A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. - lines: str, list - List of lines to anchor the wavelength solution fit. List items are coerced to - CalibrationLine objects if passed as tuples of (wavelength, pixel) with - default_refinement_method and default_refinement_kwargs as defaults. + line_list: list, array, `~astropy.table.QTable` + List or array of line pixel locations to anchor the wavelength solution fit. + Will be converted to an astropy table internally if a list or array was input. + Can also be input as an `~astropy.table.QTable` table with (minimally) a column + named "pixel_center" and optionally a "wavelength" column with known line + wavelengths populated. + line_wavelengths: `~astropy.units.Quantity`, `~astropy.table.QTable`, optional + `astropy.units.Quantity` array of line wavelength values corresponding to the + line_pixels. Does not have to be in the same order (the lists will be sorted) + but does currently need to be the same length as line_list. Can also be input + as an `~astropy.table.QTable` with (minimally) a "wavelength" column. + catalog: list, str, optional + The name of a catalog of line wavelengths to load and use in automated and + template-matching line matching. model: `~astropy.modeling.Model` The model to fit for the wavelength solution. Defaults to a linear model. - default_refinement_method: str, optional - words - default_refinement_kwargs: dict, optional - words """ self._input_spectrum = input_spectrum self._model = model + self._line_list = line_list self.spectral_unit = spectral_unit - self.default_refinement_method = default_refinement_method - self.default_refinement_kwargs = default_refinement_kwargs self._cached_properties = ['wcs',] self.fitter = fitter + self._potential_wavelengths = None + self._catalog = catalog + + # ToDo: Implement having line catalogs + self._available_catalogs = get_available_catalogs() + + if isinstance(line_list, (list, np.ndarray)): + self._line_list = QTable([line_list], names=["pixel_center"]) + + if self._line_list["pixel_center"].unit is None: + self._line_list["pixel_center"].unit = u.pix + + # Make sure our pixel locations are sorted + self._line_list.sort("pixel_center") + + if (line_wavelengths is None and catalog is None + and "wavelength" not in self._line_list.columns): + raise ValueError("You must specify at least one of line_wavelengths, " + "catalog, or 'wavelength' column in line_list.") + + # Sanity checks on line_wavelengths value + if line_wavelengths is not None: + if "wavelength" in line_list: + raise ValueError("Cannot specify line_wavelengths separately if there is" + "a 'wavelength' column in line_list.") + if len(line_wavelengths) != len(line_list): + raise ValueError("If line_wavelengths is specified, it must have the same " + "length as line_pixels") + if not isinstance(line_wavelengths, (u.Quantity, QTable)): + raise ValueError("line_wavelengths must be specified as an astropy.units.Quantity" + "array or as an astropy.table.QTable") + if isinstance(line_wavelengths, u.Quantity): + line_wavelengths.value.sort() + self._line_list["wavelength"] = line_wavelengths + elif isinstance(line_wavelengths, QTable): + line_wavelengths.sort("wavelength") + self._line_list = hstack(self._line_list, line_wavelengths) + + if catalog is not None: + if isinstance(catalog, list): + self._catalog = catalog + else: + self._catalog = [catalog] + for cat in catalog: + if isinstance(cat, str): + if cat not in self._available_catalogs: + raise ValueError(f"Line list '{cat}' is not an available catalog.") - if self.default_refinement_method in ("gaussian", "min", "max"): - if 'range' not in self.default_refinement_kwargs: - # We may want to adjust this default based on real calibration spectra - self.default_refinement_method['range'] = 10 + # Get the potential lines from any specified catalogs to use in matching + self._potential_wavelengths = concatenate_catalogs(self._catalog) - self._lines = [] - if isinstance(lines, str): - if lines in self.available_line_lists: - raise ValueError(f"Line list '{lines}' is not an available line list.") - else: - for line in lines: - if isinstance(line, CalibrationLine): - self._lines.append(line) - else: - self._lines.append(CalibrationLine(self.input_spectrum, line[0], line[1], - self.default_refinement_method, - self.default_refinement_kwargs)) + def identify_lines(self): + """ + ToDo: Code matching algorithm between line pixel locations and potential line + wavelengths from catalogs. + """ + pass def _clear_cache(self, *attrs): """ @@ -189,6 +125,10 @@ def _clear_cache(self, *attrs): if attr in self.__dict__: del self.__dict__[attr] + @property + def available_catalogs(self): + return self._available_catalogs + @property def input_spectrum(self): return self._input_spectrum @@ -197,19 +137,8 @@ def input_spectrum(self): def input_spectrum(self, new_spectrum): # We want to clear the refined locations if a new calibration spectrum is provided self._clear_cache() - for line in self.lines: - line.input_spectrum = new_spectrum self._input_spectrum = new_spectrum - @property - def lines(self): - return self._lines - - @lines.setter - def lines(self, new_lines): - self._clear_cache() - self._lines = new_lines - @property def model(self): return self._model @@ -219,20 +148,11 @@ def model(self, new_model): self._clear_cache() self._model = new_model - @property - def refined_lines(self): - return [line.with_refined_pixel for line in self.lines] - - @property - def refined_pixels(self): - # useful for plotting over spectrum to change input lines/refinement options - return [(line.wavelength, line.refined_pixel) for line in self.lines] - @cached_property def wcs(self): # computes and returns WCS after fitting self.model to self.refined_pixels - x = [line[1] for line in self.refined_pixels] * u.pix - y = [line[0].value for line in self.refined_pixels] * self.spectral_unit + x = self._line_list["pixel_center"] + y = self._line_list["wavelength"] if self.fitter is None: # Flexible defaulting if self.fitter is None From c3ba30d97b25f73cf2caf973dcaf5af1924f96bc Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 6 Apr 2023 17:09:43 -0400 Subject: [PATCH 21/31] Codestyle --- specreduce/tests/test_wavelength_calibration.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 521c373f..40beb899 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -17,6 +17,7 @@ def test_linear_from_list(spec1d): assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) + def test_linear_from_table(spec1d): centers = [0, 10, 20, 30] w = [5000, 5100, 5198, 5305]*u.AA @@ -27,6 +28,7 @@ def test_linear_from_table(spec1d): assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) + def test_poly_from_table(spec1d): # This test is mostly to prove that you can use other models centers = [0, 10, 20, 30, 40] From 62c802587758199fd0174845f5eeca9689880c36 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen <39831871+rosteen@users.noreply.github.com> Date: Fri, 7 Apr 2023 11:29:45 -0400 Subject: [PATCH 22/31] Apply suggestions from code review Co-authored-by: Kyle Conroy --- specreduce/wavelength_calibration.py | 44 ++-------------------------- 1 file changed, 2 insertions(+), 42 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index d0a9ab22..38fd3f53 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -42,7 +42,7 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non wavelengths populated. line_wavelengths: `~astropy.units.Quantity`, `~astropy.table.QTable`, optional `astropy.units.Quantity` array of line wavelength values corresponding to the - line_pixels. Does not have to be in the same order (the lists will be sorted) + line pixels defined in ``line_list``. Does not have to be in the same order (the lists will be sorted) but does currently need to be the same length as line_list. Can also be input as an `~astropy.table.QTable` with (minimally) a "wavelength" column. catalog: list, str, optional @@ -100,7 +100,7 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non self._catalog = catalog else: self._catalog = [catalog] - for cat in catalog: + for cat in self._catalog: if isinstance(cat, str): if cat not in self._available_catalogs: raise ValueError(f"Line list '{cat}' is not an available catalog.") @@ -186,43 +186,3 @@ def apply_to_spectrum(self, spectrum=None): return updated_spectrum -''' -# WavelengthCalibration2D is a planned future feature -class WavelengthCalibration2D(): - # input_spectrum must be 2-dimensional - - # lines are coerced to CalibrationLine objects if passed as tuples with - # default_refinement_method and default_refinement_kwargs as defaults - def __init__(input_spectrum, trace, lines, model=Linear1D, - default_refinement_method=None, default_refinement_kwargs={}): - return NotImplementedError("2D wavelength calibration is not yet implemented") - - @classmethod - def autoidentify(cls, input_spectrum, trace, line_list, model=Linear1D): - # does this do 2D identification, or just search a 1d spectrum exactly at the trace? - # or should we require using WavelengthCalibration1D.autoidentify? - return cls(...) - - @property - def refined_lines(self): - return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=False)) - for l in self.lines] for row in rows] - - @property - def refined_pixels(self): - return [[(l.wavelength, l.refine(input_spectrum.get_row(row), return_object=True)) - for l in self.lines] for row in rows] - - @cached_property - def wcs(self): - # computes and returns WCS after fitting self.model to self.refined_pixels - pass - - def __call__(self, apply_to_spectrum=None): - # returns spectrum1d with wavelength calibration applied - # actual line refinement and WCS solution should already be done so that this - # can be called on multiple science sources - apply_to_spectrum = self.input_spectrum if apply_to_spectrum is None else apply_to_spectrum - apply_to_spectrum.wcs = apply_to_spectrum # might need a deepcopy! - return apply_to_spectrum -''' From 12a8661856b985587f0f75161a467bf77423b5e0 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 7 Apr 2023 11:51:04 -0400 Subject: [PATCH 23/31] Remove speactral_units arg, add docstring for fitter arg --- specreduce/wavelength_calibration.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 38fd3f53..f406c900 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -30,7 +30,7 @@ def concatenate_catalogs(): class WavelengthCalibration1D(): def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=None, - model=Linear1D(), spectral_unit=u.Angstrom, fitter=None): + model=Linear1D(), fitter=None): """ input_spectrum: `~specutils.Spectrum1D` A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. @@ -50,11 +50,14 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non template-matching line matching. model: `~astropy.modeling.Model` The model to fit for the wavelength solution. Defaults to a linear model. + fitter: `~astropy.modeling.fitting.Fitter`, optional + The fitter to use in optimizing the model fit. Defaults to + `~astropy.modeling.fitting.LinearLSQFitter` if the model to fit is linear + or `~astropy.modeling.fitting.LMLSQFitter` if the model to fit is non-linear. """ self._input_spectrum = input_spectrum self._model = model self._line_list = line_list - self.spectral_unit = spectral_unit self._cached_properties = ['wcs',] self.fitter = fitter self._potential_wavelengths = None @@ -168,7 +171,8 @@ def wcs(self): # Build a GWCS pipeline from the fitted model pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) - spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], unit=[self.spectral_unit,]) + spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], + unit=[self._line_list["wavelength"].unit,]) pipeline = [(pixel_frame, self.model), (spectral_frame, None)] From 5063b5ec3d13e7070836ff749cdd9779158c4c06 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 7 Apr 2023 13:34:48 -0400 Subject: [PATCH 24/31] Update documentation, add NotImplementedError for catalog input --- docs/wavelength_calibration.rst | 67 ++++++++++------------------ specreduce/wavelength_calibration.py | 5 +++ 2 files changed, 29 insertions(+), 43 deletions(-) diff --git a/docs/wavelength_calibration.rst b/docs/wavelength_calibration.rst index 5b05e4ed..28a5caa1 100644 --- a/docs/wavelength_calibration.rst +++ b/docs/wavelength_calibration.rst @@ -7,53 +7,25 @@ Wavelength calibration is currently supported for 1D spectra. Given a list of sp lines with known wavelengths and estimated pixel positions on an input calibration spectrum, you can currently use ``specreduce`` to: -#. Refine the pixel position estimates of the lines based on features in the input - calibration spectrum. #. Fit an ``astropy`` model to the wavelength/pixel pairs to generate a spectral WCS solution for the dispersion. #. Apply the generated spectral WCS to other `~specutils.Spectrum1D` objects. -Calibration Lines ------------------ - -``specreduce`` provides a `~specreduce.wavelength_calibration.CalibrationLine` class for -encoding the information about each spectral line to be used in the dispersion model -solution fitting. The minimum information required is the wavelength and estimated pixel -location on the calibration spectrum of each line. By default, no refinement of the pixel -position is done; to enable this feature you must provide a refinement method and any -associated keywords (currently ``range`` defines the number of pixels on either side -of the estimated location to use for all available refinement options). For example, -to define a line and use ``specreduce`` to update the pixel location the the centroid -of a gaussian fit to the line, you would do the following:: - - import astropy.units as u - from specreduce import CalibrationLine - test_line = CalibrationLine(calibration_spectrum, 6562.81*u.AA, 500, - refinement_method="gaussian", - refinement_kwargs={"range": 10}) - refined_line = test_line.with_refined_pixel() - -Note that in this example and in the example below, ``calibration_spectrum`` is a -`~specutils.Spectrum1D` spectrum to be used to refine the line locations (for example -a lamp spectrum). - 1D Wavelength Calibration ------------------------- -Of course, refining the location of a line is only really useful when done for multiple -lines to be fit as part of the dispersion model. This can be accomplished with the -`~specreduce.wavelength_calibration.WavelengthCalibration1D` class as follows below. -Note that the lines can be input as simple (wavelength, pixel) pairs rather than as -`~specreduce.wavelength_calibration.CalibrationLine` objects - in this case the -``default_refinement_method`` and ``default_refinement_kwargs`` will be applied to -each input line:: - - import astropy.units as u - from specreduce import WavelengthCalibration1D - test_cal = WavelengthCalibration1D(calibration_spectrum, - [(5000*u.AA, 0), (5100*u.AA, 10), - (5198*u.AA, 20), (5305*u.AA, 30)], - default_refinement_method="gaussian") +The `~specreduce.wavelength_calibration.WavelengthCalibration1D` class can be used +to fit a dispersion model to a list of line positions and wavelengths. Future development +will implement catalogs of known lamp spectra for use in matching observed lines. In the +example below, the line positions (``line_list``) have already been extracted from +``lamp_spectrum``:: + + import astropy.units as u + from specreduce import WavelengthCalibration1D + line_list = [10, 22, 31, 43] + wavelengths = [5340, 5410, 5476, 5543]*u.AA + test_cal = WavelengthCalibration1D(lamp_spectrum, line_list, + line_wavelengths=wavelengths) calibrated_spectrum = test_cal.apply_to_spectrum(science_spectrum) The example above uses the default model (`~astropy.modeling.functional_models.Linear1D`) @@ -64,8 +36,17 @@ In the above example, the model fit and WCS construction is all done as part of ``apply_to_spectrum()`` call, but you could also access the `~gwcs.wcs.WCS` object itself by calling:: - test_cal.wcs + test_cal.wcs -The calculated WCS is a cached property that will be cleared if the ``lines``, ``model``, +The calculated WCS is a cached property that will be cleared if the ``line_list``, ``model``, or ``input_spectrum`` properties are updated, since these will alter the calculated dispersion -fit. \ No newline at end of file +fit. + +You can also provide the input pixel locations and wavelengths of the lines as an +`~astropy.table.QTable`, with columns ``pixel_center`` and ``wavelength``:: + + from astropy.table import QTable + pixels = [10, 20, 30, 40]*u.pix + wavelength = [5340, 5410, 5476, 5543]*u.AA + line_list = QTable([pixels, wavelength], names=["pixel_center", "wavelength"]) + test_cal = WavelengthCalibration1D(lamp_spectrum, line_list) \ No newline at end of file diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index f406c900..e26a4021 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -98,7 +98,12 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non line_wavelengths.sort("wavelength") self._line_list = hstack(self._line_list, line_wavelengths) + # Parse desired catalogs of lines for matching. if catalog is not None: + # For now we avoid going into the later logic and just throw an error + raise NotImplementedError("No catalogs are available yet, please input " + "wavelengths with line_wavelengths or as a " + "column in line_list") if isinstance(catalog, list): self._catalog = catalog else: From b6a92f84a9bf928a9b24ddbf2598d615b2e71ee4 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 7 Apr 2023 16:42:05 -0400 Subject: [PATCH 25/31] Reverse sort frequencies --- specreduce/wavelength_calibration.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index e26a4021..04ff382f 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -92,7 +92,11 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non raise ValueError("line_wavelengths must be specified as an astropy.units.Quantity" "array or as an astropy.table.QTable") if isinstance(line_wavelengths, u.Quantity): - line_wavelengths.value.sort() + # Ensure frequency is descending or wavelength is ascending + if str(line_wavelengths.unit.physical_type) == "frequency": + line_wavelengths[::-1].sort() + else: + line_wavelengths.sort() self._line_list["wavelength"] = line_wavelengths elif isinstance(line_wavelengths, QTable): line_wavelengths.sort("wavelength") From 1aa62c6d9f218c8e583fe9bd3c407522c9cb7a5f Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 7 Apr 2023 16:43:59 -0400 Subject: [PATCH 26/31] Codestyle --- specreduce/wavelength_calibration.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 04ff382f..912286de 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -42,9 +42,10 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non wavelengths populated. line_wavelengths: `~astropy.units.Quantity`, `~astropy.table.QTable`, optional `astropy.units.Quantity` array of line wavelength values corresponding to the - line pixels defined in ``line_list``. Does not have to be in the same order (the lists will be sorted) - but does currently need to be the same length as line_list. Can also be input - as an `~astropy.table.QTable` with (minimally) a "wavelength" column. + line pixels defined in ``line_list``. Does not have to be in the same order] + (the lists will be sorted) but does currently need to be the same length as + line_list. Can also be input as an `~astropy.table.QTable` with (minimally) + a "wavelength" column. catalog: list, str, optional The name of a catalog of line wavelengths to load and use in automated and template-matching line matching. @@ -197,5 +198,3 @@ def apply_to_spectrum(self, spectrum=None): updated_spectrum = Spectrum1D(spectrum.flux, wcs=self.wcs, mask=spectrum.mask, uncertainty=spectrum.uncertainty) return updated_spectrum - - From 7f31b61fe991838530e19c3e0b93643ae8ec00d2 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Mon, 10 Apr 2023 15:52:20 -0400 Subject: [PATCH 27/31] Increase test coverage --- .../tests/test_wavelength_calibration.py | 31 ++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 40beb899..07b74d6c 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -1,8 +1,10 @@ from numpy.testing import assert_allclose +import pytest from astropy.table import QTable import astropy.units as u from astropy.modeling.models import Polynomial1D +from astropy.modeling.fitting import LinearLSQFitter from astropy.tests.helper import assert_quantity_allclose from specreduce import WavelengthCalibration1D @@ -18,6 +20,13 @@ def test_linear_from_list(spec1d): assert_quantity_allclose(spec2.spectral_axis[-1], 5495.169999*u.AA) +def test_wavelength_from_table(spec1d): + centers = [0, 10, 20, 30] + w = [5000, 5100, 5198, 5305]*u.AA + table = QTable([w], names=["wavelength"]) + WavelengthCalibration1D(spec1d, centers, line_wavelengths=table) + + def test_linear_from_table(spec1d): centers = [0, 10, 20, 30] w = [5000, 5100, 5198, 5305]*u.AA @@ -35,7 +44,7 @@ def test_poly_from_table(spec1d): w = [5005, 5110, 5214, 5330, 5438]*u.AA table = QTable([centers, w], names=["pixel_center", "wavelength"]) - test = WavelengthCalibration1D(spec1d, table, model=Polynomial1D(2)) + test = WavelengthCalibration1D(spec1d, table, model=Polynomial1D(2), fitter=LinearLSQFitter()) test.apply_to_spectrum(spec1d) assert_allclose(test.model.parameters, [5.00477143e+03, 1.03457143e+01, 1.28571429e-02]) @@ -52,3 +61,23 @@ def test_replace_spectrum(spec1d, spec1d_with_emission_line): # Replace the input spectrum, which should clear the cached properties test.input_spectrum = spec1d_with_emission_line assert "wcs" not in test.__dict__ + + +def test_expected_errors(spec1d): + centers = [0, 10, 20, 30, 40] + w = [5005, 5110, 5214, 5330, 5438]*u.AA + table = QTable([centers, w], names=["pixel_center", "wavelength"]) + + with pytest.raises(ValueError, match="Cannot specify line_wavelengths separately"): + WavelengthCalibration1D(spec1d, table, line_wavelengths=w) + + with pytest.raises(ValueError, match="must have the same length"): + w2 = [5005, 5110, 5214, 5330, 5438, 5500]*u.AA + WavelengthCalibration1D(spec1d, table, line_wavelengths=w2) + + with pytest.raises(ValueError, match="must have the same length"): + w2 = [5005, 5110, 5214, 5330, 5438, 5500]*u.AA + WavelengthCalibration1D(spec1d, table, line_wavelengths=w2) + + with pytest.raises(ValueError, match="specify at least one"): + WavelengthCalibration1D(spec1d, centers) From 45939a9d0594ee05717e73a64d255baf052f740c Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Mon, 10 Apr 2023 16:07:49 -0400 Subject: [PATCH 28/31] Fix tests and bugs found by tests --- specreduce/tests/test_wavelength_calibration.py | 9 +++++---- specreduce/wavelength_calibration.py | 10 +++++----- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 07b74d6c..3caf7cc6 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -73,11 +73,12 @@ def test_expected_errors(spec1d): with pytest.raises(ValueError, match="must have the same length"): w2 = [5005, 5110, 5214, 5330, 5438, 5500]*u.AA - WavelengthCalibration1D(spec1d, table, line_wavelengths=w2) + WavelengthCalibration1D(spec1d, centers, line_wavelengths=w2) - with pytest.raises(ValueError, match="must have the same length"): - w2 = [5005, 5110, 5214, 5330, 5438, 5500]*u.AA - WavelengthCalibration1D(spec1d, table, line_wavelengths=w2) + with pytest.raises(ValueError, match="astropy.units.Quantity array or" + " as an astropy.table.QTable"): + w2 = [5005, 5110, 5214, 5330, 5438] + WavelengthCalibration1D(spec1d, centers, line_wavelengths=w2) with pytest.raises(ValueError, match="specify at least one"): WavelengthCalibration1D(spec1d, centers) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 912286de..91781b7d 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -46,7 +46,7 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non (the lists will be sorted) but does currently need to be the same length as line_list. Can also be input as an `~astropy.table.QTable` with (minimally) a "wavelength" column. - catalog: list, str, optional + catalog: list, str, `~astropy.table.QTable`, optional The name of a catalog of line wavelengths to load and use in automated and template-matching line matching. model: `~astropy.modeling.Model` @@ -83,15 +83,15 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non # Sanity checks on line_wavelengths value if line_wavelengths is not None: - if "wavelength" in line_list: + if isinstance(line_list, QTable) and "wavelength" in line_list.columns: raise ValueError("Cannot specify line_wavelengths separately if there is" - "a 'wavelength' column in line_list.") + " a 'wavelength' column in line_list.") if len(line_wavelengths) != len(line_list): raise ValueError("If line_wavelengths is specified, it must have the same " "length as line_pixels") if not isinstance(line_wavelengths, (u.Quantity, QTable)): raise ValueError("line_wavelengths must be specified as an astropy.units.Quantity" - "array or as an astropy.table.QTable") + " array or as an astropy.table.QTable") if isinstance(line_wavelengths, u.Quantity): # Ensure frequency is descending or wavelength is ascending if str(line_wavelengths.unit.physical_type) == "frequency": @@ -101,7 +101,7 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non self._line_list["wavelength"] = line_wavelengths elif isinstance(line_wavelengths, QTable): line_wavelengths.sort("wavelength") - self._line_list = hstack(self._line_list, line_wavelengths) + self._line_list = hstack([self._line_list, line_wavelengths]) # Parse desired catalogs of lines for matching. if catalog is not None: From f92bac589e698dda3148d23ccb82ed117a228a12 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Mon, 10 Apr 2023 16:12:46 -0400 Subject: [PATCH 29/31] Codestyle --- specreduce/wavelength_calibration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 91781b7d..19646ef7 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -85,7 +85,7 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non if line_wavelengths is not None: if isinstance(line_list, QTable) and "wavelength" in line_list.columns: raise ValueError("Cannot specify line_wavelengths separately if there is" - " a 'wavelength' column in line_list.") + " a 'wavelength' column in line_list.") if len(line_wavelengths) != len(line_list): raise ValueError("If line_wavelengths is specified, it must have the same " "length as line_pixels") From d573219b14b2285c064d6da9df3349570479e3fb Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Mon, 10 Apr 2023 16:18:32 -0400 Subject: [PATCH 30/31] Add QTable to catalog options --- specreduce/wavelength_calibration.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 19646ef7..2954ffae 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -109,14 +109,21 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non raise NotImplementedError("No catalogs are available yet, please input " "wavelengths with line_wavelengths or as a " "column in line_list") - if isinstance(catalog, list): + + if isinstance(catalog, QTable): + if "wavelength" not in catalog.columns: + raise ValueError("Catalog table must have a 'wavelength' column.") self._catalog = catalog else: - self._catalog = [catalog] - for cat in self._catalog: - if isinstance(cat, str): - if cat not in self._available_catalogs: - raise ValueError(f"Line list '{cat}' is not an available catalog.") + # This will need to be updated to match up with Tim's catalog code + if isinstance(catalog, list): + self._catalog = catalog + else: + self._catalog = [catalog] + for cat in self._catalog: + if isinstance(cat, str): + if cat not in self._available_catalogs: + raise ValueError(f"Line list '{cat}' is not an available catalog.") # Get the potential lines from any specified catalogs to use in matching self._potential_wavelengths = concatenate_catalogs(self._catalog) From 82170a5de0771594b4d1b1dd89991d67f441177d Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Tue, 11 Apr 2023 13:25:09 -0400 Subject: [PATCH 31/31] Finish addressing review comments --- docs/wavelength_calibration.rst | 15 ++-- .../tests/test_wavelength_calibration.py | 19 ++--- specreduce/wavelength_calibration.py | 69 ++++++++++++------- 3 files changed, 64 insertions(+), 39 deletions(-) diff --git a/docs/wavelength_calibration.rst b/docs/wavelength_calibration.rst index 28a5caa1..f34d7dec 100644 --- a/docs/wavelength_calibration.rst +++ b/docs/wavelength_calibration.rst @@ -17,20 +17,20 @@ spectrum, you can currently use ``specreduce`` to: The `~specreduce.wavelength_calibration.WavelengthCalibration1D` class can be used to fit a dispersion model to a list of line positions and wavelengths. Future development will implement catalogs of known lamp spectra for use in matching observed lines. In the -example below, the line positions (``line_list``) have already been extracted from +example below, the line positions (``pixel_centers``) have already been extracted from ``lamp_spectrum``:: import astropy.units as u from specreduce import WavelengthCalibration1D - line_list = [10, 22, 31, 43] + pixel_centers = [10, 22, 31, 43] wavelengths = [5340, 5410, 5476, 5543]*u.AA - test_cal = WavelengthCalibration1D(lamp_spectrum, line_list, - line_wavelengths=wavelengths) + test_cal = WavelengthCalibration1D(lamp_spectrum, line_pixels=pixel_centers, + line_wavelengths=wavelengths) calibrated_spectrum = test_cal.apply_to_spectrum(science_spectrum) The example above uses the default model (`~astropy.modeling.functional_models.Linear1D`) to fit the input spectral lines, and then applies the calculated WCS solution to a second -spectrum (``science_spectrum``). Any other ``astropy`` model can be provided as the +spectrum (``science_spectrum``). Any other 1D ``astropy`` model can be provided as the input ``model`` parameter to the `~specreduce.wavelength_calibration.WavelengthCalibration1D`. In the above example, the model fit and WCS construction is all done as part of the ``apply_to_spectrum()`` call, but you could also access the `~gwcs.wcs.WCS` object itself @@ -43,10 +43,11 @@ or ``input_spectrum`` properties are updated, since these will alter the calcula fit. You can also provide the input pixel locations and wavelengths of the lines as an -`~astropy.table.QTable`, with columns ``pixel_center`` and ``wavelength``:: +`~astropy.table.QTable` with (at minimum) columns ``pixel_center`` and ``wavelength``, +using the ``matched_line_list`` input argument:: from astropy.table import QTable pixels = [10, 20, 30, 40]*u.pix wavelength = [5340, 5410, 5476, 5543]*u.AA line_list = QTable([pixels, wavelength], names=["pixel_center", "wavelength"]) - test_cal = WavelengthCalibration1D(lamp_spectrum, line_list) \ No newline at end of file + test_cal = WavelengthCalibration1D(lamp_spectrum, matched_line_list=line_list) \ No newline at end of file diff --git a/specreduce/tests/test_wavelength_calibration.py b/specreduce/tests/test_wavelength_calibration.py index 3caf7cc6..93048f67 100644 --- a/specreduce/tests/test_wavelength_calibration.py +++ b/specreduce/tests/test_wavelength_calibration.py @@ -13,7 +13,7 @@ def test_linear_from_list(spec1d): centers = [0, 10, 20, 30] w = [5000, 5100, 5198, 5305]*u.AA - test = WavelengthCalibration1D(spec1d, centers, line_wavelengths=w) + test = WavelengthCalibration1D(spec1d, line_pixels=centers, line_wavelengths=w) spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) @@ -24,14 +24,14 @@ def test_wavelength_from_table(spec1d): centers = [0, 10, 20, 30] w = [5000, 5100, 5198, 5305]*u.AA table = QTable([w], names=["wavelength"]) - WavelengthCalibration1D(spec1d, centers, line_wavelengths=table) + WavelengthCalibration1D(spec1d, line_pixels=centers, line_wavelengths=table) def test_linear_from_table(spec1d): centers = [0, 10, 20, 30] w = [5000, 5100, 5198, 5305]*u.AA table = QTable([centers, w], names=["pixel_center", "wavelength"]) - test = WavelengthCalibration1D(spec1d, table) + test = WavelengthCalibration1D(spec1d, matched_line_list=table) spec2 = test.apply_to_spectrum(spec1d) assert_quantity_allclose(spec2.spectral_axis[0], 4998.8*u.AA) @@ -44,7 +44,8 @@ def test_poly_from_table(spec1d): w = [5005, 5110, 5214, 5330, 5438]*u.AA table = QTable([centers, w], names=["pixel_center", "wavelength"]) - test = WavelengthCalibration1D(spec1d, table, model=Polynomial1D(2), fitter=LinearLSQFitter()) + test = WavelengthCalibration1D(spec1d, matched_line_list=table, + model=Polynomial1D(2), fitter=LinearLSQFitter()) test.apply_to_spectrum(spec1d) assert_allclose(test.model.parameters, [5.00477143e+03, 1.03457143e+01, 1.28571429e-02]) @@ -53,7 +54,7 @@ def test_poly_from_table(spec1d): def test_replace_spectrum(spec1d, spec1d_with_emission_line): centers = [0, 10, 20, 30]*u.pix w = [5000, 5100, 5198, 5305]*u.AA - test = WavelengthCalibration1D(spec1d, centers, line_wavelengths=w) + test = WavelengthCalibration1D(spec1d, line_pixels=centers, line_wavelengths=w) # Accessing this property causes fits the model and caches the resulting WCS test.wcs assert "wcs" in test.__dict__ @@ -69,16 +70,16 @@ def test_expected_errors(spec1d): table = QTable([centers, w], names=["pixel_center", "wavelength"]) with pytest.raises(ValueError, match="Cannot specify line_wavelengths separately"): - WavelengthCalibration1D(spec1d, table, line_wavelengths=w) + WavelengthCalibration1D(spec1d, matched_line_list=table, line_wavelengths=w) with pytest.raises(ValueError, match="must have the same length"): w2 = [5005, 5110, 5214, 5330, 5438, 5500]*u.AA - WavelengthCalibration1D(spec1d, centers, line_wavelengths=w2) + WavelengthCalibration1D(spec1d, line_pixels=centers, line_wavelengths=w2) with pytest.raises(ValueError, match="astropy.units.Quantity array or" " as an astropy.table.QTable"): w2 = [5005, 5110, 5214, 5330, 5438] - WavelengthCalibration1D(spec1d, centers, line_wavelengths=w2) + WavelengthCalibration1D(spec1d, line_pixels=centers, line_wavelengths=w2) with pytest.raises(ValueError, match="specify at least one"): - WavelengthCalibration1D(spec1d, centers) + WavelengthCalibration1D(spec1d, line_pixels=centers) diff --git a/specreduce/wavelength_calibration.py b/specreduce/wavelength_calibration.py index 2954ffae..d870f06e 100644 --- a/specreduce/wavelength_calibration.py +++ b/specreduce/wavelength_calibration.py @@ -29,17 +29,20 @@ def concatenate_catalogs(): class WavelengthCalibration1D(): - def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=None, - model=Linear1D(), fitter=None): + def __init__(self, input_spectrum, matched_line_list=None, line_pixels=None, + line_wavelengths=None, catalog=None, model=Linear1D(), fitter=None): """ input_spectrum: `~specutils.Spectrum1D` A one-dimensional Spectrum1D calibration spectrum from an arc lamp or similar. - line_list: list, array, `~astropy.table.QTable` + matched_line_list: `~astropy.table.QTable`, optional + An `~astropy.table.QTable` table with (minimally) columns named + "pixel_center" and "wavelength" with known corresponding line pixel centers + and wavelengths populated. + line_pixels: list, array, `~astropy.table.QTable`, optional List or array of line pixel locations to anchor the wavelength solution fit. Will be converted to an astropy table internally if a list or array was input. Can also be input as an `~astropy.table.QTable` table with (minimally) a column - named "pixel_center" and optionally a "wavelength" column with known line - wavelengths populated. + named "pixel_center". line_wavelengths: `~astropy.units.Quantity`, `~astropy.table.QTable`, optional `astropy.units.Quantity` array of line wavelength values corresponding to the line pixels defined in ``line_list``. Does not have to be in the same order] @@ -55,10 +58,13 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non The fitter to use in optimizing the model fit. Defaults to `~astropy.modeling.fitting.LinearLSQFitter` if the model to fit is linear or `~astropy.modeling.fitting.LMLSQFitter` if the model to fit is non-linear. + + Note that either ``matched_line_list`` or ``line_pixels`` must be specified, + and if ``matched_line_list`` is not input, at least one of ``line_wavelengths`` + or ``catalog`` must be specified. """ self._input_spectrum = input_spectrum self._model = model - self._line_list = line_list self._cached_properties = ['wcs',] self.fitter = fitter self._potential_wavelengths = None @@ -67,28 +73,45 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non # ToDo: Implement having line catalogs self._available_catalogs = get_available_catalogs() - if isinstance(line_list, (list, np.ndarray)): - self._line_list = QTable([line_list], names=["pixel_center"]) + # We use either line_pixels or matched_line_list to create self._matched_line_list, + # and check that various requirements are fulfilled by the input args. + if matched_line_list is not None: + pixel_arg = "matched_line_list" + if not isinstance(matched_line_list, QTable): + raise ValueError("matched_line_list must be an astropy.table.QTable.") + self._matched_line_list = matched_line_list + elif line_pixels is not None: + pixel_arg = "line_pixels" + if isinstance(line_pixels, (list, np.ndarray)): + self._matched_line_list = QTable([line_pixels], names=["pixel_center"]) + elif isinstance(line_pixels, QTable): + self._matched_line_list = line_pixels + else: + raise ValueError("Either matched_line_list or line_pixels must be specified.") + + if "pixel_center" not in self._matched_line_list.columns: + raise ValueError(f"{pixel_arg} must have a 'pixel_center' column.") - if self._line_list["pixel_center"].unit is None: - self._line_list["pixel_center"].unit = u.pix + if self._matched_line_list["pixel_center"].unit is None: + self._matched_line_list["pixel_center"].unit = u.pix # Make sure our pixel locations are sorted - self._line_list.sort("pixel_center") + self._matched_line_list.sort("pixel_center") if (line_wavelengths is None and catalog is None - and "wavelength" not in self._line_list.columns): + and "wavelength" not in self._matched_line_list.columns): raise ValueError("You must specify at least one of line_wavelengths, " - "catalog, or 'wavelength' column in line_list.") + "catalog, or 'wavelength' column in matched_line_list.") # Sanity checks on line_wavelengths value if line_wavelengths is not None: - if isinstance(line_list, QTable) and "wavelength" in line_list.columns: + if (isinstance(self._matched_line_list, QTable) and + "wavelength" in self._matched_line_list.columns): raise ValueError("Cannot specify line_wavelengths separately if there is" - " a 'wavelength' column in line_list.") - if len(line_wavelengths) != len(line_list): + " a 'wavelength' column in matched_line_list.") + if len(line_wavelengths) != len(self._matched_line_list): raise ValueError("If line_wavelengths is specified, it must have the same " - "length as line_pixels") + f"length as {pixel_arg}") if not isinstance(line_wavelengths, (u.Quantity, QTable)): raise ValueError("line_wavelengths must be specified as an astropy.units.Quantity" " array or as an astropy.table.QTable") @@ -98,17 +121,17 @@ def __init__(self, input_spectrum, line_list, line_wavelengths=None, catalog=Non line_wavelengths[::-1].sort() else: line_wavelengths.sort() - self._line_list["wavelength"] = line_wavelengths + self._matched_line_list["wavelength"] = line_wavelengths elif isinstance(line_wavelengths, QTable): line_wavelengths.sort("wavelength") - self._line_list = hstack([self._line_list, line_wavelengths]) + self._matched_line_list = hstack([self._matched_line_list, line_wavelengths]) # Parse desired catalogs of lines for matching. if catalog is not None: # For now we avoid going into the later logic and just throw an error raise NotImplementedError("No catalogs are available yet, please input " "wavelengths with line_wavelengths or as a " - "column in line_list") + f"column in {pixel_arg}") if isinstance(catalog, QTable): if "wavelength" not in catalog.columns: @@ -171,8 +194,8 @@ def model(self, new_model): @cached_property def wcs(self): # computes and returns WCS after fitting self.model to self.refined_pixels - x = self._line_list["pixel_center"] - y = self._line_list["wavelength"] + x = self._matched_line_list["pixel_center"] + y = self._matched_line_list["wavelength"] if self.fitter is None: # Flexible defaulting if self.fitter is None @@ -189,7 +212,7 @@ def wcs(self): # Build a GWCS pipeline from the fitted model pixel_frame = cf.CoordinateFrame(1, "SPECTRAL", [0,], axes_names=["x",], unit=[u.pix,]) spectral_frame = cf.SpectralFrame(axes_names=["wavelength",], - unit=[self._line_list["wavelength"].unit,]) + unit=[self._matched_line_list["wavelength"].unit,]) pipeline = [(pixel_frame, self.model), (spectral_frame, None)]