Issue with specutils correlation #1283
Description
The section on Template cross-correlation in the docs (https://specutils.readthedocs.io/en/stable/analysis.html) returns the following error. I used Spectrum1D instead of Spectrum. I'm using specutils version 2.2.0 and astropy version 7.1.1. Error occurred in a Jupyter notebook.
UnitTypeError Traceback (most recent call last)
Cell In[26], line 31
27 ospec = Spectrum1D(spectral_axis=spec_axis, flux=flux1, uncertainty=uncertainty, velocity_convention='optical', rest_value=rest_value)
29 tspec = Spectrum1D(spectral_axis=spec_axis, flux=flux2, uncertainty=uncertainty)
---> 31 corr, lag = correlation.template_correlate(ospec, tspec)
File ~\anaconda3\envs\astropy\lib\site-packages\specutils\analysis\correlation.py:72, in template_correlate(observed_spectrum, template_spectrum, lag_units, apodization_window, resample)
69 log_template = template_spectrum
71 # apodize (might be a no-op if apodization_window is None)
---> 72 observed_log_spectrum, template_log_spectrum = _apodize(log_spectrum,
73 log_template,
74 apodization_window)
75 # Normalize template
76 normalization = _normalize(observed_log_spectrum, template_log_spectrum)
File ~\anaconda3\envs\astropy\lib\site-packages\specutils\analysis\correlation.py:118, in _apodize(spectrum, template, apodization_window)
116 def window(wlen):
117 return tukey(wlen, alpha=apodization_window)
--> 118 clean_spectrum = spectrum * window(len(spectrum.spectral_axis))
119 clean_template = template * window(len(template.spectral_axis))
121 return clean_spectrum, clean_template
File ~\anaconda3\envs\astropy\lib\site-packages\specutils\spectra\spectrum1d.py:510, in Spectrum1D.mul(self, other)
507 if not isinstance(other, NDDataRef):
508 other = u.Quantity(other)
--> 510 return self.multiply(other)
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\nddata\mixins\ndarithmetic.py:621, in NDArithmeticMixin.multiply(self, operand, operand2, **kwargs)
618 @sharedmethod
619 @format_doc(_arit_doc, name="multiplication", op="*")
620 def multiply(self, operand, operand2=None, **kwargs):
--> 621 return self._prepare_then_do_arithmetic(
622 np.multiply, operand, operand2, **kwargs
623 )
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\nddata\mixins\ndarithmetic.py:731, in NDArithmeticMixin._prepare_then_do_arithmetic(self_or_cls, operation, operand, operand2, **kwargs)
723 # At this point operand, operand2, kwargs and cls are determined.
724 if operand2 is not None and not issubclass(
725 operand2.class, NDArithmeticMixin
726 ):
(...)
729 # arrays, astropy quantities, masked quantities and of other subclasses
730 # of NDData.
--> 731 operand2 = cls(operand2)
733 # Now call the _arithmetics method to do the arithmetic.
734 result, init_kwds = operand._arithmetic(operation, operand2, **kwargs)
File ~\anaconda3\envs\astropy\lib\site-packages\specutils\spectra\spectrum1d.py:252, in Spectrum1D.init(self, flux, spectral_axis, wcs, velocity_convention, rest_value, redshift, radial_velocity, bin_specification, **kwargs)
247 # If no spectral_axis was provided, create a SpectralCoord based on
248 # the WCS
249 if spectral_axis is None:
250 # If spectral_axis wasn't provided, set _spectral_axis based on
251 # the WCS
--> 252 spec_axis = self.wcs.pixel_to_world(np.arange(self.flux.shape[-1]))
254 if spec_axis.unit.is_equivalent(u.one):
255 spec_axis = spec_axis * u.pixel
File ~\anaconda3\envs\astropy\lib\site-packages\specutils\utils\wcs_utils.py:217, in gwcs_from_array..SpectralGWCS.pixel_to_world(self, *args, **kwargs)
216 def pixel_to_world(self, *args, **kwargs):
--> 217 return super().pixel_to_world(*args, **kwargs).to(
218 orig_array.unit, equivalencies=u.spectral())
File ~\anaconda3\envs\astropy\lib\site-packages\gwcs\api.py:299, in GWCSAPIMixin.pixel_to_world(self, *pixel_arrays)
295 """
296 Convert pixel values to world coordinates.
297 """
298 pixels = self._sanitize_pixel_inputs(*pixel_arrays)
--> 299 return self(*pixels, with_units=True)
File ~\anaconda3\envs\astropy\lib\site-packages\gwcs\wcs.py:379, in WCS.call(self, *args, **kwargs)
377 if with_units:
378 if self.output_frame.naxes == 1:
--> 379 result = self.output_frame.coordinates(result)
380 else:
381 result = self.output_frame.coordinates(*result)
File ~\anaconda3\envs\astropy\lib\site-packages\gwcs\coordinate_frames.py:463, in SpectralFrame.coordinates(self, *args)
461 else:
462 if hasattr(args[0], 'unit'):
--> 463 return coord.SpectralCoord(*args).to(self.unit[0])
464 else:
465 return coord.SpectralCoord(*args, self.unit[0])
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\units\decorators.py:313, in QuantityInput.call..wrapper(*func_args, **func_kwargs)
311 # Call the original function with any equivalencies in force.
312 with add_enabled_equivalencies(self.equivalencies):
--> 313 return_ = wrapped_function(*func_args, **func_kwargs)
315 # Return
316 ra = wrapped_signature.return_annotation
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\coordinates\spectral_coordinate.py:193, in SpectralCoord.new(cls, value, unit, observer, target, radial_velocity, redshift, **kwargs)
182 @u.quantity_input(radial_velocity=u.km / u.s)
183 def new(
184 cls,
(...)
191 **kwargs,
192 ):
--> 193 obj = super().new(cls, value, unit=unit, **kwargs)
195 # There are two main modes of operation in this class. Either the
196 # observer and target are both defined, in which case the radial
197 # velocity and redshift are automatically computed from these, or
(...)
200 # observer are both specified, we can't also accept a radial velocity
201 # or redshift.
202 if target is not None and observer is not None:
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\coordinates\spectral_quantity.py:57, in SpectralQuantity.new(cls, value, unit, doppler_rest, doppler_convention, **kwargs)
54 def new(
55 cls, value, unit=None, doppler_rest=None, doppler_convention=None, **kwargs
56 ):
---> 57 obj = super().new(cls, value, unit=unit, **kwargs)
59 # If we're initializing from an existing SpectralQuantity, keep any
60 # parameters that aren't being overridden
61 if doppler_rest is None:
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\units\quantity.py:450, in Quantity.new(cls, value, unit, dtype, copy, order, subok, ndmin)
447 copy = False
449 if type(value) is not cls and not (subok and isinstance(value, cls)):
--> 450 value = value.view(cls)
452 if float_default and value.dtype.kind in "iu":
453 dtype = float
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\coordinates\spectral_coordinate.py:244, in SpectralCoord.array_finalize(self, obj)
243 def array_finalize(self, obj):
--> 244 super().array_finalize(obj)
245 self._radial_velocity = getattr(obj, "_radial_velocity", None)
246 self._observer = getattr(obj, "_observer", None)
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\coordinates\spectral_quantity.py:72, in SpectralQuantity.array_finalize(self, obj)
71 def array_finalize(self, obj):
---> 72 super().array_finalize(obj)
73 self._doppler_rest = getattr(obj, "_doppler_rest", None)
74 self._doppler_convention = getattr(obj, "_doppler_convention", None)
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\units\quantity.py:592, in Quantity.array_finalize(self, obj)
590 unit = getattr(obj, "_unit", None)
591 if unit is not None:
--> 592 self._set_unit(unit)
594 # Copy info if the original had info defined. Because of the way the
595 # DataInfo works, 'info' in obj.__dict__ is False until the
596 # info attribute is accessed or set.
597 if "info" in obj.dict:
File ~\anaconda3\envs\astropy\lib\site-packages\astropy\units\quantity.py:2121, in SpecificTypeQuantity._set_unit(self, unit)
2119 def _set_unit(self, unit):
2120 if unit is None or not unit.is_equivalent(self._equivalent_unit):
-> 2121 raise UnitTypeError(
2122 "{} instances require units equivalent to '{}'".format(
2123 type(self).name, self._equivalent_unit
2124 )
2125 + (
2126 ", but no unit was given."
2127 if unit is None
2128 else f", so cannot set it to '{unit}'."
2129 )
2130 )
2132 super()._set_unit(unit)
UnitTypeError: SpectralCoord instances require units equivalent to '(Unit("Hz"), Unit("m"), Unit("J"), Unit("1 / m"), Unit("km / s"))', so cannot set it to ''.