Provide wrapper to instantiate model spectrum

[biprateep]

It will be helpful to provide a wrapper to take a data table with the continuum coefficients as inputs and instantiate the best fit model spectrum for both photometry only (FASTPHOT) and spectra only (FASTSPEC) fits.

This can help with science cases like (but not limited to):

• Synthesize photometry in any arbitrary band
• apply aperture corrections
• change distance modulus to any arbitrary cosmology

[moustakas]

Related to this ticket, @Ragadeepika-Pucha has asked how to generate the standard QA for fastspecfit. Although the code base is run within a dedicated Docker container, the standard DESI software stack satisfies all the dependencies, so this is trivial.

First, clone a local copy of the repository (be sure to clone by SSH if you want to develop/contribute):

cd /path/to/code
git clone https://github.com/desihub/fastspecfit.git

then do the following (which can be put in a simple shell script):

source /global/cfs/cdirs/desi/software/desi_environment.sh master
package=fastspecfit
export PATH=/path/to/code/$package/bin:${PATH}
export PYTHONPATH=/path/to/code/$package/py:${PYTHONPATH}

and then run it! E.g., to reproduce the plots from #50 do:

fastspecfit-qa --fastspecfile /global/cfs/cdirs/desi/spectro/fastspecfit/everest/catalogs/fastspec-everest-main-bright.fits \
  --targetids 39633006721237481,39632971702995257 -o ./

Note that you must specify the full-path location to either a fastspecfile (here, fastspec-everest-main-bright.fits) or a fastphotfile (e.g., fastphot-everest-main-bright.fits) and then code will make the appropriate QA. It will also figure out (and go read) the original (reduced) DESI spectra from the metadata stored in the fastspecfit files.

Also: if the TARGETIDS of the object(s) you're interested in aren't in the file you specify, then the code will exit gracefully. But remember that it will make QA for every row in the input file, so be careful about handing it a file with millions of galaxies! (And you can also invoke parallelism with --mp 32, e.g., on an interactive node.)

For the record, here are all the options of the QA code (which needs to be documented better, as always):

ioannis@cori02:~% fastspecfit-qa -h
usage: fastspecfit-qa [-h] [--hpxpixel [HPXPIXEL [HPXPIXEL ...]]] [--tile [TILE [TILE ...]]] [--night [NIGHT [NIGHT ...]]] [--targetids TARGETIDS] [-n NTARGETS]
                      [--firsttarget FIRSTTARGET] [--mp MP] [--outprefix OUTPREFIX] [-o OUTDIR] [--fastphotfile FASTPHOTFILE] [--fastspecfile FASTSPECFILE]

optional arguments:
  -h, --help            show this help message and exit
  --hpxpixel [HPXPIXEL [HPXPIXEL ...]]
                        Generate QA for all objects with this healpixels (only defined for coadd-type 'healpix'). (default: None)
  --tile [TILE [TILE ...]]
                        Generate QA for all objects on this tile. (default: None)
  --night [NIGHT [NIGHT ...]]
                        Generate QA for all objects observed on this night (only defined for coadd-type 'pernight' and 'perexp'). (default: None)
  --targetids TARGETIDS
                        Comma-separated list of target IDs to process. (default: None)
  -n NTARGETS, --ntargets NTARGETS
                        Number of targets to process in each file. (default: None)
  --firsttarget FIRSTTARGET
                        Index of first object to to process in each file (0-indexed). (default: 0)
  --mp MP               Number of multiprocessing processes per MPI rank or node. (default: 1)
  --outprefix OUTPREFIX
                        Optional prefix for output filename. (default: None)
  -o OUTDIR, --outdir OUTDIR
                        Full path to desired output directory. (default: .)
  --fastphotfile FASTPHOTFILE
                        Full path to fastphot fitting results. (default: None)
  --fastspecfile FASTSPECFILE
                        Full path to fastphot fitting results. (default: None)

