From c07d845fb0a5c09bbdfb7e8746de8acb2dff28b2 Mon Sep 17 00:00:00 2001 From: Stephen Bailey Date: Wed, 9 Aug 2023 13:38:33 -0700 Subject: [PATCH] coadd_fibermap handle RA wraparound --- py/desispec/coaddition.py | 69 +++++++++++++++++++++++++------- py/desispec/test/test_coadd.py | 72 +++++++++++++++++++++++++++++++++- 2 files changed, 125 insertions(+), 16 deletions(-) diff --git a/py/desispec/coaddition.py b/py/desispec/coaddition.py index 66e80b7fa..a73281ba4 100644 --- a/py/desispec/coaddition.py +++ b/py/desispec/coaddition.py @@ -130,6 +130,38 @@ 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, does not try to handle dec= +/-90 poles correctly + """ + 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 @@ -178,7 +210,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 @@ -189,8 +220,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': @@ -208,7 +240,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 @@ -266,7 +307,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) @@ -302,6 +343,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: @@ -312,19 +354,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) diff --git a/py/desispec/test/test_coadd.py b/py/desispec/test/test_coadd.py index 6fc44a32d..fc1041dd9 100644 --- a/py/desispec/test/test_coadd.py +++ b/py/desispec/test/test_coadd.py @@ -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 @@ -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 @@ -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