astropy / specreduce

Tools for the reduction of spectroscopic observations from Optical and NIR instruments
https://specreduce.readthedocs.io
62 stars 38 forks source link

Background and Extract not respecting disp_axis correctly #223

Open jradavenport opened 4 months ago

jradavenport commented 4 months ago

When trying to use the non-default disp_axis and crossdisp_axis, @fchabourbarra and I are getting a variety of broadcasting errors, and results that are nonsense. Our workaround right now is transposing our data, which isn't ideal. The biggest issue is the case where it apparently extracts correctly (i.e. without complaint), but extracts along the wrong axis.

I've attached a (lengthy) minimum working example that tests I think most of the failure modes

import numpy as np
from astropy.nddata import CCDData
import astropy.units as u
from specreduce.tracing import ArrayTrace
from specreduce.background import Background
from specreduce.extract import BoxcarExtract 

# Making some fake data that is "vertical" (swapping disp/cross axes from default)
data = np.zeros((300,100)) # disp_axis=0, crossdisp_axis=1
data[[120,239],:] = 1 # put 2 lines at fixed wavelengths
img = CCDData(data, unit=u.adu)

print('img:', img.shape)

trace_array = np.ones(300)*75 # a line at crossdisp_axis = 75, has length=300
tr_bad = ArrayTrace(img, trace_array)
print('bad trace:', len(tr_bad.trace)) # should = 300

# if we transpose the data, ArrayTrace does it right
tr = ArrayTrace(img.data.T, trace_array)
print('trace:', len(tr.trace)) # should = 300

# this raises the error: 
#   ValueError: could not broadcast input array from shape (100,) into shape (300,)
# bk = Background(img, tr, disp_axis=0, crossdisp_axis=1)

# this also raises the error:
#   ValueError: could not broadcast input array from shape (100,) into shape (300,)
# bk = Background(img.data.T, tr, disp_axis=0, crossdisp_axis=1)

# setting axes back to default (which is WRONG for image, and doesn't match trace)
# this... works, but gives wrong length (100)
bk = Background(img, tr, disp_axis=1, crossdisp_axis=0)
print('bad bkgd_spectrum:',bk.bkg_spectrum().spectral_axis.shape)

# transposing the image fixes it
bk = Background(img.data.T, tr, disp_axis=1, crossdisp_axis=0)
print('bkgd_spectrum:',bk.bkg_spectrum().spectral_axis.shape)

# this appears to work
extract = BoxcarExtract(img, tr, disp_axis=0, crossdisp_axis=1)
# but this raises error:
#  ValueError: could not broadcast input array from shape (100,) into shape (300,)
# print(extract.spectrum.flux.shape)

# this appears to work, with WRONG axes that don't match trace
extract = BoxcarExtract(img, tr, disp_axis=1, crossdisp_axis=0)
# and gives wrong length (100)
print('bad extract:', extract.spectrum.flux.shape)

# transposing works
extract = BoxcarExtract(img.data.T, tr, disp_axis=1, crossdisp_axis=0)
# and gives wrong length (100)
print('extract:', extract.spectrum.flux.shape)

which outputs this:

img: (300, 100)
bad trace: 100
trace: 300
bad bkgd_spectrum: (100,)
bkgd_spectrum: (300,)
bad extract: (100,)
extract: (300,)
cshanahan1 commented 3 months ago

Looking into this, will update with a solution as soon as possible. Thanks!