Skip to content

Commit

Permalink
Merge pull request #2090 from desihub/rawrap
Browse files Browse the repository at this point in the history 
coadd_fibermap handle RA wraparound
  • Loading branch information
@sbailey
sbailey authored Aug 10, 2023
2 parents 93d3e34 + 4b1cf9f commit 98044b2
Show file tree
Hide file tree
Showing 2 changed files with 128 additions and 16 deletions.
72 changes: 58 additions & 14 deletions py/desispec/coaddition.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,41 @@
fibermap_fiberassign_cols + \
fibermap_coords_cols + fibermap_cframe_cols

def calc_mean_std_ra_dec(ras, decs):
"""
Calculate mean/std of ras, decs accounting for RA wraparound and cos(dec)
Args:
ras (array): input RA values in degrees
decs (array): input declination values in degrees
Returns: mean_ra, std_ra, mean_dec, std_dec
where the means are in degrees and the standard deviations are in arcsec,
including cos(dec) correction.
For efficiency, this does not try to handle dec= +/-90 poles correctly,
nor arbitrarily large spreads of angles. i.e. this is ok for a mean
of fiber positions scattered about a single target, but not for e.g.
a science analysis of the central location of a cluster of galaxies.
"""
ras = np.asarray(ras)
decs = np.asarray(decs)
if np.max(ras) - np.min(ras) > 180:
offset = 180.0
ras = (ras + offset) % 360
else:
offset = 0.0

mean_dec = np.mean(decs)
std_dec = np.std(decs) * 3600

mean_ra = (np.mean(ras) - offset + 360) % 360
std_ra = np.std(ras) * np.cos(np.radians(mean_dec)) * 3600

return mean_ra, std_ra, mean_dec, std_dec


