Source code for specutils.utils.wcs_utils

import copy

import numpy as np
from astropy import units as u
from astropy.modeling.models import Shift
from astropy.modeling.tabular import Tabular1D
from gwcs import WCS as GWCS
from gwcs import coordinate_frames as cf


class SpectralGWCS(GWCS):
    """
    This is a placeholder lookup-table GWCS created when a :class:`~specutils.Spectrum1D` is
    instantiated with a ``spectral_axis`` and no WCS.
    """
    def __init__(self, *args, **kwargs):
        self.original_unit = kwargs.pop("original_unit", "")
        super().__init__(*args, **kwargs)

    def copy(self):
        """
        Return a shallow copy of the object.

        Convenience method so user doesn't have to import the
        :mod:`copy` stdlib module.

        .. warning::
            Use `deepcopy` instead of `copy` unless you know why you need a
            shallow copy.
        """
        return copy.copy(self)

    def deepcopy(self):
        """
        Return a deep copy of the object.

        Convenience method so user doesn't have to import the
        :mod:`copy` stdlib module.
        """
        return copy.deepcopy(self)

    def pixel_to_world(self, *args, **kwargs):
        if self.original_unit == '':
            return u.Quantity(super().pixel_to_world_values(*args, **kwargs))
        return super().pixel_to_world(*args, **kwargs).to(
            self.original_unit, equivalencies=u.spectral())


