diff --git a/docs/spectral_regions.rst b/docs/spectral_regions.rst index dd98e30db..15cae1b20 100644 --- a/docs/spectral_regions.rst +++ b/docs/spectral_regions.rst @@ -291,8 +291,13 @@ finding functions: (7.690954773869347 um, 8.690954773869347 um) This can be fed into the ``exclude_regions`` argument of the `~specutils.fitting.fit_generic_continuum` or -`~specutils.fitting.fit_continuum` functions to avoid fitting regions that contain line features. Conversely, users can -also invert the spectral region +`~specutils.fitting.fit_continuum` functions to avoid fitting regions that contain line features. Note that, +by default, this uses pythonic slicing, i.e., spectral values greater than or equal to the lower bound and +less than the upper bound of the region will be excluded from the fit. For convenience in some cases, the +``exclude_region_upper_bounds`` keyword can be set to ``True`` to exlude spectral values less than *or equal* +to the upper bound instead. + +Conversely, users can also invert the spectral region .. code-block:: python diff --git a/specutils/fitting/continuum.py b/specutils/fitting/continuum.py index fc0457e11..31512108d 100644 --- a/specutils/fitting/continuum.py +++ b/specutils/fitting/continuum.py @@ -57,7 +57,8 @@ def fit_generic_continuum(spectrum, median_window=3, model=Chebyshev1D(3), def fit_continuum(spectrum, model=Chebyshev1D(3), fitter=LevMarLSQFitter(), - exclude_regions=None, window=None, weights=None): + exclude_regions=None, exclude_region_upper_bounds=False, + window=None, weights=None): """ Entry point for fitting using the `~astropy.modeling.fitting` machinery. @@ -74,6 +75,11 @@ def fit_continuum(spectrum, model=Chebyshev1D(3), fitter=LevMarLSQFitter(), exclude_regions : list of 2-tuples List of regions to exclude in the fitting. Passed through to the fitmodels routine. + exclude_region_upper_bounds : bool, optional + Set to True to override the default pythonic slicing on the excluded regions (inclusive + lower bound, exclusive upper bound) and remove spectral values less than or equal to + the upper bound of the region rather than strictly less than the upper bound. Passed + through to the fitmodels routine. window : tuple of wavelengths Start and end wavelengths used for fitting. weights : list (NOT IMPLEMENTED YET) @@ -93,7 +99,7 @@ def fit_continuum(spectrum, model=Chebyshev1D(3), fitter=LevMarLSQFitter(), # Fit the flux to the model continuum_spectrum = fit_lines(spectrum, model, fitter, exclude_regions, - weights, w) + exclude_region_upper_bounds, weights, w) return continuum_spectrum diff --git a/specutils/fitting/fitmodels.py b/specutils/fitting/fitmodels.py index 27d09a731..a44c58945 100644 --- a/specutils/fitting/fitmodels.py +++ b/specutils/fitting/fitmodels.py @@ -259,8 +259,8 @@ def _generate_line_list_table(spectrum, emission_inds, absorption_inds): def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties=True), - exclude_regions=None, weights=None, window=None, get_fit_info=False, - **kwargs): + exclude_regions=None, exclude_region_upper_bounds=False, weights=None, + window=None, get_fit_info=False, **kwargs): """ Fit the input models to the spectrum. The parameter values of the input models will be used as the initial conditions for the fit. @@ -275,8 +275,12 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties Fitter instance to be used when fitting model to spectrum. exclude_regions : list of `~specutils.SpectralRegion` List of regions to exclude in the fitting. + exclude_region_upper_bounds : bool, optional + Set to True to override the default pythonic slicing on the excluded regions (inclusive + lower bound, exclusive upper bound) and remove spectral values less than or equal to + the upper bound of the region rather than strictly less than the upper bound. weights : array-like or 'unc', optional - If 'unc', the unceratinties from the spectrum object are used to + If 'unc', the uncertainties from the spectrum object are used to to calculate the weights. If array-like, represents the weights to use in the fitting. Note that if a mask is present on the spectrum, it will be applied to the ``weights`` as it would be to the spectrum @@ -310,7 +314,8 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties # if exclude_regions is not None: - spectrum = excise_regions(spectrum, exclude_regions) + spectrum = excise_regions(spectrum, exclude_regions, + inclusive_upper=exclude_region_upper_bounds) # # Make the model a list if not already @@ -358,6 +363,7 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties model_guess, fitter=fitter, exclude_regions=exclude_regions, + exclude_region_upper_bounds=exclude_region_upper_bounds, weights=weights, window=model_window, get_fit_info=get_fit_info, @@ -377,8 +383,8 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties def _fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties=True), - exclude_regions=None, weights=None, window=None, get_fit_info=False, - ignore_units=False, **kwargs): + exclude_regions=None, exclude_region_upper_bounds=False, weights=None, + window=None, get_fit_info=False, ignore_units=False, **kwargs): """ Fit the input model (initial conditions) to the spectrum. Output will be the same model with the parameters set based on the fitting. @@ -390,7 +396,8 @@ def _fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertaintie # if exclude_regions is not None: - spectrum = excise_regions(spectrum, exclude_regions) + spectrum = excise_regions(spectrum, exclude_regions, + inclusive_upper=exclude_region_upper_bounds) if isinstance(weights, str): if weights == 'unc': diff --git a/specutils/manipulation/utils.py b/specutils/manipulation/utils.py index d8eab1b5d..b8420b67c 100644 --- a/specutils/manipulation/utils.py +++ b/specutils/manipulation/utils.py @@ -8,7 +8,7 @@ __all__ = ['excise_regions', 'linear_exciser', 'spectrum_from_model'] -def true_exciser(spectrum, region): +def true_exciser(spectrum, region, inclusive_upper=False): """ Basic spectral excise method where the array elements in the spectral region defined by the parameter ``region`` (a `~specutils.SpectralRegion`) @@ -26,6 +26,10 @@ def true_exciser(spectrum, region): region : `~specutils.SpectralRegion` The region of the spectrum to remove. + inclusive_upper : bool, optional + Set to True to override the default pythonic slicing on the regions (inclusive + lower bound, exclusive upper bound) to include the specified upper bound instead. + Returns ------- spectrum : `~specutils.Spectrum1D` @@ -43,7 +47,11 @@ def true_exciser(spectrum, region): for subregion in region: # Find the indices of the spectral_axis array corresponding to the subregion - region_mask = (spectral_axis >= region.lower) & (spectral_axis < region.upper) + if inclusive_upper: + region_mask = (spectral_axis >= subregion.lower) & (spectral_axis <= subregion.upper) + else: + region_mask = (spectral_axis >= subregion.lower) & (spectral_axis < subregion.upper) + temp_indices = np.nonzero(region_mask)[0] if excise_indices is None: excise_indices = temp_indices @@ -76,7 +84,7 @@ def true_exciser(spectrum, region): radial_velocity=spectrum.radial_velocity if not isinstance(new_spectral_axis, SpectralCoord) else None) -def linear_exciser(spectrum, region): +def linear_exciser(spectrum, region, inclusive_upper=False): """ Basic spectral excise method where the spectral region defined by the parameter ``region`` (a `~specutils.SpectralRegion`) will result @@ -125,7 +133,11 @@ def linear_exciser(spectrum, region): for subregion in region: # Find the indices of the spectral_axis array corresponding to the subregion - region_mask = (spectral_axis >= subregion.lower) & (spectral_axis < subregion.upper) + if inclusive_upper: + region_mask = (spectral_axis >= subregion.lower) & (spectral_axis <= subregion.upper) + else: + region_mask = (spectral_axis >= subregion.lower) & (spectral_axis < subregion.upper) + inclusive_indices = np.nonzero(region_mask)[0] # Now set the flux values for these indices to be a # linear range @@ -150,7 +162,7 @@ def linear_exciser(spectrum, region): radial_velocity=spectrum.radial_velocity if not isinstance(spectral_axis, SpectralCoord) else None) -def excise_regions(spectrum, regions, exciser=true_exciser): +def excise_regions(spectrum, regions, exciser=true_exciser, inclusive_upper=False): """ Method to remove or replace the flux in the defined regions of the spectrum depending on the function provided in the ``exciser`` argument. @@ -173,6 +185,10 @@ def excise_regions(spectrum, regions, exciser=true_exciser): methods could be defined and used by this routine. default: true_exciser + inclusive_upper : bool, optional + Set to True to override the default pythonic slicing on the regions (inclusive + lower bound, exclusive upper bound) to include the specified upper bound instead. + Returns ------- spectrum : `~specutils.Spectrum1D` @@ -190,12 +206,12 @@ def excise_regions(spectrum, regions, exciser=true_exciser): raise ValueError('The spectrum parameter must be Spectrum1D object.') for region in regions: - spectrum = excise_region(spectrum, region, exciser) + spectrum = excise_region(spectrum, region, exciser, inclusive_upper=inclusive_upper) return spectrum -def excise_region(spectrum, region, exciser=true_exciser): +def excise_region(spectrum, region, exciser=true_exciser, inclusive_upper=False): """ Method to remove or replace the flux in the defined regions of the spectrum depending on the function provided in the ``exciser`` argument. @@ -217,6 +233,10 @@ def excise_region(spectrum, region, exciser=true_exciser): methods could be defined and used by this routine. default: true_exciser + inclusive_upper : bool, optional + Set to True to override the default pythonic slicing on the regions (inclusive + lower bound, exclusive upper bound) to include the specified upper bound instead. + Returns ------- spectrum : `~specutils.Spectrum1D` @@ -240,7 +260,7 @@ def excise_region(spectrum, region, exciser=true_exciser): # Call the exciser method # - return exciser(spectrum, region) + return exciser(spectrum, region, inclusive_upper=inclusive_upper) def spectrum_from_model(model_input, spectrum): diff --git a/specutils/tests/test_manipulation.py b/specutils/tests/test_manipulation.py index dd20b7ea4..841490d43 100644 --- a/specutils/tests/test_manipulation.py +++ b/specutils/tests/test_manipulation.py @@ -11,7 +11,7 @@ def test_true_exciser(): np.random.seed(84) - spectral_axis = np.linspace(5000,5100,num=100)*u.AA + spectral_axis = np.linspace(5000, 5099, num=100)*u.AA flux = (np.random.randn(100) + 3) * u.Jy spec = Spectrum1D(flux=flux, spectral_axis=spectral_axis) region = SpectralRegion([(5005,5010), (5060,5065)]*u.AA) @@ -21,10 +21,14 @@ def test_true_exciser(): assert len(excised_spec.flux) == len(spec.flux)-10 assert np.isclose(excised_spec.flux.sum(), 243.2617*u.Jy, atol=0.001*u.Jy) + excised_spec = excise_regions(spec, region, inclusive_upper=True) + assert len(excised_spec.spectral_axis) == len(spec.spectral_axis)-12 + assert len(excised_spec.flux) == len(spec.flux)-12 + def test_linear_exciser(): np.random.seed(84) - spectral_axis = np.linspace(5000,5100,num=100)*u.AA + spectral_axis = np.linspace(5000,5099,num=100)*u.AA flux = (np.random.rand(100)*100) * u.Jy spec = Spectrum1D(flux=flux, spectral_axis = spectral_axis) region = SpectralRegion([(5020,5030)]*u.AA)