diff --git a/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb b/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb new file mode 100644 index 00000000..1ec8e072 --- /dev/null +++ b/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb @@ -0,0 +1,281 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Boxcar extraction with specreduce" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook demonstrates a simplified boxcar extraction scenario, using a rectified spectrum in a 2D spectral image where the dispersion axis is the second (X) axis. \n", + "\n", + "The extraction algorithm in `specreduce` is a plain adaptation of the algorithm presented in the MIRI LRS extraction notebook at https://github.com/spacetelescope/jdat_notebooks/tree/main/notebooks/MIRI_LRS_spectral_extraction This algorithm is demonstrated separately from the `specreduce` package, in the acompanying notebook http://localhost:8888/notebooks/notebook_sandbox/jwst_boxcar/jwst_boxcar_algorithm.ipynb\n", + "\n", + "Note that the extraction algorithm in the MIRI LRS extraction notebook performs simultaneous background subtraction when extracting from the source. Although it can only subtract the average of the two background extraction boxes. The original extraction code in `specreduce`'s `BoxcarExtract` class was capable of modeling the background with a polynomial along the cross-dispersion direction. This feature was lost in the conversion. \n", + "\n", + "Note also that we cannot yet perform offline, after-the-fact background subtractions from spectra extracted with `BoxcarExtract`. The reason is that `BoxcarExtract` delivers its `Spectrum1D` product with `spectral_axis` in units of `pixel`. Subsequent operations with such spectra, such as subtraction, are severely limited due to the `pixel` units." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "from matplotlib.colors import LogNorm\n", + "%matplotlib inline\n", + "\n", + "from pathlib import Path\n", + "from zipfile import ZipFile\n", + "\n", + "from astropy.io import fits\n", + "from astropy.table import Table\n", + "from astropy.visualization import simple_norm\n", + "from astropy.utils.data import download_file\n", + "\n", + "from jwst import datamodels\n", + "\n", + "from specreduce.extract import BoxcarExtract\n", + "from specreduce.tracing import FlatTrace\n", + "\n", + "import os\n", + "import tempfile\n", + "\n", + "import numpy as np\n", + "import ccdproc\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib as mpl" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Ingest s2d data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# data is taken from s2d file. x1d is used for comparison with pipeline extraction.\n", + "zipped_datapath = Path(download_file('https://stsci.box.com/shared/static/qdj79y7rv99wuui2hot0l3kg5ohq0ah9.zip', cache=True))\n", + "\n", + "data_dir = Path(tempfile.gettempdir())\n", + "\n", + "with ZipFile(zipped_datapath, 'r') as sample_data_zip:\n", + " sample_data_zip.extractall(data_dir)\n", + "\n", + "s2dfile = str(data_dir / \"nirspec_fssim_d1_s2d.fits\")\n", + "x1dfile = str(data_dir / \"nirspec_fssim_d1_x1d.fits\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# use a jwst datamodel to provide a good interface to the data and wcs info\n", + "s2d = datamodels.open(s2dfile)\n", + "image = s2d.slits[0].data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# display s2d image\n", + "norm_data = simple_norm(image, \"sqrt\")\n", + "plt.figure(figsize=(15, 15))\n", + "plt.imshow(image, norm=norm_data, origin=\"lower\")\n", + "plt.title(\"slit[0]\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# pipeline 1d extraction (for comparison)\n", + "jpipe_x1d = Table.read(x1dfile, hdu=1)\n", + "print(jpipe_x1d.columns)\n", + "# plot\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "ax.plot(jpipe_x1d['WAVELENGTH'], jpipe_x1d['FLUX'], 'k-', label=\"jpipe_x1d\")\n", + "ax.set_title(\"JWST Pipeline x1d extracted spectrum\")\n", + "ax.set_xlabel(\"wavelength\")\n", + "ax.set_ylabel(\"Flux Density [Jy]\")\n", + "ax.set_yscale(\"log\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define region to be extracted" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At this point, `specreduce` doesn't provide a tool for finding the trace. Since this is rectified data, we can use the `FlatTrace` class and find the spectrum position and width by eye." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# blow up of the region to be extracted\n", + "plt.figure(figsize=(15, 15))\n", + "plt.imshow(image[::,0:100], norm=norm_data, origin=\"lower\")\n", + "plt.title(\"slit[0] slice\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# extraction parameters based on image above\n", + "ext_center = 27\n", + "ext_width = 4" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot along cross-disperion cut showing the extraction parameters\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "y = np.arange(image.shape[0])\n", + "ax.plot(y, image[:,70], 'k-')\n", + "mm = np.array([ext_center, ext_center])\n", + "mm_y = ax.get_ylim()\n", + "\n", + "ax.plot(mm, mm_y, 'b--')\n", + "ax.plot(mm - ext_width/2., mm_y, 'g:')\n", + "ax.plot(mm + ext_width/2., mm_y, 'g:')\n", + "\n", + "ax.set_title(\"Cross-dispersion Cut at Pixel=70\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extract" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# define the Trace\n", + "trace = FlatTrace(image, ext_center)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# extract\n", + "boxcar = BoxcarExtract()\n", + "spectrum = boxcar(image, trace, width=ext_width)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot \n", + "f, ax = plt.subplots(figsize=(10, 6))\n", + "ax.plot(spectrum.flux.value, 'k-')\n", + "ax.set_title(\"Boxcar 1D extracted spectrum\")\n", + "ax.set_xlabel(r\"pixel\")\n", + "ax.set_ylabel(\"Flux Density [Jy]\")\n", + "ax.set_yscale(\"log\")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## About this notebook\n", + "\n", + "**Author:** Ivo Busko, JWST\n", + "**Updated On:** 2022-02-18" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Top of Page](#top)\n", + "\"Space " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/specreduce/extract.py b/specreduce/extract.py index 68113228..0fec112b 100644 --- a/specreduce/extract.py +++ b/specreduce/extract.py @@ -3,12 +3,11 @@ from dataclasses import dataclass import numpy as np -import matplotlib.pyplot as plt from astropy import units as u -from astropy.nddata import StdDevUncertainty from specreduce.core import SpecreduceOperation +from specreduce.tracing import FlatTrace from specutils import Spectrum1D __all__ = ['BoxcarExtract'] @@ -21,198 +20,126 @@ class BoxcarExtract(SpecreduceOperation): Parameters ---------- - img : nddata-compatible image - The input image - trace_object : - The trace of the spectrum to be extracted TODO: define - apwidth : int - The width of the extraction aperture in pixels - skysep : int - The spacing between the aperture and the sky regions - skywidth : int - The width of the sky regions in pixels - skydeg : int - The degree of the polynomial that's fit to the sky + image : nddata-compatible image + image with 2-D spectral image data + trace_object : Trace + trace object + width : float + width of extraction aperture in pixels + disp_axis : int + dispersion axis + crossdisp_axis : int + cross-dispersion axis Returns ------- spec : `~specutils.Spectrum1D` - The extracted spectrum - skyspec : `~specutils.Spectrum1D` - The sky spectrum used in the extraction process + The extracted 1d spectrum expressed in DN and pixel units """ - apwidth: int = 8 - skysep: int = 3 - skywidth: int = 7 - skydeg: int = 0 - - def __call__(self, img, trace_object): - self.last_trace = trace_object - self.last_img = img - - if self.apwidth < 1: - raise ValueError('apwidth must be >= 1') - if self.skysep < 1: - raise ValueError('skysep must be >= 1') - if self.skywidth < 1: - raise ValueError('skywidth must be >= 1') - - trace_line = trace_object.trace - - onedspec = np.zeros_like(trace_line) - skysubflux = np.zeros_like(trace_line) - fluxerr = np.zeros_like(trace_line) - mask = np.zeros_like(trace_line, dtype=bool) - - for i in range(0, len(trace_line)): - # if the trace isn't defined at a position (e.g. if it is out of the image boundaries), - # it will be masked. so we propagate that into the output mask and move on. - if np.ma.is_masked(trace_line[i]): - mask[i] = True - continue - - # first do the aperture flux - # juuuust in case the trace gets too close to an edge - widthup = self.apwidth / 2. - widthdn = self.apwidth / 2. - if (trace_line[i] + widthup > img.shape[0]): - widthup = img.shape[0] - trace_line[i] - 1. - if (trace_line[i] - widthdn < 0): - widthdn = trace_line[i] - 1. - - # extract from box around the trace line - low_end = trace_line[i] - widthdn - high_end = trace_line[i] + widthdn - - self._extract_from_box(img, i, low_end, high_end, onedspec) - - # now do the sky fit - # Note that we are not including fractional pixels, since we are doing - # a polynomial fit over the sky values. - j1 = self._find_nearest_int(trace_line[i] - self.apwidth/2. - - self.skysep - self.skywidth) - j2 = self._find_nearest_int(trace_line[i] - self.apwidth/2. - self.skysep) - sky_y_1 = np.arange(j1, j2) - - j1 = self._find_nearest_int(trace_line[i] + self.apwidth/2. + self.skysep) - j2 = self._find_nearest_int(trace_line[i] + self.apwidth/2. + - self.skysep + self.skywidth) - sky_y_2 = np.arange(j1, j2) - - sky_y = np.append(sky_y_1, sky_y_2) - - # sky can't be outside image - np_indices = np.indices(img[::, i].shape) - sky_y = np.intersect1d(sky_y, np_indices) - - sky_flux = img[sky_y, i] - if (self.skydeg > 0): - # fit a polynomial to the sky in this column - pfit = np.polyfit(sky_y, sky_flux, self.skydeg) - # define the aperture in this column - ap = np.arange( - self._find_nearest_int(trace_line[i] - self.apwidth/2.), - self._find_nearest_int(trace_line[i] + self.apwidth/2.) - ) - # evaluate the polynomial across the aperture, and sum - skysubflux[i] = np.nansum(np.polyval(pfit, ap)) - elif (self.skydeg == 0): - skysubflux[i] = np.nanmean(sky_flux) * self.apwidth - - # finally, compute the error in this pixel - sigma_bkg = np.nanstd(sky_flux) # stddev in the background data - n_bkg = np.float(len(sky_y)) # number of bkgd pixels - n_ap = self.apwidth # number of aperture pixels - - # based on aperture phot err description by F. Masci, Caltech: - # http://wise2.ipac.caltech.edu/staff/fmasci/ApPhotUncert.pdf - fluxerr[i] = np.sqrt( - np.nansum(onedspec[i] - skysubflux[i]) + (n_ap + n_ap**2 / n_bkg) * (sigma_bkg**2) - ) - - img_unit = u.DN - if hasattr(img, 'unit'): - img_unit = img.unit - - spec = Spectrum1D( - spectral_axis=np.arange(len(onedspec)) * u.pixel, - flux=onedspec * img_unit, - uncertainty=StdDevUncertainty(fluxerr), - mask=mask - ) - skyspec = Spectrum1D( - spectral_axis=np.arange(len(onedspec)) * u.pixel, - flux=skysubflux * img_unit, - mask=mask - ) - - return spec, skyspec - - def _extract_from_box(self, image, wave_index, low_end, high_end, extracted_result): - - # compute nearest integer endpoints defining an internal interval, - # and fractional pixel areas that remain outside this interval. - # (taken from the HST STIS pipeline code: - # https://github.com/spacetelescope/hstcal/blob/master/pkg/stis/calstis/cs6/x1dspec.c) - # - # This assumes that the pixel coordinates represent the center of the pixel. - # E.g. pixel at y=15.0 covers the image from y=14.5 to y=15.5 - - # nearest integer endpoints - j1 = self._find_nearest_int(low_end) - j2 = self._find_nearest_int(high_end) - - # fractional pixel areas at the end points - s1 = 0.5 - (low_end - j1) - s2 = 0.5 + high_end - j2 - - # add up the total flux around the trace_line - extracted_result[wave_index] = np.nansum(image[j1 + 1:j2, wave_index]) - extracted_result[wave_index] += np.nansum(image[j1, wave_index]) * s1 - extracted_result[wave_index] += np.nansum(image[j2, wave_index]) * s2 - - def _find_nearest_int(self, end_point): - if (end_point % 1) < 0.5: - return int(end_point) - else: - return int(end_point + 1) - - def get_checkplot(self): - trace_line = self.last_trace.line - - fig = plt.figure() - plt.imshow(self.last_img, origin='lower', aspect='auto', cmap=plt.cm.Greys_r) - plt.clim(np.percentile(self.last_img, (5, 98))) - - plt.plot(np.arange(len(trace_line)), trace_line, c='C0') - plt.fill_between( - np.arange(len(trace_line)), - trace_line + self.apwidth, - trace_line - self.apwidth, - color='C0', - alpha=0.5 - ) - plt.fill_between( - np.arange(len(trace_line)), - trace_line + self.apwidth + self.skysep, - trace_line + self.apwidth + self.skysep + self.skywidth, - color='C1', - alpha=0.5 - ) - plt.fill_between( - np.arange(len(trace_line)), - trace_line - self.apwidth - self.skysep, - trace_line - self.apwidth - self.skysep - self.skywidth, - color='C1', - alpha=0.5 - ) - plt.ylim( - np.min( - trace_line - (self.apwidth + self.skysep + self.skywidth) * 2 - ), - np.max( - trace_line + (self.apwidth + self.skysep + self.skywidth) * 2 - ) - ) - - return fig + + # TODO: should disp_axis and crossdisp_axis be defined in the Trace object? + + def __call__(self, image, trace_object, width=5, + disp_axis=1, crossdisp_axis=0): + """ + Extract the 1D spectrum using the boxcar method. + + Parameters + ---------- + image : nddata-compatible image + image with 2-D spectral image data + trace_object : Trace + trace object + width : float + width of extraction aperture in pixels + disp_axis : int + dispersion axis + crossdisp_axis : int + cross-dispersion axis + + + Returns + ------- + spec : `~specutils.Spectrum1D` + The extracted 1d spectrum with flux expressed in the same + units as the input image, or u.DN, and pixel units + """ + def _get_boxcar_weights(center, hwidth, npix): + """ + Compute weights given an aperture center, half width, + and number of pixels + """ + weights = np.zeros((npix)) + + # pixels with full weight + fullpixels = [max(0, int(center - hwidth + 1)), + min(int(center + hwidth), npix)] + weights[fullpixels[0]:fullpixels[1]] = 1.0 + + # pixels at the edges of the boxcar with partial weight + if fullpixels[0] > 0: + w = hwidth - (center - fullpixels[0] + 0.5) + if w >= 0: + weights[fullpixels[0] - 1] = w + else: + weights[fullpixels[0]] = 1. + w + if fullpixels[1] < npix: + weights[fullpixels[1]] = hwidth - (fullpixels[1] - center - 0.5) + + return weights + + def _ap_weight_image(trace, width, disp_axis, crossdisp_axis, image_shape): + + """ + Create a weight image that defines the desired extraction aperture. + + Parameters + ---------- + trace : Trace + trace object + width : float + width of extraction aperture in pixels + disp_axis : int + dispersion axis + crossdisp_axis : int + cross-dispersion axis + image_shape : tuple with 2 elements + size (shape) of image + + Returns + ------- + wimage : 2D image + weight image defining the aperture + """ + wimage = np.zeros(image_shape) + hwidth = 0.5 * width + image_sizes = image_shape[crossdisp_axis] + + # loop in dispersion direction and compute weights. + for i in range(image_shape[disp_axis]): + # TODO trace must handle transposed data (disp_axis == 0) + wimage[:, i] = _get_boxcar_weights(trace[i], hwidth, image_sizes) + + return wimage + + # TODO: this check can be removed if/when implemented as a check in FlatTrace + if isinstance(trace_object, FlatTrace): + if trace_object.trace_pos < 1: + raise ValueError('trace_object.trace_pos must be >= 1') + + # weight image to use for extraction + wimage = _ap_weight_image( + trace_object, + width, + disp_axis, + crossdisp_axis, + image.shape) + + # extract + ext1d = np.sum(image * wimage, axis=crossdisp_axis) + + # TODO: add wavelenght units, uncertainty and mask to spectrum1D object + spec = Spectrum1D(spectral_axis=np.arange(len(ext1d)) * u.pixel, + flux=ext1d * getattr(image, 'unit', u.DN)) + + return spec diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py index 531bd0d2..55eb5af6 100644 --- a/specreduce/tests/test_extract.py +++ b/specreduce/tests/test_extract.py @@ -4,7 +4,7 @@ from astropy.nddata import CCDData from specreduce.extract import BoxcarExtract -from specreduce.tracing import FlatTrace +from specreduce.tracing import FlatTrace, ArrayTrace # Test image is comprised of 30 rows with 10 columns each. Row content @@ -13,7 +13,7 @@ image = np.ones(shape=(30, 10)) for j in range(image.shape[0]): image[j, ::] *= j -image = CCDData(image, unit=u.count) +image = CCDData(image, unit=u.Jy) def test_extraction(): @@ -22,99 +22,58 @@ def test_extraction(): # extraction aperture sizes. # boxcar = BoxcarExtract() - - boxcar.apwidth = 5 - trace = FlatTrace(image, 15.0) - spectrum, bkg_spectrum = boxcar(image, trace) + + spectrum = boxcar(image, trace) assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 75.)) + assert spectrum.unit is not None and spectrum.unit == u.Jy trace.set_position(14.5) - spectrum, bkg_spectrum = boxcar(image, trace) + spectrum = boxcar(image, trace) assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 72.5)) trace.set_position(14.7) - spectrum, bkg_spectrum = boxcar(image, trace) + spectrum = boxcar(image, trace) assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 73.5)) - boxcar.apwidth = 6 - trace.set_position(15.0) - spectrum, bkg_spectrum = boxcar(image, trace) + spectrum = boxcar(image, trace, width=6) assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 90.)) trace.set_position(14.5) - spectrum, bkg_spectrum = boxcar(image, trace) + spectrum = boxcar(image, trace, width=6) assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 87.)) - boxcar.apwidth = 4.5 - trace.set_position(15.0) - spectrum, bkg_spectrum = boxcar(image, trace) + spectrum = boxcar(image, trace, width=4.5) assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.5)) - -def test_sky_extraction(): - # - # Try combinations of sky extraction parameters - # - boxcar = BoxcarExtract() - - boxcar.apwidth = 5. - boxcar.skysep = int(2) - boxcar.skywidth = 5. - - trace = FlatTrace(image, 15.0) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 75.)) - - trace.set_position(14.5) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 70.)) - - boxcar.skydeg = 1 - trace.set_position(15.0) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 75.)) + spectrum = boxcar(image, trace, width=4.7) + assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 70.5)) - trace.set_position(14.5) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 70.)) - - boxcar.skydeg = 2 - - trace.set_position(15.0) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 75.)) - - trace.set_position(14.5) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 70.)) - - boxcar.apwidth = 7. - boxcar.skysep = int(3) - boxcar.skywidth = 8. - - trace.set_position(15.0) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 105.)) - - trace.set_position(14.5) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 98.)) + trace.set_position(14.3) + spectrum = boxcar(image, trace, width=4.7) + assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.0)) def test_outside_image_condition(): # - # Trace is such that one of the sky regions lays partially outside the image + # Trace is such that extraction aperture lays partially outside the image # boxcar = BoxcarExtract() + trace = FlatTrace(image, 3.0) + + spectrum = boxcar(image, trace, width=10.) + assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 32.0)) + + +def test_array_trace(): + boxcar = BoxcarExtract() - boxcar.apwidth = 5. - boxcar.skysep = int(2) - boxcar.skywidth = 5. + trace_array = np.ones_like(image[1]) * 15. - trace = FlatTrace(image, 22.0) - spectrum, bkg_spectrum = boxcar(image, trace) - assert np.allclose(bkg_spectrum.flux.value, np.full_like(bkg_spectrum.flux.value, 99.375)) + trace = ArrayTrace(image, trace_array) + + spectrum = boxcar(image, trace) + assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 75.))