[docs] def refraction_index(wavelength, method='Morton2000', co2=None): """ Calculates the index of refraction of dry air at standard temperature and pressure, at different wavelengths, using different methods. Parameters ---------- wavelength : `Quantity` object (number or sequence) Vacuum wavelengths with an astropy.unit. method : str, optional Method used to convert wavelengths. Options are: 'Morton2000' (default) - from Morton (2000, ApJS, 130, 403), eqn 8. Used by VALD, the Vienna Atomic Line Database. Very similar to Edlen (1966). 'Griesen2006' - from Greisen et al. (2006, A&A 446, 747), eqn. 65, standard used by International Union of Geodesy and Geophysics 'Edlen1953' - from Edlen (1953, J. Opt. Soc. Am, 43, 339). Standard adopted by IAU (resolution No. C15, Commission 44, XXI GA, 1991), which refers to Oosterhoff (1957) that uses Edlen (1953). Also used by Morton (1991, ApJS, 77, 119), which is frequently cited as IAU source. 'Edlen1966' - from Edlen (1966, Metrologia 2, 71), rederived constants from optical and near UV data. 'PeckReeder1972' - from Peck & Reeder (1972, J. Opt. Soc. 62), derived from additional infrared measurements (up to 1700 nm). 'Ciddor1996' - from Ciddor (1996, Appl. Opt. 35, 1566). Based on Peck & Reeder (1972), but updated to account for the changes in the international temperature scale and adjust the results for CO2 concentration. Arguably most accurate conversion available. Note that all options except for 'Griesen2006' have singularities in the far UV. 'Griesen2006' gives values that are slightly inconsistent with the other methods (~0.07 Angstrom difference at visible wavelengths), but it is the best option in the FUV due to the mathematical singularities in the others. See https://specutils.readthedocs.io/en/latest/wcs_utils.html for more detail, or https://github.com/astropy/specutils/issues/1162 for additional context. co2 : number, optional CO2 concentration in ppm. Only used for method='Ciddor1996'. If not given, a default concentration of 450 ppm is used. Returns ------- refr : number or sequence Index of refraction at each given air wavelength. """ VALID_METHODS = ['Griesen2006', 'Edlen1953', 'Edlen1966', 'Morton2000', 'PeckReeder1972', 'Ciddor1996'] assert isinstance(method, str), 'method must be a string' if method != 'Griesen2006' and wavelength.min() < 200 * u.nm: raise ValueError("The chosen method is invalid for wavelengths below 250 nm." " 'Griesen2006' is the only option for this wavelength range -" " see the specutils.utils.wcs_utils.refraction_index docstring" " for more detail.") method = method.lower() sigma2 = (1 / wavelength.to(u.um).value)**2 if method == 'griesen2006': refr = 1e-6 * (287.6155 + 1.62887 * sigma2 + 0.01360 * sigma2**2) elif method == 'edlen1953': refr = 6.4328e-5 + 2.94981e-2 / (146 - sigma2) + 2.5540e-4 / (41 - sigma2) elif method == 'edlen1966': refr = 8.34213e-5 + 2.406030e-2 / (130 - sigma2) + 1.5997e-4 / (38.9 - sigma2) elif method == 'morton2000': refr = 8.34254e-5 + 2.406147e-2 / (130 - sigma2) + 1.5998e-4 / (38.9 - sigma2) elif method == 'peckreeder1972': refr = 5.791817e-2 / (238.0185 - sigma2) + 1.67909e-3 / (57.362 - sigma2) elif method == 'ciddor1996': refr = 5.792105e-2 / (238.0185 - sigma2) + 1.67917e-3 / (57.362 - sigma2) if co2: refr *= 1 + 0.534e-6 * (co2 - 450) else: raise ValueError("Method must be one of " + ", ".join(VALID_METHODS)) return refr + 1
[docs] def vac_to_air(wavelength, method='Morton2000', co2=None): """ Converts vacuum to air wavelengths using different methods. Parameters ---------- wavelength : `Quantity` object (number or sequence) Vacuum wavelengths with an astropy.unit. method : str, optional One of the methods in refraction_index(), default is 'Morton2000'. co2 : number, optional Atmospheric CO2 concentration in ppm. Only used for method='Ciddor1996'. If not given, a default concentration of 450 ppm is used. Returns ------- air_wavelength : `Quantity` object (number or sequence) Air wavelengths with the same unit as wavelength. """ refr = refraction_index(wavelength, method=method, co2=co2) return wavelength / refr
[docs] def air_to_vac(wavelength, scheme='inversion', method='Morton2000', co2=None, precision=1e-12, maxiter=30): """ Converts air to vacuum wavelengths using different methods. Parameters ---------- wavelength : `Quantity` object (number or sequence) Air wavelengths with an astropy.unit. scheme : str, optional How to convert from vacuum to air wavelengths. Options are: * 'inversion' (default) - result is simply the inversion (1 / n) of the refraction index of air. Griesen et al. (2006) report that the error in naively inverting is less than 10^-9. * 'Piskunov' - uses an analytical solution derived by Nikolai Piskunov and used by the Vienna Atomic Line Database (VALD). * 'iteration' - uses an iterative scheme to invert the index of refraction. method : str, optional Only used if scheme is 'inversion' or 'iteration'. One of the methods in `~specutils.utils.wcs_utils.refraction_index`, default is 'Morton2000' co2 : number, optional Atmospheric CO2 concentration in ppm. Only used if scheme='inversion' and method='Ciddor1996'. If not given, a default concentration of 450 ppm is used. precision : float Maximum fractional value in refraction conversion beyond at which iteration will be stopped. Only used if scheme='iteration'. maxiter : integer Maximum number of iterations to run. Only used if scheme='iteration'. Returns ------- vac_wavelength : `Quantity` object (number or sequence) Vacuum wavelengths with the same unit as wavelength. """ VALID_SCHEMES = ['inversion', 'iteration', 'piskunov'] assert isinstance(scheme, str), 'scheme must be a string' scheme = scheme.lower() if scheme == 'inversion': refr = refraction_index(wavelength, method=method, co2=co2) elif scheme == 'piskunov': wlum = wavelength.to(u.angstrom).value sigma2 = (1e4 / wlum)**2 refr = (8.336624212083e-5 + 2.408926869968e-2 / (130.1065924522 - sigma2) + 1.599740894897e-4 / (38.92568793293 - sigma2)) + 1 elif scheme == 'iteration': # Refraction index is a function of vacuum wavelengths. # Iterate to get index of refraction that gives air wavelength that # is consistent with the reverse transformation. counter = 0 result = wavelength.copy() refr = refraction_index(wavelength, method=method, co2=co2) while True: counter += 1 diff = wavelength * refr - result if abs(diff.max().value) < precision: break if counter > maxiter: raise RuntimeError("Reached maximum number of iterations " "without reaching desired precision level.") result += diff refr = refraction_index(result, method=method, co2=co2) else: raise ValueError("Method must be one of " + ", ".join(VALID_SCHEMES)) return wavelength * refr
[docs] def air_to_vac_deriv(wavelength, method='Griesen2006'): """ Calculates the derivative d(wave_vacuum) / d(wave_air) using different methods. Parameters ---------- wavelength : `Quantity` object (number or sequence) Air wavelengths with an astropy.unit. method : str, optional Method used to convert wavelength derivative. Options are: 'Griesen2006' (default) - from Greisen et al. (2006, A&A 446, 747), eqn. 66. Returns ------- wave_deriv : `Quantity` object (number or sequence) Derivative d(wave_vacuum) / d(wave_air). """ assert method.lower() == 'griesen2006', "Only supported method is 'Griesen2006'" wlum = wavelength.to(u.um).value return (1 + 1e-6 * (287.6155 - 1.62887 / wlum**2 - 0.04080 / wlum**4))
[docs] def gwcs_from_array(array): """ Create a new WCS from provided tabular data. This defaults to being a GWCS object. """ orig_array = u.Quantity(array) coord_frame = cf.CoordinateFrame(naxes=1, axes_type=('SPECTRAL',), axes_order=(0,)) spec_frame = cf.SpectralFrame(unit=array.unit, axes_order=(0,)) # In order for the world_to_pixel transformation to automatically convert # input units, the equivalencies in the look up table have to be extended # with spectral unit information. SpectralTabular1D = type("SpectralTabular1D", (Tabular1D,), {'input_units_equivalencies': {'x0': u.spectral()}}) forward_transform = SpectralTabular1D(np.arange(len(array)), lookup_table=array) # If our spectral axis is in descending order, we have to flip the lookup # table to be ascending in order for world_to_pixel to work. if len(array) == 0 or array[-1] > array[0]: forward_transform.inverse = SpectralTabular1D( array, lookup_table=np.arange(len(array))) else: forward_transform.inverse = SpectralTabular1D( array[::-1], lookup_table=np.arange(len(array))[::-1]) tabular_gwcs = SpectralGWCS(original_unit = orig_array.unit, forward_transform=forward_transform, input_frame=coord_frame, output_frame=spec_frame) # Store the intended unit from the origin input array # tabular_gwcs._input_unit = orig_array.unit return tabular_gwcs
[docs] def gwcs_slice(self, item): """ This is a bit of a hack in order to fix the slicing of the WCS in the spectral dispersion direction. The NDData slices properly but the spectral dispersion result was not. There is code slightly downstream that sets the *number* of entries in the dispersion axis, this is just needed to shift to the correct starting element. When WCS gets the ability to do slicing then we might be able to remove this code. """ # Create shift of x-axis if isinstance(item, int): shift = item elif isinstance(item, slice): shift = item.start else: raise TypeError('Unknown index type {}, must be int or slice.'.format(item)) # Create copy as we need to modify this and return it. new_wcs = copy.deepcopy(self) if shift == 0: return new_wcs shifter = Shift(shift) # Get the current forward transform forward = new_wcs.forward_transform # Set the new transform new_wcs.set_transform(new_wcs.input_frame, new_wcs.output_frame, shifter | forward) return new_wcs