Skip to content

Commit

Permalink
support coadd cameras without resolution
Browse files Browse the repository at this point in the history
  • Loading branch information
Stephen Bailey authored and Stephen Bailey committed Nov 3, 2023
1 parent d34aa51 commit 4a8795e
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 12 deletions.
30 changes: 20 additions & 10 deletions py/desispec/coaddition.py
Original file line number Diff line number Diff line change
Expand Up @@ -656,10 +656,10 @@ def coadd_cameras(spectra, cosmics_nsig=0., onetile=False) :

# ndiag = max of all cameras
ndiag=0
for b in sbands :
ndiag=max(ndiag,spectra.resolution_data[b].shape[1])
log.debug("ndiag= {}".format(ndiag))

if spectra.resolution_data is not None:
for b in sbands :
ndiag=max(ndiag,spectra.resolution_data[b].shape[1])
log.debug("ndiag=%d", ndiag)

b = sbands[0]
flux=np.zeros((ntarget,nwave),dtype=spectra.flux[b].dtype)
Expand All @@ -671,15 +671,21 @@ def coadd_cameras(spectra, cosmics_nsig=0., onetile=False) :
ivar_unmasked=ivar
mask=None

rdata=np.zeros((ntarget,ndiag,nwave),dtype=spectra.resolution_data[b].dtype)
if spectra.resolution_data is not None:
rdata=np.zeros((ntarget,ndiag,nwave),dtype=spectra.resolution_data[b].dtype)
else:
rdata = None

for b in spectra.bands :
log.debug("coadding band '{}'".format(b))

# indices
windices = windict[b]

band_ndiag = spectra.resolution_data[b].shape[1]
if spectra.resolution_data is not None:
band_ndiag = spectra.resolution_data[b].shape[1]
else:
band_ndiag = 0

if 'FIBERSTATUS' in spectra.fibermap.dtype.names:
fiberstatus = spectra.fibermap['FIBERSTATUS']
Expand Down Expand Up @@ -749,8 +755,9 @@ def coadd_cameras(spectra, cosmics_nsig=0., onetile=False) :

ivar[i,windices] += np.sum(ivarjj,axis=0)
flux[i,windices] += np.sum(ivarjj*spectra.flux[b][jj],axis=0)
for r in range(band_ndiag) :
rdata[i,r+(ndiag-band_ndiag)//2,windices] += np.sum((spectra.ivar[b][jj]*spectra.resolution_data[b][jj,r]),axis=0)
if spectra.resolution_data is not None:
for r in range(band_ndiag) :
rdata[i,r+(ndiag-band_ndiag)//2,windices] += np.sum((spectra.ivar[b][jj]*spectra.resolution_data[b][jj,r]),axis=0)
if spectra.mask is not None :
# this deserves some attention ...

Expand All @@ -770,7 +777,7 @@ def coadd_cameras(spectra, cosmics_nsig=0., onetile=False) :
if np.sum(ok)>0 :
flux[i][ok] /= ivar[i][ok]
ok=(ivar_unmasked[i]>0)
if np.sum(ok)>0 :
if np.sum(ok)>0 and rdata is not None:
rdata[i][:,ok] /= ivar_unmasked[i][ok]

if 'COADD_NUMEXP' in spectra.fibermap.colnames:
Expand All @@ -788,13 +795,16 @@ def coadd_cameras(spectra, cosmics_nsig=0., onetile=False) :
else :
dmask=None

if rdata is not None:
rdata = {bands:rdata}

res = Spectra(
bands=[bands,],
wave={bands:wave,},
flux={bands:flux,},
ivar={bands:ivar,},
mask=dmask,
resolution_data={bands:rdata,},
resolution_data=rdata,
fibermap=fibermap,
exp_fibermap=exp_fibermap,
meta=spectra.meta,
Expand Down
14 changes: 12 additions & 2 deletions py/desispec/specscore.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,14 @@ def compute_coadd_scores(coadd, specscores=None, update_coadd=True):

if coadd.bands == ['brz']:
#- i.e. this is a coadd across cameras
if coadd.resolution_data is not None:
rdat = coadd.resolution_data['brz']
else:
rdat = None

fr = Frame(coadd.wave['brz'], coadd.flux['brz'], coadd.ivar['brz'],
fibermap=coadd.fibermap, fibers=fibers, meta=coadd.meta,
resolution_data=coadd.resolution_data['brz'])
resolution_data=rdat)
for band in ['b', 'r', 'z']:
bandscores, bandcomments = compute_frame_scores(fr, band=band,
suffix='COADD', flux_per_angstrom=True)
Expand All @@ -63,9 +68,14 @@ def compute_coadd_scores(coadd, specscores=None, update_coadd=True):
#- otherwise try individual bands, upper or lowercase
for band in ['b', 'r', 'z', 'B', 'R', 'Z']:
if band in coadd.bands:
if coadd.resolution_data is not None:
rdat = coadd.resolution_data[band]
else:
rdat = None

fr = Frame(coadd.wave[band], coadd.flux[band], coadd.ivar[band],
fibermap=coadd.fibermap, fibers=fibers, meta=coadd.meta,
resolution_data=coadd.resolution_data[band])
resolution_data=rdat)
bandscores, bandcomments = compute_frame_scores(fr, band=band,
suffix='COADD', flux_per_angstrom=True)
scores.update(bandscores)
Expand Down
7 changes: 7 additions & 0 deletions py/desispec/test/test_coadd.py
Original file line number Diff line number Diff line change
Expand Up @@ -885,6 +885,13 @@ def test_coadd_cameras(self):
coadds = coadd_cameras(self.spectra)
self.assertEqual(len(coadds.wave['brz']), 7781)

# Test coadding without resolution or mask data
spec = self.spectra[:] #- copy before modifying
spec.resolution_data = None
spec.R = None
spec.mask = None
coadd = coadd_cameras(spec)


def test_suite():
"""Allows testing of only this module with the command::
Expand Down

0 comments on commit 4a8795e

Please sign in to comment.