def coadd_fibermap(fibermap, onetile=False):
"""
Coadds fibermap
Expand Down Expand Up @@ -178,7 +213,6 @@ def coadd_fibermap(fibermap, onetile=False):
# some cols get combined into mean or rms
mean_cols = [
'DELTA_X', 'DELTA_Y',
'FIBER_RA', 'FIBER_DEC',
'PSF_TO_FIBER_SPECFLUX', 'MJD'
]
# Note: treat the fiber coordinates separately because of missing coordinate problem
Expand All @@ -189,8 +223,9 @@ def coadd_fibermap(fibermap, onetile=False):

#- rms_cols and std_cols must also be in mean_cols
rms_cols = ['DELTA_X', 'DELTA_Y']
std_cols = ['FIBER_RA', 'FIBER_DEC']
std_cols = [] # currently none; RA/dec handled separately

#- Add other MEAN/RMS/STD columns
for k in mean_cols:
if k in fibermap.colnames :
if k.endswith('_RA') or k.endswith('_DEC') or k=='MJD':
Expand All @@ -208,7 +243,16 @@ def coadd_fibermap(fibermap, onetile=False):
tfmap.add_column(xx,name='STD_'+k)

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:
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')
tfmap.add_column(Column(np.zeros(ntarget, dtype=np.float32)), name='STD_FIBER_DEC')
tfmap.remove_column('FIBER_RA')
tfmap.remove_column('FIBER_DEC')

#- MIN_, MAX_MJD over exposures used in coadd
if 'MJD' in fibermap.colnames :
dtype = np.float64
Expand Down Expand Up @@ -266,7 +310,7 @@ def coadd_fibermap(fibermap, onetile=False):
#- 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:
if 'FIBER_RA' in fibermap.colnames and 'FIBER_DEC' in fibermap.colnames:
good_coords = (fibermap['FIBER_RA'][jj]!=0)|(fibermap['FIBER_DEC'][jj]!=0)

#- Check whether entries with good coordinates exist (if not use all coordinates)
Expand Down Expand Up @@ -302,6 +346,7 @@ def coadd_fibermap(fibermap, onetile=False):
vals=fibermap[k][jj][compute_coadds_coords]
else:
vals=fibermap[k][jj][compute_coadds]

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

for k in rms_cols:
Expand All @@ -312,19 +357,18 @@ def coadd_fibermap(fibermap, onetile=False):

#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:
dec = fibermap['FIBER_DEC'][jj][compute_coadds_coords][0]
vals = fibermap['FIBER_RA'][jj][compute_coadds_coords]
std = np.std(vals+360.0) * 3600 * np.cos(np.radians(dec))
tfmap['STD_FIBER_RA'][i] = std

if 'FIBER_DEC' in fibermap.colnames:
vals = fibermap['FIBER_DEC'][jj][compute_coadds_coords]
tfmap['STD_FIBER_DEC'][i] = np.std(vals) * 3600
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]
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)
tfmap['MEAN_FIBER_DEC'][i] = mean_dec
tfmap['STD_FIBER_DEC'][i] = np.float32(std_dec)

#- future proofing possibility of other STD cols
for k in std_cols:
if k in fibermap.colnames and k not in ('FIBER_RA', 'FIBER_DEC') :
if k in fibermap.colnames:
vals=fibermap[k][jj][compute_coadds]
# STD removes mean offset, not same as RMS
tfmap['STD_'+k][i] = np.std(vals).astype(np.float32)
Expand Down
72 changes: 70 additions & 2 deletions py/desispec/test/test_coadd.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,13 +163,18 @@ def test_coadd_fibermap_onetile(self):
fm['DESI_TARGET'] = [4,4,8,8,16,16]
fm['TILEID'] = [1,1,1,1,1,1]
fm['NIGHT'] = [20201220,20201221]*3
fm['MJD'] = 55555.0 + np.arange(6)
fm['EXPID'] = [10,20,11,21,12,22]
fm['FIBER'] = [5,6,]*3
fm['FIBERSTATUS'] = [0,0,0,0,0,0]
fm['FIBER_X'] = [1.1, 2.1]*3
fm['FIBER_Y'] = [10.2, 5.3]*3
fm['DELTA_X'] = [1.1, 2.1]*3
fm['DELTA_Y'] = [10.2, 5.3]*3
fm['FIBER_RA'] = [5.5, 6.6]*3
fm['FIBER_DEC'] = [7.7, 8.8]*3
fm['FLUX_R'] = np.ones(6)

fm['PSF_TO_FIBER_SPECFLUX'] = np.ones(6)
cofm, expfm = coadd_fibermap(fm, onetile=True)

#- Single tile coadds include these in the coadded fibermap
Expand All @@ -179,7 +184,9 @@ def test_coadd_fibermap_onetile(self):
self.assertIn(col, cofm.colnames)

#- but these columns should not be in the coadd
for col in ['NIGHT', 'EXPID', 'FIBERSTATUS', 'FIBER_X', 'FIBER_Y']:
for col in ['NIGHT', 'EXPID', 'FIBERSTATUS', 'FIBER_X', 'FIBER_Y',
'FIBER_RA', 'FIBER_DEC', 'DELTA_X', 'DELTA_Y',
'PSF_TO_FIBER_SPECFLUX']:
self.assertNotIn(col, cofm.colnames)

#- the exposure-level fibermap has columns specific to individual
Expand Down Expand Up @@ -354,6 +361,67 @@ def test_coadd_fibermap_badradec(self):
self.assertTrue(np.allclose(cofm['MEAN_FIBER_RA'], 0.0))
self.assertTrue(np.allclose(cofm['MEAN_FIBER_DEC'], 0.0))

def test_coadd_fibermap_ra_wrap(self):
"""Test coadding fibermap near RA=0 boundary"""
#- differences for (FIBER_RA - TARGET_RA)
delta_ra = np.array([-0.2, 0.0, 0.2])
ref_std_ra = np.float32(np.std(delta_ra) * 3600)
n = len(delta_ra)
dec = 60.0

#- one tile, 1 target
fm = Table()
fm['TARGETID'] = [111,] * n
fm['DESI_TARGET'] = [4,] * n
fm['TILEID'] = [1,] * n
fm['NIGHT'] = [20201220,] * n
fm['EXPID'] = 10 + np.arange(n)
fm['FIBER'] = [5,] * n
fm['FIBERSTATUS'] = [0,] * n
fm['TARGET_DEC'] = [dec,]*n
fm['FIBER_DEC'] = [dec,]*n

for ra in (359.9, 0, 0.1, 10):
fm['TARGET_RA'] = [ra,] * n
fm['FIBER_RA'] = ra + delta_ra / np.cos(np.radians(dec))

cofm, expfm = coadd_fibermap(fm, onetile=True)
self.assertAlmostEqual(cofm['MEAN_FIBER_RA'][0], ra, msg=f'mean(RA) at {ra=} {dec=}')
self.assertAlmostEqual(cofm['STD_FIBER_RA'][0], ref_std_ra, msg=f'std(RA) at {ra=} {dec=}')

def test_mean_std_ra_dec(self):
"""Test calc_mean_std_ra"""
from desispec.coaddition import calc_mean_std_ra_dec
delta = np.array([-0.2, 0.0, 0.2])
std_delta = np.std(delta) * 3600

for ra in (0.0, 0.1, 1.0, 179.9, 180.0, 180.1, 359.0, 359.9):
for dec in (-60, 0, 60):
ras = (ra + delta/np.cos(np.radians(dec)) + 360) % 360
decs = dec * np.ones(len(ras))
mean_ra, std_ra, mean_dec, std_dec = calc_mean_std_ra_dec(ras, decs)
self.assertAlmostEqual(mean_ra, ra,
msg=f'mean RA at {ra=} {dec=}')
self.assertAlmostEqual(std_ra, std_delta,
msg=f'std RA at {ra=} {dec=}')
self.assertAlmostEqual(mean_dec, dec,
msg=f'mean dec at {ra=} {dec=}')
self.assertAlmostEqual(std_dec, 0.0,
msg=f'std dec at {ra=} {dec=}')

#- Also check that std_dec doesn't depend upon RA
mean_ra, std_ra, mean_dec, std_dec = calc_mean_std_ra_dec(delta, delta)
self.assertAlmostEqual(std_dec, std_delta)
mean_ra, std_ra, mean_dec, std_dec = calc_mean_std_ra_dec(180+delta, delta)
self.assertAlmostEqual(std_dec, std_delta)

#- Confirm that 0 <= RA < 360
ras = [359.8, 0.1] # should average to 359.95, not -0.05
decs = [0, 0]
mean_ra, std_ra, mean_dec, std_dec = calc_mean_std_ra_dec(ras, decs)
self.assertAlmostEqual(mean_ra, 359.95) # not -0.05


def test_coadd_fibermap_mjd_night(self):
"""Test adding MIN/MAX/MEAN_MJD and FIRST/LASTNIGHT columns"""
nspec = 5
Expand Down

0 comments on commit 98044b2

Please sign in to comment.