Fujilupe release candidate

.github/workflows/python-package.yml

@@ -87,7 +87,7 @@ jobs:
               env:
                 COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }}
                 GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-              run: coveralls
+              run: coveralls --service=github
 
     docs:
         name: Doc test

bin/CKC2fits

@@ -9,14 +9,18 @@ DESI at ~1 micron) but to save space every just fourth pixel longward of 1
 micron (so an effective pixel spacing of ~100 km/s) is written out.  Each SSP
 has ~83,000 pixels.
 
+version 1.0: 25 km/s pixels
+version 1.1: 10 km/s pixels
+
 """
 import os, pdb
 import numpy as np
 
 from scipy import constants
 C_LIGHT = constants.c / 1000.0 # [km/s]
 
-ssppath = '/global/cfs/cdirs/desi/science/gqp/templates/SSP-CKC14z-loglam'
+ssppath = '/global/cfs/cdirs/desi/science/gqp/templates/SSP-CKC14z'
+version = 'v1.1'
 
 def _CKC2fits(metallicity='Z0.0190'):
     from astropy.io import fits
@@ -30,11 +34,11 @@ def _CKC2fits(metallicity='Z0.0190'):
 
     sspfile = os.path.join(ssppath, 'spec', 'SSP_{}_{}_{}_{}.out.spec'.format(
         isochrone, library, imf, metallicity))
-    outfile = os.path.join(ssppath, 'SSP_{}_{}_{}_{}.fits'.format(
+    outfile = os.path.join(ssppath, version, 'SSP_{}_{}_{}_{}.fits'.format(
         isochrone, library, imf, metallicity))
 
     skiprow = 8 # skip the informational header
-    nage, _npix = np.loadtxt(sspfile, skiprows=skiprow, max_rows=1, dtype=np.int, unpack=True)
+    nage, _npix = np.loadtxt(sspfile, skiprows=skiprow, max_rows=1, dtype=int, unpack=True)
     skiprow += 1
 
     wave = np.loadtxt(sspfile, skiprows=skiprow, max_rows=1) # vacuum wavelengths
@@ -43,7 +47,7 @@ def _CKC2fits(metallicity='Z0.0190'):
 
     # Resample to constant log-lambda / velocity. In the IR (starting at ~1
     # micron), take every fourth sampling, to save space.
-    pixkms = 25.0                            # pixel size [km/s]
+    pixkms = 10.0                            # pixel size [km/s]
     irfactor = 4
     wavesplit = 1e4
     dlogwave = pixkms / C_LIGHT / np.log(10) # pixel size [log-lambda]
@@ -89,6 +93,7 @@ def _CKC2fits(metallicity='Z0.0190'):
 
     hduflux = fits.PrimaryHDU(flux.value)
     hduflux.header['EXTNAME'] = 'FLUX'
+    hduflux.header['VERSION'] = version
     hduflux.header['BUNIT'] = 'erg/(s cm2 Angstrom)'
 
     hduwave = fits.ImageHDU(newwave.value)
@@ -113,4 +118,3 @@ if __name__ == '__main__':
     #for metallicity in [allZ[6]]:
     for metallicity in allZ:
         _CKC2fits(metallicity)
-

bin/fastspecfit-deredshift

@@ -279,8 +279,8 @@ def main(args=None, comm=None):
                     petal, tile, thrunights[indx[0]]))
                 zbestfiles.append(zbestfile)
 
-            Spec.find_specfiles(zbestfiles=zbestfiles, specprod=specprod,
-                                coadd_type='cumulative', targetids=targetids)
+            Spec.select(zbestfiles=zbestfiles, specprod=specprod,
+                        coadd_type='cumulative', targetids=targetids)
             #dataphot = Spec.read_and_unpack(CFit, fastphot=True, synthphot=False)
             dataspec = Spec.read_and_unpack(CFit, fastphot=False, synthphot=False,
                                             remember_coadd=True)

bin/fastspecfit-html

@@ -5,7 +5,7 @@ fastspecfit-html --fastspecfile ./fastspect.fits
 
 fastspecfit-html --fastphotfile /global/cfs/cdirs/desi/spectro/fastspecfit/everest/catalogs/fastphot-everest-sv3-dark.fits \
   --fastspecfile /global/cfs/cdirs/desi/spectro/fastspecfit/everest/catalogs/fastspec-everest-sv3-dark.fits \
-  -o /global/cfs/cdirs/desi/users/ioannis/fastspecfit/everest --hpxpixel 10016
+  -o /global/cfs/cdirs/desi/users/ioannis/fastspecfit/everest --healpix 10016
 
 fastspecfit-html --fastphotfile /global/cfs/cdirs/desi/spectro/fastspecfit/everest/tiles/merged/fastphot-everest-cumulative.fits \
   --fastspecfile /global/cfs/cdirs/desi/spectro/fastspecfit/everest/tiles/merged/fastspec-everest-cumulative.fits \
@@ -30,7 +30,7 @@ def parse(options=None):
 
     parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-    parser.add_argument('--hpxpixel', default=None, type=str, nargs='*', help="""Generate QA for all objects
