Merge 7166f5f into 7a17923

bin/desi_zcatalog

@@ -401,6 +401,25 @@ def main():
     else:
         expfm = None
 
+    #- Add FIRSTNIGHT for tile-based cumulative catalogs
+    #- (LASTNIGHT was added while reading from NIGHT header keyword)
+    if args.group == 'cumulative' and expfm is not None and 'FIRSTNIGHT' not in zcat.colnames:
+        log.info('Adding FIRSTNIGHT per tile')
+        icol = zcat.colnames.index('LASTNIGHT')
+        zcat.add_column(np.zeros(len(zcat), dtype=np.int32),
+                    index=icol, name='FIRSTNIGHT')
+        for tilefm in Table(expfm[['TILEID', 'NIGHT']]).group_by('TILEID').groups:
+            tileid = tilefm['TILEID'][0]
+            iitile = zcat['TILEID'] == tileid
+            zcat['FIRSTNIGHT'][iitile] = np.min(tilefm['NIGHT'])
+
+        #- all FIRSTNIGHT entries should be filled (no more zeros)
+        bad = zcat['FIRSTNIGHT'] == 0
+        if np.any(bad):
+            badtiles = np.unique(zcat['TILEID'][bad])
+            raise ValueError(f'FIRSTNIGHT not set for tiles {badtiles}')
+
+
     #- if TARGETIDs appear more than once, which one is best within this catalog?
     if 'TSNR2_LRG' in zcat.colnames and 'ZWARN' in zcat.colnames:
         log.info('Finding best spectrum for each target')

py/desispec/coaddition.py

@@ -245,10 +245,11 @@ def coadd_fibermap(fibermap, onetile=False):
     tfmap['COADD_NUMNIGHT'] = np.zeros(len(tfmap), dtype=np.int16) - 1
     tfmap['COADD_NUMTILE'] = np.zeros(len(tfmap), dtype=np.int16) - 1
 
-    # some cols get combined into mean or rms
+    # some cols get combined into mean or rms;
+    # MJD handled separately to get min/max/mean
     mean_cols = [
         'DELTA_X', 'DELTA_Y',
-        'PSF_TO_FIBER_SPECFLUX', 'MJD'
+        'PSF_TO_FIBER_SPECFLUX',
         ]
     # Note: treat the fiber coordinates separately because of missing coordinate problem
     # that require an additional "good_coords" condition relative to other mean cols
@@ -263,10 +264,7 @@ def coadd_fibermap(fibermap, onetile=False):
     #- Add other MEAN/RMS/STD columns
     for k in mean_cols:
         if k in exp_fibermap.colnames :
-            if k.endswith('_RA') or k.endswith('_DEC') or k=='MJD':
-                dtype = np.float64
-            else:
-                dtype = np.float32
+            dtype = np.float32
             if k in mean_cols:
                 xx = Column(np.zeros(ntarget, dtype=dtype))
                 tfmap.add_column(xx,name='MEAN_'+k)
@@ -288,26 +286,19 @@ def coadd_fibermap(fibermap, onetile=False):
         tfmap.remove_column('FIBER_RA')
         tfmap.remove_column('FIBER_DEC')
 
-    #- MIN_, MAX_MJD over exposures used in coadd
+    #- MIN_, MAX_, MEAN_MJD over exposures used in coadd
     if 'MJD' in exp_fibermap.colnames :
         dtype = np.float64
         if not 'MIN_MJD' in tfmap.dtype.names :
-            xx = Column(np.arange(ntarget, dtype=dtype))
+            xx = Column(np.zeros(ntarget, dtype=dtype))
             tfmap.add_column(xx,name='MIN_MJD')
         if not 'MAX_MJD' in tfmap.dtype.names :
-            xx = Column(np.arange(ntarget, dtype=dtype))
+            xx = Column(np.zeros(ntarget, dtype=dtype))
             tfmap.add_column(xx,name='MAX_MJD')
+        if not 'MEAN_MJD' in tfmap.dtype.names :
+            xx = Column(np.zeros(ntarget, dtype=dtype))
+            tfmap.add_column(xx,name='MEAN_MJD')
 
-    #- FIRSTNIGHT, LASTNIGHT over all exposures (not just good_coadds)
-    if 'NIGHT' in exp_fibermap.colnames :
-        dtype = np.int32
-        if not 'FIRSTNIGHT' in tfmap.dtype.names :
-            xx = Column(np.arange(ntarget, dtype=dtype))
-            tfmap.add_column(xx,name='FIRSTNIGHT')
-        if not 'LASTNIGHT' in tfmap.dtype.names :
-            xx = Column(np.arange(ntarget, dtype=dtype))
-            tfmap.add_column(xx,name='LASTNIGHT')
-
     if 'FIBERSTATUS' in tfmap.dtype.names :
         tfmap.rename_column('FIBERSTATUS', 'COADD_FIBERSTATUS')
     if not  'COADD_FIBERSTATUS' in tfmap.dtype.names :
@@ -376,16 +367,11 @@ def coadd_fibermap(fibermap, onetile=False):
         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)
+        # Calc MEAN_* and RMS_* columns using the same exposures as the coadd
+        # Note RA/DEC/MJD are handled separately
         for k in mean_cols:
             if k in exp_fibermap.colnames :
-                if k.endswith('_RA') or k.endswith('_DEC'):
-                    vals=exp_fibermap[k][jj][compute_coadds_coords]
-                else:
-                    vals=exp_fibermap[k][jj][compute_coadds]
-
+                vals=exp_fibermap[k][jj][compute_coadds]
                 tfmap['MEAN_'+k][i] = np.mean(vals)
 
         for k in rms_cols:
