Skip to content

Commit

Permalink
rename fibermap -> exp_fibermap for clarity
Browse files Browse the repository at this point in the history
  • Loading branch information
Stephen Bailey authored and Stephen Bailey committed Aug 25, 2023
1 parent 1500dac commit 0fba222
Showing 1 changed file with 52 additions and 46 deletions.
98 changes: 52 additions & 46 deletions py/desispec/coaddition.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,19 +219,25 @@ def coadd_fibermap(fibermap, onetile=False):
log.error(msg)
raise ValueError(msg)

#- make a copy of input fibermap that we can modify with new columns
fibermap = Table(fibermap, copy=True)
#- make a copy of input fibermap that we can modify with new columns.
#- This will become the per-exposure fibermap for the EXP_FIBERMAP HDU
exp_fibermap = Table(fibermap, copy=True)

#- Remove the "fibermap" input from the current namespace so that we don't accidentally use it
#- NOTE: does not actually delete/modify the original input
del fibermap

#- Get TARGETIDs, preserving order in which they first appear
targets, ii = ordered_unique(fibermap["TARGETID"], return_index=True)
tfmap = fibermap[ii]
#- tfmap = "Target Fiber Map", i.e. one row per target instead of one row per exposure
targets, ii = ordered_unique(exp_fibermap["TARGETID"], return_index=True)
tfmap = exp_fibermap[ii]
assert np.all(targets == tfmap['TARGETID'])
ntarget = targets.size

#- New columns to fill in for whether exposure was used in coadd
fibermap['IN_COADD_B'] = np.zeros(len(fibermap), dtype=bool)
fibermap['IN_COADD_R'] = np.zeros(len(fibermap), dtype=bool)
fibermap['IN_COADD_Z'] = np.zeros(len(fibermap), dtype=bool)
exp_fibermap['IN_COADD_B'] = np.zeros(len(exp_fibermap), dtype=bool)
exp_fibermap['IN_COADD_R'] = np.zeros(len(exp_fibermap), dtype=bool)
exp_fibermap['IN_COADD_Z'] = np.zeros(len(exp_fibermap), dtype=bool)

#- initialize NUMEXP=-1 to check that they all got filled later
tfmap['COADD_NUMEXP'] = np.zeros(len(tfmap), dtype=np.int16) - 1
Expand All @@ -256,7 +262,7 @@ def coadd_fibermap(fibermap, onetile=False):

#- Add other MEAN/RMS/STD columns
for k in mean_cols:
if k in fibermap.colnames :
if k in exp_fibermap.colnames :
if k.endswith('_RA') or k.endswith('_DEC') or k=='MJD':
dtype = np.float64
else:
Expand All @@ -274,7 +280,7 @@ def coadd_fibermap(fibermap, onetile=False):
tfmap.remove_column(k)

#- FIBER_RA/DEC handled differently due to RA wraparound and cos(dec)
if 'FIBER_RA' in fibermap.colnames and 'FIBER_DEC' in fibermap.colnames:
if 'FIBER_RA' in exp_fibermap.colnames and 'FIBER_DEC' in exp_fibermap.colnames:
tfmap.add_column(Column(np.zeros(ntarget, dtype=np.float64)), name='MEAN_FIBER_RA')
tfmap.add_column(Column(np.zeros(ntarget, dtype=np.float32)), name='STD_FIBER_RA')
tfmap.add_column(Column(np.zeros(ntarget, dtype=np.float64)), name='MEAN_FIBER_DEC')
Expand All @@ -283,7 +289,7 @@ def coadd_fibermap(fibermap, onetile=False):
tfmap.remove_column('FIBER_DEC')

#- MIN_, MAX_MJD over exposures used in coadd
if 'MJD' in fibermap.colnames :
if 'MJD' in exp_fibermap.colnames :
dtype = np.float64
if not 'MIN_MJD' in tfmap.dtype.names :
xx = Column(np.arange(ntarget, dtype=dtype))
Expand All @@ -293,7 +299,7 @@ def coadd_fibermap(fibermap, onetile=False):
tfmap.add_column(xx,name='MAX_MJD')

#- FIRSTNIGHT, LASTNIGHT over all exposures (not just good_coadds)
if 'NIGHT' in fibermap.colnames :
if 'NIGHT' in exp_fibermap.colnames :
dtype = np.int32
if not 'FIRSTNIGHT' in tfmap.dtype.names :
xx = Column(np.arange(ntarget, dtype=dtype))
Expand All @@ -307,44 +313,44 @@ def coadd_fibermap(fibermap, onetile=False):
if not 'COADD_FIBERSTATUS' in tfmap.dtype.names :
raise KeyError("no COADD_FIBERSTATUS column in tfmap")