+    parser.add_argument('--healpix', default=None, type=str, nargs='*', help="""Generate QA for all objects
         with this healpixels (only defined for coadd-type 'healpix').""")
     parser.add_argument('--tile', default=None, type=str, nargs='*', help='Tile(s) to process.')
     parser.add_argument('--night', default=None, type=str, nargs='*', help='Night(s) to process (ignored if coadd-type is cumulative).')
@@ -101,9 +101,9 @@ def main(args=None, comm=None):
             fastphot['NIGHT'] = metaphot['NIGHT']
 
         if coadd_type == 'healpix':
-            joinkeys += ['HPXPIXEL']
-            fastspec['HPXPIXEL'] = metaspec['HPXPIXEL']
-            fastphot['HPXPIXEL'] = metaphot['HPXPIXEL']
+            joinkeys += ['HEALPIX']
+            fastspec['HEALPIX'] = metaspec['HEALPIX']
+            fastphot['HEALPIX'] = metaphot['HEALPIX']
 
         fastfit = join(fastspec, fastphot, keys=joinkeys, join_type='outer',
                        table_names=['SPEC', 'PHOT'])
@@ -138,7 +138,7 @@ def main(args=None, comm=None):
         metadata = metadata[keep]
 
     fastfit, metadata = select(fastfit, metadata, coadd_type,
-                               hpxpixels=args.hpxpixel,
+                               healpixs=args.healpix,
                                tiles=args.tile, nights=args.night)
 
     # read the tile info file for this production
@@ -151,7 +151,7 @@ def main(args=None, comm=None):
     # add tileid --> healpix mapping
     if coadd_type == 'healpix':
         hpxinfo = Table.read(os.path.join(specproddir, 'healpix', 'tilepix.fits'))
-        metadata['TILEID'] = hpxinfo[np.isin(hpxinfo['HEALPIX'], metadata['HPXPIXEL'])]['TILEID']
+        metadata['TILEID'] = hpxinfo[np.isin(hpxinfo['HEALPIX'], metadata['HEALPIX'])]['TILEID']
         #info['HEALPIX'] = hpxinfo[np.isin(info['TILEID'], hpxinfo['TILEID'])]
 
     itile = np.isin(info['TILEID'], metadata['TILEID'])

bin/fastspecfit-qa

