Merge pull request #1647 from pypeit/telpca

PCA method for telluric corrections

CHANGES.rst

@@ -58,8 +58,7 @@ RELEASE FILE IN doc/releases
       files that have frametype None (this prevent ``run_pypeit`` to crash)
     - Added a function ``check_spectrograph()`` (currently only defined for LRIS),
       that checks (during ``pypeit_setup``) if the selected spectrograph is the
-      corrected one for the data used.
-
+      corrected one for the data used. 
 
 1.13.0 (2 June 2023)
 --------------------

doc/releases/1.14.1dev.rst

@@ -43,6 +43,7 @@ Functionality/Performance Improvements and Additions
   to reflect the *type* of polynomial being fit.  Both names point to the same
   code, but the name ``'polynomial'`` is deprecated and will be removed in the
   future.
+- Introduced PCA method for telluric corrections
 - Added a new GUI for creating and editing PypeIt input files: ``pypeit_setup_gui``
 
 Instrument-specific Updates

pypeit/par/pypeitpar.py

@@ -2209,7 +2209,7 @@ class TelluricPar(ParSet):
     """
 
     def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_bounds=None, pix_shift_bounds=None,
-                 delta_coeff_bounds=None, minmax_coeff_bounds=None, maxiter=None,
+                 delta_coeff_bounds=None, minmax_coeff_bounds=None, maxiter=None, tell_npca=None, teltype=None,
                  sticky=None, lower=None, upper=None, seed=None, tol=None, popsize=None, recombination=None, polish=None,
                  disp=None, objmodel=None, redshift=None, delta_redshift=None, pca_file=None, npca=None,
                  bal_wv_min_max=None, bounds_norm=None, tell_norm_thresh=None, only_orders=None, pca_lower=None,
@@ -2224,6 +2224,7 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
         # Initialize the other used specifications for this parameter
         # set
         defaults = OrderedDict.fromkeys(pars.keys())
+        options = OrderedDict.fromkeys(pars.keys())
         dtypes = OrderedDict.fromkeys(pars.keys())
         descr = OrderedDict.fromkeys(pars.keys())
 
@@ -2234,6 +2235,18 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
                                'must be downloaded from the GoogleDrive and installed in your PypeIt installation via ' \
                                'the pypeit_install_telluric script. NOTE: This parameter no longer includes the full ' \
                                'pathname to the Telluric Grid file, but is just the filename of the grid itself.'
+
+        defaults['tell_npca'] = 5
+        dtypes['tell_npca'] = int
+        descr['tell_npca'] = 'Number of telluric PCA components used. Can be set to any number from 1 to 10.'
+
+        defaults['teltype'] = 'pca'
+        options['teltype'] = TelluricPar.valid_teltype()
+        dtypes['teltype'] = str
+        descr['teltype'] = 'Method used to evaluate telluric models, either pca or grid. The grid option uses a ' \
+                           'fixed grid of pre-computed HITRAN+LBLRTM atmospheric transmission models for each ' \
+                           'observatory, whereas the pca option uses principal components of a larger model grid ' \
+                           'to compute an accurate pseudo-telluric model with a much lighter telgridfile.'
 
         defaults['sn_clip'] = 30.0
         dtypes['sn_clip'] = [int, float]
@@ -2254,12 +2267,12 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
 
 
         pars['resln_frac_bounds'] = tuple_force(pars['resln_frac_bounds'])
-        defaults['resln_frac_bounds'] = (0.5,1.5)
+        defaults['resln_frac_bounds'] = (0.6,1.4)
         dtypes['resln_frac_bounds'] = tuple
         descr['resln_frac_bounds'] = 'Bounds for the resolution fit optimization which is part of the telluric model. ' \
-                                     'This range is in units of the resln_guess, so the (0.5, 1.5) would bound the ' \
+                                     'This range is in units of the resln_guess, so the (0.6, 1.4) would bound the ' \
                                      'spectral resolution fit to be within the range ' \
-                                     'bounds_resln = (0.5*resln_guess, 1.5*resln_guess)'
+                                     'bounds_resln = (0.6*resln_guess, 1.4*resln_guess)'
 
         pars['pix_shift_bounds'] = tuple_force(pars['pix_shift_bounds'])
         defaults['pix_shift_bounds'] = (-5.0,5.0)
@@ -2462,7 +2475,7 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
     @classmethod
     def from_dict(cls, cfg):
         k = np.array([*cfg.keys()])
-        parkeys = ['telgridfile', 'sn_clip', 'resln_guess', 'resln_frac_bounds',
+        parkeys = ['telgridfile', 'teltype', 'sn_clip', 'resln_guess', 'resln_frac_bounds', 'tell_npca',
                    'pix_shift_bounds', 'delta_coeff_bounds', 'minmax_coeff_bounds',
                    'maxiter', 'sticky', 'lower', 'upper', 'seed', 'tol',
                    'popsize', 'recombination', 'polish', 'disp', 'objmodel','redshift', 'delta_redshift',
@@ -2479,12 +2492,27 @@ def from_dict(cls, cfg):
         for pk in parkeys:
             kwargs[pk] = cfg[pk] if pk in k else None
         return cls(**kwargs)
+
+    @staticmethod
+    def valid_teltype():
+        """
+        Return the valid telluric methods.
+        """
+        return ['pca', 'grid']
 
     def validate(self):
         """
         Check the parameters are valid for the provided method.
         """
-        pass
+        if self.data['tell_npca'] < 1 or self.data['tell_npca'] > 10:
+            raise ValueError('Invalid value {:d} for tell_npca '.format(self.data['tell_npca'])+
+                             '(must be between 1 and 10).')
+
+        self.data['teltype'] = self.data['teltype'].lower()
+        if self.data['teltype'] not in TelluricPar.valid_teltype():
+            raise ValueError('Invalid teltype "{}"'.format(self.data['teltype'])+
+                             ', valid options are: {}.'.format(TelluricPar.valid_teltype()))
+
         # JFH add something in here which checks that the recombination value provided is bewteen 0 and 1, although
         # scipy.optimize.differential_evoluiton probalby checks this.
 

pypeit/scripts/tellfit.py

@@ -155,13 +155,17 @@ def main(args):
                                            par['telluric']['pca_file'],
                                            par['telluric']['redshift'], modelfile, outfile,
                                            npca=par['telluric']['npca'],
+                                           teltype=par['telluric']['teltype'], tell_npca=par['telluric']['tell_npca'],
                                            pca_lower=par['telluric']['pca_lower'],
                                            pca_upper=par['telluric']['pca_upper'],
                                            bounds_norm=par['telluric']['bounds_norm'],
                                            tell_norm_thresh=par['telluric']['tell_norm_thresh'],
                                            only_orders=par['telluric']['only_orders'],
                                            bal_wv_min_max=par['telluric']['bal_wv_min_max'],
+                                           pix_shift_bounds=par['telluric']['pix_shift_bounds'],
                                            maxiter=par['telluric']['maxiter'],
+                                           popsize=par['telluric']['popsize'],
+                                           tol=par['telluric']['tol'],
                                            debug_init=args.debug, disp=args.debug,
                                            debug=args.debug, show=args.plot)
         elif par['telluric']['objmodel']=='star':
@@ -175,12 +179,16 @@ def main(args):
                                              model=par['telluric']['model'],
                                              polyorder=par['telluric']['polyorder'],
                                              only_orders=par['telluric']['only_orders'],
+                                             teltype=par['telluric']['teltype'], tell_npca=par['telluric']['tell_npca'],
                                              mask_hydrogen_lines=par['sensfunc']['mask_hydrogen_lines'],
                                              mask_helium_lines=par['sensfunc']['mask_helium_lines'],
                                              hydrogen_mask_wid=par['sensfunc']['hydrogen_mask_wid'],
                                              delta_coeff_bounds=par['telluric']['delta_coeff_bounds'],
                                              minmax_coeff_bounds=par['telluric']['minmax_coeff_bounds'],
+                                             pix_shift_bounds=par['telluric']['pix_shift_bounds'],
                                              maxiter=par['telluric']['maxiter'],
+                                             popsize=par['telluric']['popsize'],
+                                             tol=par['telluric']['tol'],
                                              debug_init=args.debug, disp=args.debug,
                                              debug=args.debug, show=args.plot)
         elif par['telluric']['objmodel']=='poly':
@@ -190,12 +198,16 @@ def main(args):
                                              func=par['telluric']['func'],
                                              model=par['telluric']['model'],
                                              polyorder=par['telluric']['polyorder'],
+                                             teltype=par['telluric']['teltype'], tell_npca=par['telluric']['tell_npca'],
                                              fit_wv_min_max=par['telluric']['fit_wv_min_max'],
                                              mask_lyman_a=par['telluric']['mask_lyman_a'],
                                              delta_coeff_bounds=par['telluric']['delta_coeff_bounds'],
                                              minmax_coeff_bounds=par['telluric']['minmax_coeff_bounds'],
                                              only_orders=par['telluric']['only_orders'],
+                                             pix_shift_bounds=par['telluric']['pix_shift_bounds'],
                                              maxiter=par['telluric']['maxiter'],
+                                             popsize=par['telluric']['popsize'],
+                                             tol=par['telluric']['tol'],
                                              debug_init=args.debug, disp=args.debug,
                                              debug=args.debug, show=args.plot)
         else:

pypeit/sensfunc.py

@@ -885,7 +885,10 @@ def compute_zeropoint(self):
                                                    ech_orders=self.meta_spec['ECH_ORDERS'],
                                                    resln_guess=self.par['IR']['resln_guess'],
                                                    resln_frac_bounds=self.par['IR']['resln_frac_bounds'],
+                                                   pix_shift_bounds=self.par['IR']['pix_shift_bounds'],
                                                    sn_clip=self.par['IR']['sn_clip'],
+                                                   teltype=self.par['IR']['teltype'],
+                                                   tell_npca=self.par['IR']['tell_npca'],
                                                    mask_hydrogen_lines=self.par['mask_hydrogen_lines'],
                                                    maxiter=self.par['IR']['maxiter'],
                                                    lower=self.par['IR']['lower'],

pypeit/spectrographs/gemini_flamingos.py

@@ -143,7 +143,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
         # TODO: replace the telluric grid file for Gemini-S site.
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 

pypeit/spectrographs/gemini_gmos.py

@@ -873,7 +873,7 @@ def default_pypeit_par(cls):
         """
         par = super().default_pypeit_par()
         par['sensfunc']['algorithm'] = 'IR'
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         # Bound the detector with slit edges if no edges are found. These data are often trimmed
         # so we implement this here as the default.
         par['calibrations']['slitedges']['bound_detector'] = True