if 'FIBERSTATUS' in fibermap.dtype.names :
if 'FIBERSTATUS' in exp_fibermap.dtype.names :
fiberstatus_key='FIBERSTATUS'
elif 'COADD_FIBERSTATUS' in fibermap.dtype.names :
elif 'COADD_FIBERSTATUS' in exp_fibermap.dtype.names :
fiberstatus_key='COADD_FIBERSTATUS'
else :
raise KeyError("no FIBERSTATUS nor COADD_FIBERSTATUS column in fibermap")

for i,tid in enumerate(targets) :
jj = fibermap["TARGETID"]==tid
jj = exp_fibermap["TARGETID"]==tid

#- Only a subset of "good" FIBERSTATUS flags are included in the coadd
targ_fibstatuses = fibermap[fiberstatus_key][jj]
targ_fibstatuses = exp_fibermap[fiberstatus_key][jj]
in_coadd_b = use_for_coadd(targ_fibstatuses, 'b')
in_coadd_r = use_for_coadd(targ_fibstatuses, 'r')
in_coadd_z = use_for_coadd(targ_fibstatuses, 'z')
good_coadds = (in_coadd_b | in_coadd_r | in_coadd_z)
coadd_numexp = np.count_nonzero(good_coadds)
tfmap['COADD_NUMEXP'][i] = coadd_numexp
fibermap['IN_COADD_B'][jj] = in_coadd_b
fibermap['IN_COADD_R'][jj] = in_coadd_r
fibermap['IN_COADD_Z'][jj] = in_coadd_z
exp_fibermap['IN_COADD_B'][jj] = in_coadd_b
exp_fibermap['IN_COADD_R'][jj] = in_coadd_r
exp_fibermap['IN_COADD_Z'][jj] = in_coadd_z

# Check if there are some good coadds to compute aggregate quantities;
# Otherwise just use all the (bad) exposures; will still count NUM on good_coadds
if coadd_numexp>0:
compute_coadds = good_coadds
# coadded FIBERSTATUS = bitwise AND of input FIBERSTATUS
tfmap['COADD_FIBERSTATUS'][i] = np.bitwise_and.reduce(fibermap[fiberstatus_key][jj][good_coadds])
tfmap['COADD_FIBERSTATUS'][i] = np.bitwise_and.reduce(exp_fibermap[fiberstatus_key][jj][good_coadds])
else:
compute_coadds = ~good_coadds
# if all inputs were bad, COADD_FIBERSTATUS is OR of inputs instead of AND
tfmap['COADD_FIBERSTATUS'][i] = np.bitwise_or.reduce(fibermap[fiberstatus_key][jj])
tfmap['COADD_FIBERSTATUS'][i] = np.bitwise_or.reduce(exp_fibermap[fiberstatus_key][jj])

#- For FIBER_RA/DEC quantities, only average over good coordinates.
# There is a bug that some "missing" coordinates were set to FIBER_RA=FIBER_DEC=0
# (we are assuming there are not valid targets at exactly 0,0; only missing coords)
if 'FIBER_RA' in fibermap.colnames and 'FIBER_DEC' in fibermap.colnames:
good_coords = (fibermap['FIBER_RA'][jj]!=0)|(fibermap['FIBER_DEC'][jj]!=0)
if 'FIBER_RA' in exp_fibermap.colnames and 'FIBER_DEC' in exp_fibermap.colnames:
good_coords = (exp_fibermap['FIBER_RA'][jj]!=0)|(exp_fibermap['FIBER_DEC'][jj]!=0)

#- Check whether entries with good coordinates exist (if not use all coordinates)
if np.count_nonzero(good_coords)>0:
Expand All @@ -363,36 +369,36 @@ def coadd_fibermap(fibermap, onetile=False):

# Note: NIGHT and TILEID may not be present when coadding previously
# coadded spectra.
if 'NIGHT' in fibermap.colnames:
tfmap['COADD_NUMNIGHT'][i] = len(np.unique(fibermap['NIGHT'][jj][good_coadds]))
if 'TILEID' in fibermap.colnames:
tfmap['COADD_NUMTILE'][i] = len(np.unique(fibermap['TILEID'][jj][good_coadds]))
if 'EXPTIME' in fibermap.colnames :
tfmap['COADD_EXPTIME'][i] = np.sum(fibermap['EXPTIME'][jj][good_coadds])
if 'NIGHT' in exp_fibermap.colnames:
tfmap['COADD_NUMNIGHT'][i] = len(np.unique(exp_fibermap['NIGHT'][jj][good_coadds]))
if 'TILEID' in exp_fibermap.colnames:
tfmap['COADD_NUMTILE'][i] = len(np.unique(exp_fibermap['TILEID'][jj][good_coadds]))
if 'EXPTIME' in exp_fibermap.colnames :
tfmap['COADD_EXPTIME'][i] = np.sum(exp_fibermap['EXPTIME'][jj][good_coadds])

