Add flux calibration

.github/workflows/muler-tests.yml

@@ -11,8 +11,8 @@ jobs:
         python-version: [3.9]
         specutils-version: [1.5, 1.9.1]
         astropy-version: [5.2]
-        numpy-version: [1.18, 1.24]
-
+        #numpy-version: [1.18, 1.24]
+        
     steps:
       - uses: actions/checkout@v2
       - uses: actions/checkout@v2

docs/requirements_actions.txt

@@ -16,3 +16,4 @@ bokeh # needed for gollum to install on RTD
 gollum==0.2.1
 numpydoc
 sphinx-gallery
+tynt

environment.yml

@@ -30,8 +30,11 @@ dependencies:
   - coverage
   - coveralls
   - twine
+  - dust_extinction
   - pip
   - pip:
       - sphinx-material
       - celerite2
       - specutils==1.9.1
+      - tynt
+

environment_M1.yml

@@ -30,7 +30,9 @@ dependencies:
   - coverage
   - coveralls
   - twine
+  - dust_extinction
   - pip
   - pip:
       - sphinx-material
       - celerite2
+      - tynt

src/muler/echelle.py

@@ -26,8 +26,10 @@
 from scipy.signal import savgol_filter
 from astropy.constants import R_jup, R_sun, G, M_jup, R_earth, c
 from astropy.modeling.physical_models import BlackBody
+from scipy.ndimage import median_filter, gaussian_filter1d
 import specutils
-from muler.utilities import apply_numpy_mask, is_list
+from muler.utilities import apply_numpy_mask, is_list, resample_list
+
 
 # from barycorrpy import get_BC_vel
 from astropy.coordinates import SkyCoord, EarthLocation
@@ -39,6 +41,9 @@
 import os
 import copy
 
+from specutils.manipulation import LinearInterpolatedResampler
+
+
 from specutils.spectra.spectral_region import SpectralRegion
 from specutils.analysis import equivalent_width
 
@@ -73,6 +78,7 @@ def __init__(self, *args, **kwargs):
         # self.ancillary_spectra = None
         super().__init__(*args, **kwargs)
 
+
     @property
     def snr(self):
         """The Signal-to-Noise Ratio :math:`\frac{S}{N}`, the flux divided by the uncertainty
@@ -699,9 +705,142 @@ def apply_boolean_mask(self, mask):
             )
 
         return spec
+
+    def get_slit_profile(self, lower=None, upper=None, slit_length=1.0):
+        """"For a 2D spectrum, returns the slit profile
+
+        Parameters
+        ----------
+        lower : AstroPy Quantity or float
+            The short wavelength limit at which to define the slit profile.
+            If the value is a float, it assume Angstrom units.
+        upper : AstroPy Quantity or float
+            The long wavelength limit at which to define the slit profiled.
+            If the value is a float, it assume Angstrom units.
+
+        Returns
+        -------
+        Array with the same height as the 2D spectrum of the median estimated slit profile
+        """
+        #Get the upper and lower wavelength limits in the correct units
+
+        assert len(np.shape(self.flux)) == 2, "Spectrum must be 2D to estimate slit profile." #Test to make sure this is a 2D spectrum
+
+        if lower is None:
+            lower = self.wavelength.min().value
+        if upper is None:
+            upper = self.wavelength.max().value
+
+        if type(lower) is not u.Quantity:
+            # Assume it's Angstroms
+            lower = lower * u.Angstrom
+        if type(upper) is not u.Quantity:
+            upper = upper * u.Angstrom
+
+        mask = (self.wavelength >= lower) & (self.wavelength <= upper)
+
+        flux = self.flux[:, mask].value
+        normalized_flux = flux / np.nansum(flux, axis=0)
+        median_slit_profile = np.nanmedian(normalized_flux, axis=1)
+
+        return median_slit_profile
+
+
+    def resample(self, target_spectrum):
+        """Resample spectrum onto a new spectral_axis.
+        Copied from gollum.
+
+        Parameters
+        ----------
+        target_spectrum : Spectrum1D
+            Spectrum whose wavelength grid you seek to match
+
+        Returns
+        -------
+        resampled_spec : PrecomputedSpectrum
+            Resampled spectrum
+        """
+        output = LinearInterpolatedResampler()(self, target_spectrum.wavelength)
+
+        return self._copy(
+            spectral_axis=output.wavelength.value * output.wavelength.unit,
+            flux=output.flux, uncertainty=output.uncertainty, meta=self.meta,
+            wcs=None,
+        )
+
+    def instrumental_broaden(self, resolving_power=55000):
+        r"""Instrumentally broaden the spectrum for a given instrumental resolution
+        Copied verbatim from gollum.
+
+        Known limitation: If the wavelength sampling changes with wavelength,
+          the convolution becomes inaccurate.  It may be better to FFT,
+          following Starfish.
+
+        Parameters
+        ----------
+        resolving_power : int
+            Instrumental resolving power :math:`R = \frac{\lambda}{\delta \lambda}`
+
+        Returns
+        -------
+        broadened_spec : PrecomputedSpectrum
+            Instrumentally broadened spectrum
+        """
+        # In detail the spectral resolution is wavelength dependent...
+        # For now we assume a constant resolving power
+        angstroms_per_pixel = np.median(np.diff(self.wavelength.angstrom))
+        lam0 = np.median(self.wavelength.value)
+        delta_lam = lam0 / resolving_power
+
+        scale_factor = 2.355
+        sigma = delta_lam / scale_factor / angstroms_per_pixel
+
+        convolved_flux = gaussian_filter1d(self.flux.value, sigma) * self.flux.unit
+        return self._copy(flux=convolved_flux)
+    def fill_nans(self, method=median_filter, **kwargs):
+        """Fill nans with the median of surrounding pixels using 
+        scipy.ndimage.median_filter
+
+        Parameters
+        ----------
+        method: def
+            def to apply to smooth surrounding pixels (e.g. scipy.ndimage.median_filter)
+        **kwargs:
+            Gets passed to method (e.g. size for scipy.ndimage.median_filter)
+        """
+        flux = self.flux
+        unc = self.uncertainty.array
+        filtered_flux = Quantity(method(flux.value, **kwargs), unit=self.flux.unit)
+        filtered_variance = method(unc**2, **kwargs)
+        filtered_unc = (filtered_variance**0.5)
+        found_nans = np.isnan(flux.value)
+        flux[found_nans] = filtered_flux[found_nans]
+        unc[found_nans] = filtered_unc[found_nans]
+
+        return self.__class__(
+            spectral_axis=self.spectral_axis, flux=flux, uncertainty=StdDevUncertainty(unc), meta=self.meta, wcs=None)
+    def apply(self, method=np.nansum, **kwargs):
+        """
+        Apply any method to the spectrum.  This is very general and can be used for many
+        things.  Uncertainty is propogated.
+
+        Parameters
+        ----------
+        method: def
+            def to apply to spectrum (e.g. np.nansum to collapse a multidimensional spectrum)
+        **kwargs:
+            Gets passed to method (e.g. axis for np.nansum)
+        """    
+        flux = self.flux
+        unc = self.uncertainty.array
+        flux = Quantity(method(self.flux.value, **kwargs), unit=self.flux.unit)
+        unc = method(self.uncertainty.array**2, **kwargs)**0.5
+        return self.__class__(
+            spectral_axis=self.spectral_axis, flux=flux, uncertainty=StdDevUncertainty(unc), meta=self.meta, wcs=None)
+
     def __pow__(self, power):
         """Take flux to a power while preserving the exiting flux units.