@@ -23,24 +23,21 @@ def parse(options=None):
 
     parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-    parser.add_argument('--hpxpixel', default=None, type=str, nargs='*', help="""Generate QA for all objects
+    parser.add_argument('--healpix', default=None, type=str, nargs='*', help="""Generate QA for all objects
         with this healpixels (only defined for coadd-type 'healpix').""")
     parser.add_argument('--tile', default=None, type=str, nargs='*', help='Generate QA for all objects on this tile.')
     parser.add_argument('--night', default=None, type=str, nargs='*', help="""Generate QA for all objects observed on this
         night (only defined for coadd-type 'pernight' and 'perexp').""")
+    parser.add_argument('--redux_dir', default=None, type=str, help='Top-level path to the reduced spectra.')
     parser.add_argument('--targetids', type=str, default=None, help='Comma-separated list of target IDs to process.')
     parser.add_argument('-n', '--ntargets', type=int, help='Number of targets to process in each file.')
     parser.add_argument('--firsttarget', type=int, default=0, help='Index of first object to to process in each file (0-indexed).')
     parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing processes per MPI rank or node.')
-
     parser.add_argument('--nolegend', action='store_true', help='Exclude legends from QA.')
-
     parser.add_argument('--outprefix', default=None, type=str, help='Optional prefix for output filename.')
     parser.add_argument('-o', '--outdir', default='.', type=str, help='Full path to desired output directory.')
 
-    parser.add_argument('--fastphotfile', default=None, type=str, help='Full path to fastphot fitting results.')
-    parser.add_argument('--fastspecfile', default=None, type=str, help='Full path to fastphot fitting results.')
-    #parser.add_argument('redrockdata', nargs='*', help='Full path to redrock file(s) or a directory containing redrock file(s).')
+    parser.add_argument('fastfitfile', nargs=1, help='Full path to fastspec or fastphot fitting results.')
 
     if options is None:
         args = parser.parse_args()
@@ -65,21 +62,18 @@ def main(args=None, comm=None):
     if isinstance(args, (list, tuple, type(None))):
         args = parse(args)
 
+    if args.redux_dir is None:
+        args.redux_dir = os.path.join(os.environ.get('DESI_ROOT', DESI_ROOT_NERSC), 'spectro', 'redux')
+        if not os.path.isdir(args.redux_dir):
+            errmsg = 'Data reduction directory {} not found.'.format(args.redux_dir)
+            log.critical(errmsg)
+            raise IOError(errmsg)
+
     # Read the fitting results.
-    if ((args.fastphotfile is None and args.fastspecfile is None) or
-        (args.fastphotfile is not None and args.fastspecfile is not None)):
-        log.warning('Must provide one (but not both) of --fastphotfile or --fastspecfile.')
+    if not os.path.isfile(args.fastfitfile[0]):
+        log.warning('File {} not found.'.format(args.fastfitfile[0]))
         return
-    elif args.fastphotfile is None and args.fastspecfile is not None:
-        fastphot = False
-        fastfitfile = args.fastspecfile
-    elif args.fastphotfile is not None and args.fastspecfile is None:
-        fastphot = True
-        fastfitfile = args.fastphotfile
-    else:
-        pass # should never get here
-
-    fastfit, metadata, coadd_type = read_fastspecfit(fastfitfile, fastphot=fastphot)
+    fastfit, metadata, coadd_type, fastphot = read_fastspecfit(args.fastfitfile[0])
 
     # parse the targetids optional input
     if args.targetids:
@@ -97,32 +91,23 @@ def main(args=None, comm=None):
         fastfit = fastfit[keep]
         metadata = metadata[keep]
 
-    fastfit, metadata = select(fastfit, metadata, coadd_type, hpxpixels=args.hpxpixel,
+    fastfit, metadata = select(fastfit, metadata, coadd_type, healpixels=args.healpix,
                                tiles=args.tile, nights=args.night)
 
     # Initialize the continuum- and emission-line fitting classes.
     t0 = time.time()
     CFit = ContinuumFit(nolegend=args.nolegend)
     EMFit = EMLineFit(nolegend=args.nolegend)
-    Spec = DESISpectra()
+    Spec = DESISpectra(redux_dir=args.redux_dir)
     log.info('Initializing the classes took: {:.2f} sec'.format(time.time()-t0))
 
     if args.outdir:
         if not os.path.isdir(args.outdir):
             os.makedirs(args.outdir, exist_ok=True)
 
-    print('Hacking specprod!!!')
-    specprod = 'fuji'
-    specproddir = os.path.join(os.environ.get('DESI_ROOT', DESI_ROOT_NERSC), 'spectro', 'redux', specprod)
-
-    #if fastphot:
-    #    synthphot = False
-    #else:
-    #    synthphot = True
-
     def _wrap_qa(redrockfile, indx):
         targetids = fastfit['TARGETID'][indx]
-        Spec.find_specfiles(redrockfiles=redrockfile, targetids=targetids)
+        Spec.select(redrockfiles=redrockfile, targetids=targetids)
         data = Spec.read_and_unpack(CFit, fastphot=fastphot, synthphot=True, remember_coadd=True)
         qaargs = [(CFit, EMFit, data[igal], fastfit[indx[igal]], metadata[indx[igal]],
                    coadd_type, fastphot, args.outdir, args.outprefix) for igal in np.arange(len(indx))]                
@@ -136,58 +121,68 @@ def main(args=None, comm=None):
 
     t0 = time.time()
     if coadd_type == 'healpix':
+        allspecprods = metadata['SPECPROD'].data
         allsurveys = metadata['SURVEY'].data
-        allprograms = metadata['FAPRGRM'].data
-        allpixels = metadata['HPXPIXEL'].data
-        for survey in set(allsurveys):
-            for program in set(allprograms):
-                for pixel in set(allpixels):
-                    indx = np.where((survey == allsurveys) * (program == allprograms) * (pixel == allpixels))[0]
-                    if len(indx) == 0:
-                        log.warning('No object found with survey={}, program={}, and healpixel={}!'.format(
-                            survey, program, pixel))
-                        continue
-                    redrockfile = os.path.join(specproddir, 'healpix', str(survey), str(program), str(pixel // 100), str(pixel), 
-                                               'redrock-{}-{}-{}.fits'.format(survey, program, pixel))
-                    _wrap_qa(redrockfile, indx)
+        allprograms = metadata['PROGRAM'].data
+        allpixels = metadata['HEALPIX'].data
+        for specprod in set(allspecprods):
+            for survey in set(allsurveys):
+                for program in set(allprograms):
+                    for pixel in set(allpixels):
+                        indx = np.where((specprod == allspecprods) * (survey == allsurveys) *
+                                        (program == allprograms) * (pixel == allpixels))[0]
+                        if len(indx) == 0:
+                            #log.warning('No object found with specprod={}, survey={}, program={}, and healpixel={}!'.format(
+                            #    specprod, survey, program, pixel))
+                            continue
+                        redrockfile = os.path.join(args.redux_dir, specprod, 'healpix', str(survey), str(program), str(pixel // 100), 
+                                                   str(pixel), 'redrock-{}-{}-{}.fits'.format(survey, program, pixel))
+                        _wrap_qa(redrockfile, indx)
     else:
+        allspecprods = metadata['SPECPROD'].data
         alltiles = metadata['TILEID'].astype(str).data
+        allnights = metadata['NIGHT'].astype(str).data
         allpetals = metadata['FIBER'].data // 500
         if coadd_type == 'cumulative':
-            thrunights = metadata['THRUNIGHT'].astype(str).data
-            for tile in set(alltiles):
-                for petal in set(allpetals):
-                    indx = np.where((tile == alltiles) * (petal == allpetals))[0]
-                    if len(indx) == 0:
-                        log.warning('No object found with tileid={} and petal={}!'.format(
-                            tile, petal))
-                        continue
-                    redrockfile = os.path.join(specproddir, 'tiles', 'cumulative', str(tile), thrunights[indx[0]],
-                                               'redrock-{}-{}-thru{}.fits'.format(petal, tile, thrunights[indx[0]]))
-                    _wrap_qa(redrockfile, indx)
-        elif coadd_type == 'pernight':
-            allnights = metadata['NIGHT'].astype(str).data
-            for night in set(allnights):
+            for specprod in set(allspecprods):
                 for tile in set(alltiles):
                     for petal in set(allpetals):
-                        indx = np.where((night == allnights) * (tile == alltiles) * (petal == allpetals))[0]
-                        redrockfile = os.path.join(specproddir, 'tiles', 'pernight', str(tile), str(night),
-                                                   'redrock-{}-{}-{}.fits'.format(petal, tile, night))
+                        indx = np.where((specprod == allspecprods) * (tile == alltiles) * (petal == allpetals))[0]
+                        if len(indx) == 0:
+                            #log.warning('No object found with tileid={} and petal={}!'.format(
+                            #    tile, petal))
+                            continue
+                        redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'cumulative', str(tile), allnights[indx[0]],
+                                                   'redrock-{}-{}-thru{}.fits'.format(petal, tile, allnights[indx[0]]))
                         _wrap_qa(redrockfile, indx)
-        elif coadd_type == 'perexp':
-            allnights = metadata['NIGHT'].astype(str).data
-            allexpids = metadata['EXPID'].data
-            for night in set(allnights):
-                for expid in set(allexpids):
+        elif coadd_type == 'pernight':
+            for specprod in set(allspecprods):
+                for night in set(allnights):
                     for tile in set(alltiles):
                         for petal in set(allpetals):
-                            indx = np.where((night == allnights) * (expid == allexpids) *
+                            indx = np.where((specprod == allspecprods) * (night == allnights) *
                                             (tile == alltiles) * (petal == allpetals))[0]
-                            redrockfile = os.path.join(specproddir, 'tiles', 'perexp', str(tile), '{:08d}'.format(expid), 
-                                                       'redrock-{}-{}-exp{:08d}.fits'.format(petal, tile, expid))
+                            if len(indx) == 0:
+                                continue
+                            redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'pernight', str(tile), str(night),
+                                                       'redrock-{}-{}-{}.fits'.format(petal, tile, night))
                             _wrap_qa(redrockfile, indx)
-        else:
-            pass
+        elif coadd_type == 'perexp':
+            allexpids = metadata['EXPID'].data
+            for specprod in set(allspecprods):
+                for night in set(allnights):
+                    for expid in set(allexpids):
+                        for tile in set(alltiles):
+                            for petal in set(allpetals):
+                                indx = np.where((specprod == allspecprods) * (night == allnights) *
+                                                (expid == allexpids) * (tile == alltiles) *
+                                                (petal == allpetals))[0]
+                                if len(indx) == 0:
+                                    continue
+                                redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'perexp', str(tile), '{:08d}'.format(expid), 
+                                                           'redrock-{}-{}-exp{:08d}.fits'.format(petal, tile, expid))
+                                _wrap_qa(redrockfile, indx)
+
     log.info('QA for everything took: {:.2f} sec'.format(time.time()-t0))
 
 if __name__ == '__main__':