thespacedoctor
commented
2 years ago We now have the 2D image maps giving wavelength, slit-position and order-number for each and every pixel. I have written the code to unpack these maps into a Pandas dataframe for easier (and faster) manipulation and association with data-pixel values.
We can now proceed to program the background subtraction algorithm.
I now understand the method and have a full workflow assembled. Here are my note (please excuse the non-rendered latex and citation) ...
Kelson Background Subtraction Method
The Kelson Background Subtraction method, originally outlined in Kelson, D. (2003)[#kelson2003], attempts to make full, or optimal, use of all data provided in a two-dimensional spectral image to accurately model a sky background image in the detector plane which can then be removed from the original data frame. The method allows for sky-subtraction to be performed in the early stages of the data-reduction process (just after flat-fielding) with no need to first identify and isolate cosmic-ray hits or other bad-pixel values.
The method requires an order locations table and a 2D dispersion map for the data, the output of the
soxs_spatial_solution, which can be used to transform the detector plane $(x,y)$ to a rectified wavelength-slit position plane $(\lambda,s)$.Method
The first step is to transform the detector pixel coordinates from $(x,y)$ to their rectified plane $(x_r,y_t)$ sub-pixel positions via the 2D dispersion map. This gives an array of $(x_r,y_t)$ positions and flux values.
Next we essentially flatten the $(x_r,y_t)$ array to provide an over-sampled wavelength-dependent $(x_r)$ flux array. Below (Figure 1 from Kelson, D. (2003)) we see a prominent night-skyline on the detector plane (left). In the right panel pixel values from a single row of the image stamp on the left are shown with the discrete pixel values forming a stepped profile. The thinner line is every single pixel value in the left-hand stamp plotted as a function of its wavelength-dependent position in the rectified coordinate space. This is the over-sampled wavelength-dependent $(x_r)$ flux array
This sampling of all pixel values is performed for each order (as seen below for an order sub-section in Figure 2 from Kelson, D. (2003)).
The narrow spikes seen in the bottom panel of the figure above represent the cosmic-ray hits and other bad-pixel values contaminating our measurements of the background night sky. A rolling-window median with robust 𝜎-clipping is used to remove these spikes. The smoothing routine also serendipitously removes any object-flux contaminating our measurement of the night sky (seen as a faint trace around row 65 in the top panel above). Results of the smoothing routine can be seen in the (artificially shifted) top line in the plot below from Figure 2 of Kelson, D. (2003).
Alongside the smoothed data, a rolling window 4𝜎 scatter array is compiled from the original data. This array is used later to weight the remaining data.
What essentially remains is the rectified plane ($x_r,y_t$) data with holes where deviant values have been clipped. A B-spline surface (bivariate B-spline) is now fitted to this remaining data, with data weighted by the inverse of the 4𝜎 scatter values calculated before.[^20211123031702]. Using the fitted B-spline surface, modelled values at each $(x_r,y_t)$ sub-pixel position are calculated, which are conveniently located at the exact $(x,y)$ pixel locations in the original detector plane data. We now have a detector-plane image of the modelled sky background.
In the three panel image below (Figure 4 from Kelson, D. 2003) we have original data values (top), the modelled B-spline surface of the sky-background (middle) and the sky-background subtracted image (bottom).
Workflow Diagram
[Not Cited][#kelson2003]
[#kelson2003]: Kelson, D. (2003). 'Optimal Techniques in Two‐dimensional Spectroscopy: Background Subtraction for the 21st Century.' Publications of the Astronomical Society of the Pacific, 115 (808)
[^20211123031702]: We may need to experiment with the location of knots in the B-spline surface to best fit the data ... lower density of knots makes for a smoother surface but a higher density will fit large gradients generated by bright night skylines. However, there may be no need to optimise knot locations.