@@ -394,7 +380,6 @@ def coadd_fibermap(fibermap, onetile=False):
                 # 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 exp_fibermap.colnames and 'FIBER_DEC' in exp_fibermap.colnames:
             decs = exp_fibermap['FIBER_DEC'][jj][compute_coadds_coords]
@@ -417,13 +402,8 @@ def coadd_fibermap(fibermap, onetile=False):
             vals=exp_fibermap['MJD'][jj][compute_coadds]
             tfmap['MIN_MJD'][i] = np.min(vals)
             tfmap['MAX_MJD'][i] = np.max(vals)
+            tfmap['MEAN_MJD'][i] = np.mean(vals)
 
-        # FIRST, LASTNIGHT over all exposures (not just compute_coadds)
-        if 'NIGHT' in exp_fibermap.colnames :
-            vals=exp_fibermap['NIGHT'][jj]
-            tfmap['FIRSTNIGHT'][i] = np.min(vals)
-            tfmap['LASTNIGHT'][i] = np.max(vals)
-
         # Error propagation of IVAR values when taking an unweighted MEAN 
         #- (Note 1: IVAR will be 0.0 if any of ivar[compute_coadds]=0)
         #- (Note 2: these columns are place-holder for possible future use)    

py/desispec/test/test_coadd.py

@@ -481,7 +481,7 @@ def test_mean_std_ra_dec(self):
 
 
     def test_coadd_fibermap_mjd_night(self):
-        """Test adding MIN/MAX/MEAN_MJD and FIRST/LASTNIGHT columns"""
+        """Test adding MIN/MAX/MEAN_MJD columns"""
         nspec = 5
         fm = Table()
         fm['TARGETID'] = 111 * np.ones(nspec, dtype=int)
@@ -494,8 +494,6 @@ def test_coadd_fibermap_mjd_night(self):
         #- with NIGHT but not MJD
         cofm, expfm = coadd_fibermap(fm, onetile=True)
         self.assertEqual(len(cofm), 1)  #- single target in this test
-        self.assertEqual(cofm['FIRSTNIGHT'][0], np.min(fm['NIGHT']))
-        self.assertEqual(cofm['LASTNIGHT'][0], np.max(fm['NIGHT']))
         self.assertNotIn('MIN_MJD', cofm.colnames)
 
         #- also with MJD
@@ -506,16 +504,13 @@ def test_coadd_fibermap_mjd_night(self):
         self.assertEqual(cofm['MEAN_MJD'][0], np.mean(fm['MJD']))
 
         #- with some fibers masked:
-        #- FIRSTNIGHT/LASTNIGHT still based on all exp,
-        #- while MIN/MAX/MEAN_MJD based only upon good ones
+        #- MIN/MAX/MEAN_MJD based only upon good ones
         fm['FIBERSTATUS'][0] = fibermask.BADFIBER     #- bad
         fm['FIBERSTATUS'][1] = fibermask.RESTRICTED   #- ok
         fm['FIBERSTATUS'][2] = fibermask.BADAMPR      #- ok for fibermap
         ok = np.ones(nspec, dtype=bool)
         ok[0] = False
         cofm, expfm = coadd_fibermap(fm, onetile=True)
-        self.assertEqual(cofm['FIRSTNIGHT'][0], np.min(fm['NIGHT']))
-        self.assertEqual(cofm['LASTNIGHT'][0], np.max(fm['NIGHT']))
         self.assertEqual(cofm['MIN_MJD'][0], np.min(fm['MJD'][ok]))
         self.assertEqual(cofm['MAX_MJD'][0], np.max(fm['MJD'][ok]))
         self.assertEqual(cofm['MEAN_MJD'][0], np.mean(fm['MJD'][ok]))
@@ -524,17 +519,22 @@ def test_coadd_fibermap_mjd_night(self):
         fm['TARGETID'][0:2] += 1
         fm['FIBERSTATUS'] = 0
         cofm, expfm = coadd_fibermap(fm, onetile=True)
-        self.assertEqual(cofm['FIRSTNIGHT'][0], np.min(fm['NIGHT'][0:2]))
-        self.assertEqual(cofm['LASTNIGHT'][0],  np.max(fm['NIGHT'][0:2]))
         self.assertEqual(cofm['MIN_MJD'][0],    np.min(fm['MJD'][0:2]))
         self.assertEqual(cofm['MAX_MJD'][0],    np.max(fm['MJD'][0:2]))
         self.assertEqual(cofm['MEAN_MJD'][0],  np.mean(fm['MJD'][0:2]))
-        self.assertEqual(cofm['FIRSTNIGHT'][1], np.min(fm['NIGHT'][2:]))
-        self.assertEqual(cofm['LASTNIGHT'][1],  np.max(fm['NIGHT'][2:]))
         self.assertEqual(cofm['MIN_MJD'][1],    np.min(fm['MJD'][2:]))
         self.assertEqual(cofm['MAX_MJD'][1],    np.max(fm['MJD'][2:]))
         self.assertEqual(cofm['MEAN_MJD'][1],  np.mean(fm['MJD'][2:]))
 
+        #- one target completely masked, but MJD cols still filled
+        fm['TARGETID'] = np.arange(nspec, dtype=int)//2
+        fm['FIBERSTATUS'] = 0
+        fm['FIBERSTATUS'][0:2] = fibermask.BADFIBER
+        cofm, expfm = coadd_fibermap(fm, onetile=True)
+        self.assertTrue(np.all(cofm['MEAN_MJD'] != 0))
+        self.assertTrue(np.all(cofm['MIN_MJD'] != 0))
+        self.assertTrue(np.all(cofm['MAX_MJD'] != 0))
+
 
     def test_coadd_targetmask(self):
         """Test coadding SV1/SV3/DESI_TARGET with varying bits"""