pypeit/spectrographs/gemini_gnirs.py

@@ -297,7 +297,7 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 6
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     def config_specific_par(self, scifile, inp_par=None):

pypeit/spectrographs/gtc_osiris.py

@@ -384,7 +384,7 @@ def config_specific_par(self, scifile, inp_par=None):
             par['calibrations']['wavelengths']['lamps'] = ['ArI','XeI','NeI']
             par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500I.fits'
             par['sensfunc']['algorithm'] = 'IR'
-            par['sensfunc']['IR']['telgridfile'] = "TelFit_MaunaKea_3100_26100_R20000.fits"
+            par['sensfunc']['IR']['telgridfile'] = "TellPCA_3000_26000_R10000.fits"
         else:
             msgs.warn('gtc_osiris.py: template arc missing for this grism! Trying holy-grail...')
             par['calibrations']['wavelengths']['method'] = 'holy-grail'
@@ -955,7 +955,7 @@ def config_specific_par(self, scifile, inp_par=None):
             par['calibrations']['wavelengths']['lamps'] = ['ArI','XeI','NeI']
             par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500I.fits'
             par['sensfunc']['algorithm'] = 'IR'
-            par['sensfunc']['IR']['telgridfile'] = "TelFit_MaunaKea_3100_26100_R20000.fits"
+            par['sensfunc']['IR']['telgridfile'] = "TellPCA_3000_26000_R10000.fits"
         else:
             msgs.warn('gtc_osiris.py: template arc missing for this grism! Trying holy-grail...')
             par['calibrations']['wavelengths']['method'] = 'holy-grail'