[stephjuneau]

This is great, @moustakas ! My only question so far is how to specify unique spectra when there are several healpix-spectra for a given TARGETID? The combination of identifiers is TARGETID,SURVEY,FAPRGRM to find unique spectra but having all 3 doesn't seem to be an expected input.

EDIT: I consider those three columns together as the unique identifier. So when I use the fastspecfit catalogs, I pull SURVEY ans FAPRGRM from HDU2 and copy them also in HDU1 because then I can match the table with any other healpix-based tables by joining on all three.

[Ragadeepika-Pucha]

@stephjuneau -- the above script takes the fastspecfit filename as input -- which is unique for a given SURVEY and FAPRGRM. This means, for each file the you have unique objects.

@moustakas -- can you point me to where the code for this is written? Is it okay if I use that code and manipulate it a bit to suit my requirement? Thank you.

[Ragadeepika-Pucha]

@moustakas - i tried what you mentioned above, but it is giving the following error -
If 'fastspecfit-qa' is not a typo you can use command-not-found to lookup the package that contains it, like this: cnf fastspecfit-qa

[moustakas]

@Ragadeepika-Pucha @stephjuneau @biprateep

I've decided to implement this feature request by writing out the best-fitting models with the fastspec tables themselves. I think this closes this ticket, although I'll be sure to add these instructions to the documentation (still forthcoming). I'm also eventually planning to provide instructions as to how anyone can regenerate the models from the output tables and some simple code, too, although that process is a bit more involved.

Here's one simple example based on a test run. After I wrote out these instructions I decided to add EBV to the output metadata table (in 60a97ae) so the final instructions will be even simpler. Note that the models are written out on a coadded-camera wavelength grid.

import os
import numpy as np
import fitsio 
from astropy.table import Table
import matplotlib.pyplot as plt

from desiutil.dust import SFDMap
from desiutil.dust import dust_transmission
from desispec.io import read_spectra
from desispec.coaddition import coadd_cameras

fastfile = '/global/cfs/cdirs/desi/users/ioannis/fastspecfit/vitiles/fuji/tiles/cumulative/80613/20210324/fastspec-3-80613-thru20210324.fits.gz'
specfile = os.path.join(os.environ.get('DESI_ROOT'), 'spectro', 'redux', 'fuji', 'tiles', 'cumulative', '80613', '20210324', 'coadd-3-80613-thru20210324.fits')

meta = Table(fitsio.read(fastfile, 'METADATA'))
fast = Table(fitsio.read(fastfile, 'FASTSPEC'))

models, hdr = fitsio.read(fastfile, 'MODELS', header=True)
modelwave = hdr['CRVAL1'] + np.arange(hdr['NAXIS1']) * hdr['CDELT1']

spec = read_spectra(specfile).select(targets=meta['TARGETID'])
coadd_spec = coadd_cameras(spec)
bands = coadd_spec.bands[0]

indx = 216

ebv = SFDMap(scaling=1.0).ebv(meta['RA'][indx], meta['DEC'][indx])
mw_transmission_spec = dust_transmission(coadd_spec.wave[bands], ebv)

plt.clf()
plt.plot(coadd_spec.wave[bands], coadd_spec.flux[bands][indx, :] / mw_transmission_spec, label='Data')
plt.plot(modelwave, np.sum(models[indx, :, :], axis=0), label='Final Model')
plt.plot(modelwave, models[indx, 0, :], label='Stellar Continuum Model')
plt.plot(modelwave, models[indx, 1, :], label='Smooth Continuum Correction')
plt.plot(modelwave, models[indx, 2, :], label='Emission-Line Model')
plt.legend(fontsize=8)
plt.xlabel(r'Observed-frame Wavelength ($\AA$)')
plt.ylabel(r'Flux Density ($10^{-17}~{\rm erg}~{\rm s}^{-1}~{\rm cm}^{-2}~\AA^{-1}$)')
plt.show()