#SJ: The RA needs something for the 0-360 deg wrap around probably in two steps
# with a first step grouping values close to the wrap around and the second step
# using modulo (% 360)
for k in mean_cols:
if k in fibermap.colnames :
if k in exp_fibermap.colnames :
if k.endswith('_RA') or k.endswith('_DEC'):
vals=fibermap[k][jj][compute_coadds_coords]
vals=exp_fibermap[k][jj][compute_coadds_coords]
else:
vals=fibermap[k][jj][compute_coadds]
vals=exp_fibermap[k][jj][compute_coadds]

tfmap['MEAN_'+k][i] = np.mean(vals)

for k in rms_cols:
if k in fibermap.colnames :
vals=fibermap[k][jj][compute_coadds]
if k in exp_fibermap.colnames :
vals=exp_fibermap[k][jj][compute_coadds]
# RMS includes mean offset, not same as STD
tfmap['RMS_'+k][i] = np.sqrt(np.mean(vals**2)).astype(np.float32)

#SJ: Check STD_FIBER_MAP with +360 value; doesn't do the wrap-around correctly
#- STD of FIBER_RA, FIBER_DEC in arcsec, handling cos(dec) and RA wrap
if 'FIBER_RA' in fibermap.colnames and 'FIBER_DEC' in fibermap.colnames:
decs = fibermap['FIBER_DEC'][jj][compute_coadds_coords]
ras = fibermap['FIBER_RA'][jj][compute_coadds_coords]
if 'FIBER_RA' in exp_fibermap.colnames and 'FIBER_DEC' in exp_fibermap.colnames:
decs = exp_fibermap['FIBER_DEC'][jj][compute_coadds_coords]
ras = exp_fibermap['FIBER_RA'][jj][compute_coadds_coords]
mean_ra, std_ra, mean_dec, std_dec = calc_mean_std_ra_dec(ras, decs)
tfmap['MEAN_FIBER_RA'][i] = mean_ra
tfmap['STD_FIBER_RA'][i] = np.float32(std_ra)
Expand All @@ -401,20 +407,20 @@ def coadd_fibermap(fibermap, onetile=False):

#- future proofing possibility of other STD cols
for k in std_cols:
if k in fibermap.colnames:
vals=fibermap[k][jj][compute_coadds]
if k in exp_fibermap.colnames:
vals=exp_fibermap[k][jj][compute_coadds]
# STD removes mean offset, not same as RMS
tfmap['STD_'+k][i] = np.std(vals).astype(np.float32)

# MIN_, MAX_MJD over exposures used in the coadd
if 'MJD' in fibermap.colnames :
vals=fibermap['MJD'][jj][compute_coadds]
if 'MJD' in exp_fibermap.colnames :
vals=exp_fibermap['MJD'][jj][compute_coadds]
tfmap['MIN_MJD'][i] = np.min(vals)
tfmap['MAX_MJD'][i] = np.max(vals)

# FIRST, LASTNIGHT over all exposures (not just compute_coadds)
if 'NIGHT' in fibermap.colnames :
vals=fibermap['NIGHT'][jj]
if 'NIGHT' in exp_fibermap.colnames :
vals=exp_fibermap['NIGHT'][jj]
tfmap['FIRSTNIGHT'][i] = np.min(vals)
tfmap['LASTNIGHT'][i] = np.max(vals)

Expand All @@ -423,8 +429,8 @@ def coadd_fibermap(fibermap, onetile=False):
#- (Note 2: these columns are place-holder for possible future use)
for k in ['FIBER_RA_IVAR', 'FIBER_DEC_IVAR',
'DELTA_X_IVAR', 'DELTA_Y_IVAR'] :
if k in fibermap.colnames :
tfmap[k][i]=1./np.mean(1./fibermap[k][jj][compute_coadds])
if k in exp_fibermap.colnames :
tfmap[k][i]=1./np.mean(1./exp_fibermap[k][jj][compute_coadds])

#- Remove some columns that apply to individual exp but not coadds
#- (even coadds of the same tile)
Expand All @@ -444,9 +450,9 @@ def coadd_fibermap(fibermap, onetile=False):
tfmap.remove_column(k)

#- keep exposure-specific columns that are present in the input fibermap
ii = np.isin(fibermap_perexp_cols, fibermap.dtype.names)
ii = np.isin(fibermap_perexp_cols, exp_fibermap.dtype.names)
keepcols = tuple(np.array(fibermap_perexp_cols)[ii])
exp_fibermap = fibermap[keepcols]
exp_fibermap = exp_fibermap[keepcols]

return tfmap, exp_fibermap

Expand Down

0 comments on commit 0fba222

Please sign in to comment.