pypeit/spectrographs/keck_deimos.py

@@ -301,7 +301,7 @@ def default_pypeit_par(cls):
         par['scienceframe']['process']['objlim'] = 1.5
 
         # If telluric is triggered
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R15000.fits'
         return par
 
     def config_specific_par(self, scifile, inp_par=None):

pypeit/spectrographs/keck_hires.py

@@ -157,7 +157,14 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 5 #[9, 11, 11, 9, 9, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7]
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_10500_R120000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-40.0,40.0)
+
+        # Telluric parameters
+        # HIRES is usually oversampled, so the helio shift can be large
+        par['telluric']['pix_shift_bounds'] = (-40.0,40.0)
+        # Similarly, the resolution guess is higher than it should be
+        par['telluric']['resln_frac_bounds'] = (0.25,1.25)
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'

pypeit/spectrographs/keck_lris.py

@@ -101,7 +101,7 @@ def default_pypeit_par(cls):
 
 
         # If telluric is triggered
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 

pypeit/spectrographs/keck_mosfire.py

@@ -127,7 +127,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 13
         par['sensfunc']['IR']['maxiter'] = 2
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     # NOTE: This function is used by the dev-suite

pypeit/spectrographs/keck_nires.py

@@ -140,7 +140,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
         par['sensfunc']['IR']['maxiter'] = 2
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'

