Inconsistencies between different results of region extractions

[ivkram]

The rules of including/excluding boundaries of a spectral region during region extraction might be confusing, consider this example:

import astropy.units as u
import numpy as np
from specutils import Spectrum1D, SpectralRegion
from specutils.manipulation import extract_region

flux = np.random.randn(500) * u.Jy

spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.AA, flux=flux)
extract_region(spec1d, SpectralRegion(5000 * u.AA, 5100 * u.AA))  # [ 5000.0 Angstrom, ..., 5100.0 Angstrom ]

spec1d = Spectrum1D(spectral_axis=np.arange(0, 500) * u.pix, flux=flux)
extract_region(spec1d, SpectralRegion(0 * u.pix, 100 * u.pix))  # [ 0.0 pix, ..., 99.0 pix ]

One would expect that the spectral regions are treated as closed intervals in both cases. Consider another example where the upper bound is placed between the spectral axis values:

spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.AA, flux=flux)
extract_region(spec1d, SpectralRegion(5000 * u.AA, 5100.5 * u.AA))  # [ 5000.0 Angstrom, ..., 5100.0 Angstrom ]

spec1d = Spectrum1D(spectral_axis=np.arange(0, 500) * u.pix, flux=flux)
extract_region(spec1d, SpectralRegion(0 * u.pix, 100.5 * u.pix))  # [ 0.0 pix, ..., 100.0 pix ]

The resulting spectral axis is changed only in the "pixel" case. Final example:

spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.AA, flux=flux)
extract_region(spec1d, SpectralRegion(0 * u.pix, 5499 * u.pix))  # [ 5000.0 Angstrom, ..., 5499.0 Angstrom ]

spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.pix, flux=flux)
extract_region(spec1d, SpectralRegion(5000 * u.pix, 5499 * u.pix))  # [ 5000.0 pix, ..., 5498.0 pix ]

In the first case, 5499.0 Angstrom is included in the resulting spectral axis. In the second case, 5499.0 pix is excluded.

One possible solution would be to treat all spectral regions as closed intervals during region extraction.