-        Uuseful for airmass correction.  Uncertainity is propogated by keeping the 
+        Uuseful for airmass correction.  Uncertainty is propogated by keeping the 
         singal-to-noise constant.
 
         Parameters
@@ -778,6 +917,34 @@ def trim_edges(self, limits=None):
 
         return spec_out
 
+    def trim_overlap(self, pivot=0.5):
+        """Trim all the edges that overlap with adjacent spectra (e.g. orders)
+        in the list.  Useful for running before stitch()."""
+        spec_out = copy.deepcopy(self)
+        n = len(spec_out)
+
+        for i in range(n): #Loop through each spectrum/order in list
+            #print('starting i ', i)
+            if i == 0: #Figure out where to trim the left side
+                left_limit = 0
+            elif self[i].spectral_axis[0] >  self[i-1].spectral_axis[-1]:
+                left_limit = 0
+            else:
+                mid_wave = self[i].spectral_axis[0]*(1-pivot) + self[i-1].spectral_axis[-1]*(pivot)
+                left_limit = np.where(self[i].spectral_axis > mid_wave)[-1][0] + 1
+            if i == n-1: #Figure out where to trim the right side
+                right_limit = len(self[i].spectral_axis)
+            elif self[i].spectral_axis[-1] <  self[i+1].spectral_axis[0]:
+                right_limit = len(self[i].spectral_axis)
+            else:
+                mid_wave = self[i].spectral_axis[-1]*(pivot) + self[i+1].spectral_axis[0]*(1-pivot)
+                right_limit = np.where(self[i].spectral_axis > mid_wave)[0][0] - 1
+
+            if left_limit > 0 or right_limit < len(self[i].spectral_axis):
+                spec_out[i] = spec_out[i].trim_edges((left_limit, right_limit))
+
+        return spec_out        
+
     def deblaze(self, method="spline"):
         """Remove blaze function from all orders by interpolating a spline function
 
@@ -949,7 +1116,7 @@ def __truediv__(self, other):
 
     def __pow__(self, power):
         """Take flux to a power while preserving the exiting flux units.
-        Uuseful for airmass correction.  Uncertainity is propagated by keeping the 
+        Uuseful for airmass correction.  Uncertainty is propagated by keeping the 
         singal-to-noise constant.
 
         Parameters
@@ -986,3 +1153,40 @@ def flatten(self, **kwargs):
             if "x_values" not in spec_out[i].meta:
                 spec_out[i].meta["x_values"] = self[i].meta["x_values"]
         return spec_out
+
+    def fill_nans(self, method=median_filter, **kwargs):
+        """Fill nans with the median of surrounding pixels using 
+        scipy.ndimage.median_filter
+
+        Parameters
+        ----------
+        method: def
+            def to apply to smooth surrounding pixels (e.g. scipy.ndimage.median_filter)
+        **kwargs:
+            Gets passed to method (e.g. size for scipy.ndimage.median_filter)
+        """
+        spec_out = copy.deepcopy(self)
+        for i in range(len(self)):
+            spec_out[i] = self[i].fill_nans(method=method, **kwargs)
+            if "x_values" not in spec_out[i].meta:
+                spec_out[i].meta["x_values"] = self[i].meta["x_values"]
+        return spec_out
+
+    def apply(self, method=np.nansum, **kwargs):
+        """
+        Apply any method to the spectral list.  This is very general and can be used for many
+        things.  Uncertainty is propogated.
+
+        Parameters
+        ----------
+        method: def
+            def to apply to spectrum (e.g. np.nansum to collapse a multidimensional spectrum)
+        **kwargs:
+            Gets passed to method (e.g. axis for np.nansum)
+        """    
+        spec_out = copy.deepcopy(self)
+        for i in range(len(self)):
+            spec_out[i] = self[i].apply(method=method, **kwargs)
+            if "x_values" not in spec_out[i].meta:
+                spec_out[i].meta["x_values"] = self[i].meta["x_values"]
+        return spec_out