pypeit/spectrographs/keck_nirspec.py

@@ -129,7 +129,12 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R25000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-8.0,8.0)
+
+        # Telluric parameters
+        par['telluric']['pix_shift_bounds'] = (-8.0,8.0)
+
         return par
 
     def init_meta(self):

pypeit/spectrographs/magellan_fire.py

@@ -201,7 +201,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['polyorder'] = 5
         par['sensfunc']['IR']['maxiter'] = 2
         # place holder for telgrid file
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R15000.fits'
 
         # Coadding. I'm not sure what this should be for PRISM mode?
         par['coadd1d']['wave_method'] = 'log10'
@@ -418,7 +418,7 @@ def default_pypeit_par(cls):
         par['reduce']['findobj']['find_trim_edge'] = [50,50]
         par['flexure']['spec_method'] = 'skip'
 
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         # Set the default exposure time ranges for the frame typing
         par['calibrations']['standardframe']['exprng'] = [None, 60]

pypeit/spectrographs/mmt_binospec.py

@@ -214,7 +214,7 @@ def default_pypeit_par(cls):
 
         # Sensitivity function parameters
         par['sensfunc']['polyorder'] = 7
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 

pypeit/spectrographs/mmt_mmirs.py

@@ -202,7 +202,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
         # ToDo: replace the telluric grid file for MMT site.
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 

pypeit/spectrographs/p200_dbsp.py

@@ -525,7 +525,7 @@ def default_pypeit_par(cls):
 
         par['sensfunc']['algorithm'] = 'UVIS'
         par['sensfunc']['UVIS']['polycorrect'] = False
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Lick_3100_11100_R10000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     def config_specific_par(self, scifile, inp_par=None):

pypeit/spectrographs/p200_tspec.py

@@ -215,7 +215,7 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'