From 3570d5bcb473a2c001f4b51931d01c5b0d821d37 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Thu, 25 Jul 2024 13:28:55 -1000 Subject: [PATCH 01/17] added use_std_trace parameter --- pypeit/par/pypeitpar.py | 33 ++++++++++++++++++++++----------- pypeit/pypeit.py | 7 +++++-- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index d9f42ec6dd..f47ff297a7 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -4052,7 +4052,8 @@ def __init__(self, trace_npoly=None, snr_thresh=None, find_trim_edge=None, find_maxdev=None, find_extrap_npoly=None, maxnumber_sci=None, maxnumber_std=None, find_fwhm=None, ech_find_max_snr=None, ech_find_min_snr=None, ech_find_nabove_min_snr=None, skip_second_find=None, skip_final_global=None, - skip_skysub=None, find_negative=None, find_min_max=None, std_spec1d=None, fof_link = None): + skip_skysub=None, find_negative=None, find_min_max=None, std_spec1d=None, + use_std_trace=None, fof_link = None): # Grab the parameter names and values from the function # arguments args, _, _, values = inspect.getargvalues(inspect.currentframe()) @@ -4175,14 +4176,22 @@ def __init__(self, trace_npoly=None, snr_thresh=None, find_trim_edge=None, 'detector. It only used for object finding. This parameter is helpful if your object only ' \ 'has emission lines or at high redshift and the trace only shows in part of the detector.' + # TODO do we want this to be True or False by default? + defaults['use_std_trace'] = False + dtypes['use_std_trace'] = bool + descr['use_std_trace'] = 'If True, the trace of the standard star spectrum is used as a crutch for ' \ + 'tracing the object spectra. This is useful when a direct trace is not possible ' \ + '(i.e., faint sources). Note that a standard star exposure must be included in your ' \ + 'pypeit file, or the ``std_spec1d`` parameter must be set for this to work. ' + + defaults['std_spec1d'] = None dtypes['std_spec1d'] = str - descr['std_spec1d'] = 'A PypeIt spec1d file of a previously reduced standard star. The ' \ - 'trace of the standard star spectrum is used as a crutch for ' \ - 'tracing the object spectra, when a direct trace is not possible ' \ - '(i.e., faint sources). If provided, this overrides use of any ' \ - 'standards included in your pypeit file; the standard exposures ' \ - 'will still be reduced.' + descr['std_spec1d'] = 'A PypeIt spec1d file of a previously reduced standard star. ' \ + 'This can be used to trace the object spectra, but the ``use_std_trace`` ' \ + 'parameter must be set to True. If provided, this overrides use of ' \ + 'any standards included in your pypeit file; the standard exposures ' \ + 'will still be reduced.' # Instantiate the parameter set super(FindObjPar, self).__init__(list(pars.keys()), @@ -4203,7 +4212,7 @@ def from_dict(cls, cfg): 'find_maxdev', 'find_fwhm', 'ech_find_max_snr', 'ech_find_min_snr', 'ech_find_nabove_min_snr', 'skip_second_find', 'skip_final_global', 'skip_skysub', 'find_negative', 'find_min_max', - 'std_spec1d', 'fof_link'] + 'std_spec1d', 'use_std_trace', 'fof_link'] badkeys = np.array([pk not in parkeys for pk in k]) if np.any(badkeys): @@ -4215,9 +4224,11 @@ def from_dict(cls, cfg): return cls(**kwargs) def validate(self): - if self.data['std_spec1d'] is not None \ - and not Path(self.data['std_spec1d']).absolute().exists(): - msgs.error(f'{self.data["std_spec1d"]} does not exist!') + if self.data['std_spec1d'] is not None: + if not self.data['use_std_trace']: + msgs.error('If you provide a standard star spectrum for tracing, you must set use_std_trace=True.') + elif not Path(self.data['std_spec1d']).absolute().exists(): + msgs.error(f'{self.data["std_spec1d"]} does not exist!') class SkySubPar(ParSet): diff --git a/pypeit/pypeit.py b/pypeit/pypeit.py index 6505a0ef9d..a32e85463d 100644 --- a/pypeit/pypeit.py +++ b/pypeit/pypeit.py @@ -245,7 +245,9 @@ def get_std_outfile(self, standard_frames): # isolate where the name of the standard-star spec1d file is defined. std_outfile = self.par['reduce']['findobj']['std_spec1d'] if std_outfile is not None: - if not Path(std_outfile).absolute().exists(): + if not self.par['reduce']['findobj']['use_std_trace']: + msgs.error('If you provide a standard star spectrum for tracing, you must set use_std_trace=True') + elif not Path(std_outfile).absolute().exists(): msgs.error(f'Provided standard spec1d file does not exist: {std_outfile}') return std_outfile @@ -254,7 +256,8 @@ def get_std_outfile(self, standard_frames): # standard associated with a given science frame. Below, I # just use the first standard - std_frame = None if len(standard_frames) == 0 else standard_frames[0] + std_frame = None if (len(standard_frames) == 0 or not self.par['reduce']['findobj']['use_std_trace']) \ + else standard_frames[0] # Prepare to load up standard? if std_frame is not None: std_outfile = self.spec_output_file(std_frame) \ From f79da49545b022314313eb4d439c3b8daf57188f Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Thu, 25 Jul 2024 15:46:12 -1000 Subject: [PATCH 02/17] fix None flux bug + added box extr for sensfunc --- pypeit/par/pypeitpar.py | 19 ++++++++++++++++--- pypeit/scripts/sensfunc.py | 16 ++++++++++++++++ pypeit/sensfunc.py | 6 ++++-- pypeit/specobjs.py | 8 +++++--- 4 files changed, 41 insertions(+), 8 deletions(-) diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index f47ff297a7..a57f34bf3e 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -1921,7 +1921,7 @@ class SensFuncPar(ParSet): see :ref:`parameters`. """ def __init__(self, flatfile=None, extrap_blu=None, extrap_red=None, samp_fact=None, multi_spec_det=None, algorithm=None, UVIS=None, - IR=None, polyorder=None, star_type=None, star_mag=None, star_ra=None, + IR=None, polyorder=None, star_type=None, star_mag=None, star_ra=None, extr=None, star_dec=None, mask_hydrogen_lines=None, mask_helium_lines=None, hydrogen_mask_wid=None): # Grab the parameter names and values from the function arguments args, _, _, values = inspect.getargvalues(inspect.currentframe()) @@ -1939,6 +1939,11 @@ def __init__(self, flatfile=None, extrap_blu=None, extrap_red=None, samp_fact=No 'function computed from a flat field file in the Calibrations directory, e.g.' \ 'Calibrations/Flat_A_0_DET01.fits' + defaults['extr'] = 'OPT' + dtypes['extr'] = str + descr['extr'] = 'Extraction method to use for the sensitivity function. Options are: ' \ + '\'OPT\' (optimal extraction), \'BOX\' (boxcar extraction). Default is \'OPT\'.' + defaults['extrap_blu'] = 0.1 dtypes['extrap_blu'] = float descr['extrap_blu'] = 'Fraction of minimum wavelength coverage to grow the wavelength coverage of the ' \ @@ -2028,7 +2033,7 @@ def __init__(self, flatfile=None, extrap_blu=None, extrap_red=None, samp_fact=No options=list(options.values()), dtypes=list(dtypes.values()), descr=list(descr.values())) -# self.validate() + self.validate() @classmethod def from_dict(cls, cfg): @@ -2036,7 +2041,7 @@ def from_dict(cls, cfg): # Single element parameters parkeys = ['flatfile', 'extrap_blu', 'extrap_red', 'samp_fact', 'multi_spec_det', 'algorithm', - 'polyorder', 'star_type', 'star_mag', 'star_ra', 'star_dec', + 'polyorder', 'star_type', 'star_mag', 'star_ra', 'star_dec', 'extr', 'mask_hydrogen_lines', 'mask_helium_lines', 'hydrogen_mask_wid'] # All parameters, including nested ParSets @@ -2058,6 +2063,14 @@ def from_dict(cls, cfg): return cls(**kwargs) + def validate(self): + """ + Check the parameters are valid for the provided method. + """ + allowed_extractions = ['BOX', 'OPT'] + if self.data['extr'] not in allowed_extractions: + msgs.error("'extr' must be one of:\n" + ", ".join(allowed_extractions)) + @staticmethod def valid_algorithms(): """ diff --git a/pypeit/scripts/sensfunc.py b/pypeit/scripts/sensfunc.py index b687dab2ca..c20857bb69 100644 --- a/pypeit/scripts/sensfunc.py +++ b/pypeit/scripts/sensfunc.py @@ -21,6 +21,8 @@ def get_parser(cls, width=None): parser.add_argument("spec1dfile", type=str, help='spec1d file for the standard that will be used to compute ' 'the sensitivity function') + parser.add_argument("--extr", type=str, default='OPT', choices=['OPT', 'BOX'], + help='Extraction method to use. Default is OPT. Options are: OPT or BOX') parser.add_argument("--algorithm", type=str, default=None, choices=['UVIS', 'IR'], help="R|Override the default algorithm for computing the sensitivity " "function. Note that it is not possible to set --algorithm and " @@ -120,6 +122,15 @@ def main(args): " multi_spec_det = 3,7\n" "\n") + if args.extr is not None and args.sens_file is not None: + msgs.error("It is not possible to set --extr and simultaneously use a .sens file via " + "the --sens_file option. If you are using a .sens file set the extraction " + "method there via:\n" + "\n" + " [sensfunc]\n" + " extr = BOX\n" + "\n") + # Determine the spectrograph and generate the primary FITS header with io.fits_open(args.spec1dfile) as hdul: @@ -166,6 +177,11 @@ def main(args): multi_spec_det = [int(item) for item in args.multi.split(',')] par['sensfunc']['multi_spec_det'] = multi_spec_det + # If extr was provided override defaults. Note this does undo .sens + # file since they cannot both be passed + if args.extr is not None: + par['sensfunc']['extr'] = args.extr + # TODO Add parsing of detectors here. If detectors passed from the # command line, overwrite the parset values read in from the .sens file diff --git a/pypeit/sensfunc.py b/pypeit/sensfunc.py index 9494888054..ea783dde1c 100644 --- a/pypeit/sensfunc.py +++ b/pypeit/sensfunc.py @@ -78,6 +78,7 @@ class SensFunc(datamodel.DataContainer): 'pypeline': dict(otype=str, descr='PypeIt pipeline reduction path'), 'spec1df': dict(otype=str, descr='PypeIt spec1D file used to for sensitivity function'), + 'extr': dict(otype=str, descr='Extraction method used for the standard star (OPT or BOX)'), 'std_name': dict(otype=str, descr='Type of standard source'), 'std_cal': dict(otype=str, descr='File name (or shorthand) with the standard flux data'), @@ -213,6 +214,7 @@ def __init__(self, spec1dfile, sensfile, par, par_fluxcalib=None, debug=False, # Input and Output files self.spec1df = spec1dfile + self.extr = par['extr'] self.sensfile = sensfile self.par = par self.chk_version = chk_version @@ -253,7 +255,7 @@ def __init__(self, spec1dfile, sensfile, par, par_fluxcalib=None, debug=False, # Unpack standard wave, counts, counts_ivar, counts_mask, trace_spec, trace_spat, self.meta_spec, header \ - = self.sobjs_std.unpack_object(ret_flam=False) + = self.sobjs_std.unpack_object(ret_flam=False, extract_type=self.extr) # Compute the blaze function # TODO Make the blaze function optional @@ -471,7 +473,7 @@ def flux_std(self): # Unpack the fluxed standard _wave, _flam, _flam_ivar, _flam_mask, _, _, _, _ \ - = self.sobjs_std.unpack_object(ret_flam=True) + = self.sobjs_std.unpack_object(ret_flam=True, extract_type=self.extr) # Reshape to 2d arrays wave, flam, flam_ivar, flam_mask, _, _, _ \ = utils.spec_atleast_2d(_wave, _flam, _flam_ivar, _flam_mask) diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index 62cea7d46c..c73d303298 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -221,8 +221,10 @@ def unpack_object(self, ret_flam=False, extract_type='OPT'): flux_attr = 'FLAM' if ret_flam else 'COUNTS' flux_key = '{}_{}'.format(extract_type, flux_attr) wave_key = '{}_WAVE'.format(extract_type) - if getattr(self, flux_key)[0] is None: - msgs.error("Flux not available for {}. Try the other ".format(flux_key)) + + if np.any([f is None for f in getattr(self, flux_key)]): + other = 'OPT' if extract_type == 'BOX' else 'BOX' + msgs.error(f"Some or all fluxes are not available for {extract_type} extraction. Try {other} extraction.") # nspec = getattr(self, flux_key)[0].size # Allocate arrays and unpack spectrum @@ -1068,5 +1070,5 @@ def lst_to_array(lst, mask=None): if isinstance(lst[0], units.Quantity): return units.Quantity(lst)[mask] else: - return np.array(lst)[mask] + return np.array(lst, dtype="object")[mask] From 5f0e4e6bca5d81d8f7e232f91820bcaa2c1c2fd7 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Fri, 26 Jul 2024 15:05:32 -1000 Subject: [PATCH 03/17] fix bug --- pypeit/core/findobj_skymask.py | 2 +- pypeit/specobjs.py | 12 +++++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index fe56f70bc4..eea3a622ed 100644 --- a/pypeit/core/findobj_skymask.py +++ b/pypeit/core/findobj_skymask.py @@ -2566,7 +2566,7 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ, if len(sobjs) > 0: msgs.info('Fitting the object traces') # Note the transpose is here to pass in the TRACE_SPAT correctly. - xinit_fweight = np.copy(sobjs.TRACE_SPAT.T) + xinit_fweight = np.copy(sobjs.TRACE_SPAT.T).astype(float) spec_mask = (spec_vec >= spec_min_max_out[0]) & (spec_vec <= spec_min_max_out[1]) trc_inmask = np.outer(spec_mask, np.ones(len(sobjs), dtype=bool)) xfit_fweight = fit_trace(image, xinit_fweight, ncoeff, bpm=np.invert(inmask), maxshift=1., diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index c73d303298..ca5f7a7f81 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -1070,5 +1070,15 @@ def lst_to_array(lst, mask=None): if isinstance(lst[0], units.Quantity): return units.Quantity(lst)[mask] else: - return np.array(lst, dtype="object")[mask] + # if all the elements of lst have the same type, np.array(lst)[mask] will work + if len(set(map(type, lst))) == 1: + return np.array(lst)[mask] + else: + # if not, we need to use dtype="object" + return np.array(lst, dtype="object")[mask] + # TODO: The dtype="object" is needed for the case where one element of lst is not a list but None. + # This would prevent the error "ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions..." + # However, this is may introduce other issues, where an array of objects creates problems in other parts of the code. + # I am using this workaround. Any suggestions/ideas? + From 21f7975962645f1f74952422e0a17a238c44b0db Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Fri, 26 Jul 2024 16:25:24 -1000 Subject: [PATCH 04/17] set param to True --- pypeit/par/pypeitpar.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index a57f34bf3e..b38f25cb10 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -4189,8 +4189,7 @@ def __init__(self, trace_npoly=None, snr_thresh=None, find_trim_edge=None, 'detector. It only used for object finding. This parameter is helpful if your object only ' \ 'has emission lines or at high redshift and the trace only shows in part of the detector.' - # TODO do we want this to be True or False by default? - defaults['use_std_trace'] = False + defaults['use_std_trace'] = True dtypes['use_std_trace'] = bool descr['use_std_trace'] = 'If True, the trace of the standard star spectrum is used as a crutch for ' \ 'tracing the object spectra. This is useful when a direct trace is not possible ' \ From 7db6689de71c10a6aad17116549a954dcb59dee5 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Fri, 26 Jul 2024 17:06:16 -1000 Subject: [PATCH 05/17] small things --- doc/releases/1.16.1dev.rst | 11 +++++++++++ pypeit/spectrographs/keck_hires.py | 4 +++- pypeit/spectrographs/keck_nires.py | 2 +- 3 files changed, 15 insertions(+), 2 deletions(-) diff --git a/doc/releases/1.16.1dev.rst b/doc/releases/1.16.1dev.rst index cf5b4712c0..31133d6cc2 100644 --- a/doc/releases/1.16.1dev.rst +++ b/doc/releases/1.16.1dev.rst @@ -18,6 +18,10 @@ Dependency Changes Functionality/Performance Improvements and Additions ---------------------------------------------------- +- Added the possibility to decide if the extracted standard star spectrum should be + used as a crutch for tracing the object in the science frame (before it was done as default). + This is done by setting the parameter ``use_std_trace`` in FindObjPar. + Instrument-specific Updates --------------------------- @@ -28,6 +32,9 @@ Script Changes expansion of and changes to the cache system. - Added ``pypeit_clean_cache`` script to facilitate both viewing and removing files in the cache. +- Added the ``--extr`` parameter in the ``pypeit_sensfunc`` script (also as a SensFuncPar) + to allow the user to specify the extraction method to use when computing the sensitivity + function (before only optimal extraction was used). Datamodel Changes ----------------- @@ -46,3 +53,7 @@ Bug Fixes - Fixed a fault caused when all frames in a pypeit file are identified as being part of ``all`` calibration groups. - Allow for empty 2D wavecal solution in HDU extension of WaveCalib file +- Fix error "ValueError: setting an array element with a sequence. The requested + array has an inhomogeneous shape after 1 dimensions..." occurring when unpacking + the SpecObj spectrum but having an attribute of the SpecObj object that is None. + diff --git a/pypeit/spectrographs/keck_hires.py b/pypeit/spectrographs/keck_hires.py index ba4329f0cc..13bf1fdcaf 100644 --- a/pypeit/spectrographs/keck_hires.py +++ b/pypeit/spectrographs/keck_hires.py @@ -154,7 +154,9 @@ def default_pypeit_par(cls): par['reduce']['extraction']['model_full_slit'] = True # Mask 3 edges pixels since the slit is short, insted of default (5,5) par['reduce']['findobj']['find_trim_edge'] = [3, 3] - # Continnum order for determining thresholds + # number of objects + # par['reduce']['findobj']['maxnumber_sci'] = 2 # Assume that there is only two object in each order. + # par['reduce']['findobj']['maxnumber_std'] = 1 # Assume that there is only one object in each order. # Sensitivity function parameters par['sensfunc']['algorithm'] = 'IR' diff --git a/pypeit/spectrographs/keck_nires.py b/pypeit/spectrographs/keck_nires.py index 761fcb7253..b410a10d25 100644 --- a/pypeit/spectrographs/keck_nires.py +++ b/pypeit/spectrographs/keck_nires.py @@ -121,7 +121,7 @@ def default_pypeit_par(cls): #par['reduce']['findobj']['ech_find_nabove_min_snr'] = 1 # Require detection in a single order since given only 5 orders and slitlosses for NIRES, often # things are only detected in the K-band? Decided not to make this the default. - + par['reduce']['findobj']['maxnumber_std'] = 1 # Assume that there is only one object in each order. # Flexure par['flexure']['spec_method'] = 'skip' From 888d3726b35fcc176f99edeba1800772f4dd1d7c Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Sun, 28 Jul 2024 11:07:37 -1000 Subject: [PATCH 06/17] deal with mismatched orders --- pypeit/core/findobj_skymask.py | 54 ++++++++++++++++++++++------------ pypeit/find_objects.py | 45 ++++++++++++++++------------ pypeit/specobjs.py | 25 ++++++++++------ 3 files changed, 78 insertions(+), 46 deletions(-) diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index eea3a622ed..8b0db5fb5e 100644 --- a/pypeit/core/findobj_skymask.py +++ b/pypeit/core/findobj_skymask.py @@ -11,6 +11,7 @@ import matplotlib.pyplot as plt import astropy.stats +from astropy import table from pypeit import msgs from pypeit import utils @@ -238,10 +239,11 @@ def ech_findobj_ineach_order( Good-pixel mask for input image. Must have the same shape as ``image``. If None, all pixels in ``slitmask`` with non-negative values are considered good. - std_trace (`numpy.ndarray`_, optional): - Vector with the standard star trace, which is used as a crutch for - tracing. Shape must be (nspec,). If None, the slit boundaries are - used as the crutch. + std_trace (`astropy.table.Table`_:, optional): + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. The table has two columns: + `ECH_ORDER` and `TRACE_SPAT`. The shape of each row must be (nspec,). + If None, or for missing orders, the slit boundaries are used as the crutch. ncoeff (:obj:`int`, optional): Order of polynomial fit to traces. box_radius (:obj:`float`, optional): @@ -313,7 +315,10 @@ def ech_findobj_ineach_order( specobj_dict['SLITID'] = slit_spats[iord] specobj_dict['ECH_ORDERINDX'] = iord specobj_dict['ECH_ORDER'] = iorder - std_in = None if std_trace is None else std_trace[:, iord] + std_in = None + if std_trace is not None and 'ECH_ORDER' in std_trace.keys() and np.any(std_trace['ECH_ORDER'] == iorder): + idx = np.where(std_trace['ECH_ORDER'] == iorder)[0][0] + std_in = std_trace[idx]['TRACE_SPAT'] # Get SLTIORD_ID for the objfind QA ech_objfindQA_filename = objfindQA_filename.replace('S0999', 'S{:04d}'.format(order_vec[iord])) \ @@ -444,7 +449,7 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs, slit_spat_id: np.ndarray, order_vec:np.ndarray, obj_id:np.ndarray, - std_trace:specobjs.SpecObjs=None, + std_trace:table.Table=None, show:bool=False): """ For objects which were only found on some orders, the standard (or @@ -484,9 +489,11 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs, ``np.arange(norders)`` (but this is *not* recommended). obj_id (`numpy.ndarray`_): Object IDs of the objects linked together. - std_trace (:class:`~pypeit.specobjs.SpecObjs`, optional): - Standard star objects (including the traces) - Defaults to None. + std_trace (`astropy.table.Table`_:, optional): + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. The table has two columns: + `ECH_ORDER` and `TRACE_SPAT`. The shape of each row must be (nspec,). + If None, or for missing orders, the slit boundaries are used as the crutch. show (bool, optional): Plot diagnostics related to filling the missing orders @@ -506,8 +513,9 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs, slit_width = slit_righ - slit_left # Check standard star - if std_trace is not None and std_trace.shape[1] != norders: - msgs.error('Standard star trace does not match the number of orders in the echelle data.') + if std_trace is not None and len(std_trace) != norders: + msgs.warn('Standard star trace does not match the number of orders in the echelle data.' + ' Will use the slit edges to trace the object in the missing orders.') # For traces nspec = slit_left.shape[0] @@ -607,12 +615,16 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs, #thisobj.ech_order = order_vec[iord] thisobj.SPAT_FRACPOS = uni_frac[iobj] # Assign traces using the fractional position fit above - if std_trace is not None: - x_trace = np.interp(slit_spec_pos, spec_vec, std_trace[:,iord]) + if std_trace is not None and 'ECH_ORDER' in std_trace.keys() and \ + np.any(std_trace['ECH_ORDER'] == this_order): + idx = np.where(std_trace['ECH_ORDER'] == this_order)[0][0] + # standard star trace in this order + std_in = std_trace[idx]['TRACE_SPAT'] + x_trace = np.interp(slit_spec_pos, spec_vec, std_in[:,iord]) shift = np.interp( slit_spec_pos, spec_vec, slit_left[:,iord] + slit_width[:,iord]*frac_mean_new[iord]) - x_trace - thisobj.TRACE_SPAT = std_trace[:,iord] + shift + thisobj.TRACE_SPAT = std_in[:,iord] + shift else: thisobj.TRACE_SPAT = slit_left[:, iord] + slit_width[:, iord] * frac_mean_new[iord] # new trace thisobj.trace_spec = spec_vec @@ -1124,10 +1136,11 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order Plate scale in arcsec/pix. This can either be a single float for every order, or an array with shape (norders,) providing the plate scale of each order. - std_trace (`numpy.ndarray`_, optional): - Vector with the standard star trace, which is used as a crutch for - tracing. Shape must be (nspec,). If None, the slit boundaries are - used as the crutch. + std_trace (`astropy.table.Table`_:, optional): + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. The table has two columns: + `ECH_ORDER` and `TRACE_SPAT`. The shape of each row must be (nspec,). + If None, or for missing orders, the slit boundaries are used as the crutch. ncoeff (:obj:`int`, optional): Order of polynomial fit to traces. npca (:obj:`int`, optional): @@ -1325,7 +1338,7 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order return sobjs_ech - +# TODO: THIS IS DEPRECATED def orig_ech_objfind(image, ivar, slitmask, slit_left, slit_righ, order_vec, maskslits, det='DET01', inmask=None, spec_min_max=None, fof_link=1.5, plate_scale=0.2, std_trace=None, ncoeff=5, npca=None, coeff_npoly=None, max_snr=2.0, min_snr=1.0, @@ -2537,6 +2550,9 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ, # If no standard trace is provided shift left slit boundary over to be initial trace else: # ToDO make this the average left and right boundary instead. That would be more robust. + # Print a status message for the first object + if iobj == 0: + msgs.info('Using slit edges as crutch for object tracing') sobjs[iobj].TRACE_SPAT = slit_left + xsize*sobjs[iobj].SPAT_FRACPOS sobjs[iobj].trace_spec = spec_vec diff --git a/pypeit/find_objects.py b/pypeit/find_objects.py index 6ca9dbbeed..e95284e93a 100644 --- a/pypeit/find_objects.py +++ b/pypeit/find_objects.py @@ -326,8 +326,13 @@ def run(self, std_trace=None, show_peaks=False, show_skysub_fit=False): Parameters ---------- - std_trace : `numpy.ndarray`_, optional - Trace of the standard star + std_trace : `astropy.table.Table`_:, optional + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. For MultiSlit reduction, + the table has a single column: `TRACE_SPAT`. + For Echelle reduction, the table has two columns: `ECH_ORDER` and `TRACE_SPAT`. + The shape of each row must be (nspec,). For SlicerIFU reduction, std_trace is None. + If None, the slit boundaries are used as the crutch. show_peaks : :obj:`bool`, optional Show peaks in find_objects methods show_skysub_fit : :obj:`bool`, optional @@ -400,11 +405,13 @@ def find_objects(self, image, ivar, std_trace=None, Image to search for objects from. This floating-point image has shape (nspec, nspat) where the first dimension (nspec) is spectral, and second dimension (nspat) is spatial. - std_trace : `numpy.ndarray`_, optional - This is a one dimensional float array with shape = (nspec,) containing the standard star - trace which is used as a crutch for tracing. If the no - standard star is provided the code uses the the slit - boundaries as the crutch. + std_trace : `astropy.table.Table`_:, optional + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. For MultiSlit reduction, + the table has a single column: `TRACE_SPAT`. + For Echelle reduction, the table has two columns: `ECH_ORDER` and `TRACE_SPAT`. + The shape of each row must be (nspec,). For SlicerIFU reduction, std_trace is None. + If None, the slit boundaries are used as the crutch. show_peaks : :obj:`bool`, optional Generate QA showing peaks identified by object finding show_fits : :obj:`bool`, optional @@ -724,11 +731,12 @@ def find_objects_pypeline(self, image, ivar, std_trace=None, Image to search for objects from. This floating-point image has shape (nspec, nspat) where the first dimension (nspec) is spectral, and second dimension (nspat) is spatial. - std_trace : `numpy.ndarray`_, optional - This is a one dimensional float array with shape = (nspec,) containing the standard star - trace which is used as a crutch for tracing. If the no - standard star is provided the code uses the the slit - boundaries as the crutch. + std_trace : `astropy.table.Table`_:, optional + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. For MultiSlit reduction, + the table has a single column: `TRACE_SPAT`. + The shape of each row must be (nspec,). If None, + the slit boundaries are used as the crutch. manual_extract_dict : :obj:`dict`, optional Dict guiding the manual extraction show_peaks : :obj:`bool`, optional @@ -798,7 +806,7 @@ def find_objects_pypeline(self, image, ivar, std_trace=None, self.slits_right[:,slit_idx], inmask=inmask, ncoeff=self.par['reduce']['findobj']['trace_npoly'], - std_trace=std_trace, + std_trace=std_trace[0]['TRACE_SPAT'], snr_thresh=snr_thresh, hand_extract_dict=manual_extract_dict, specobj_dict=specobj_dict, show_peaks=show_peaks, @@ -874,11 +882,12 @@ def find_objects_pypeline(self, image, ivar, std_trace=None, Image to search for objects from. This floating-point image has shape (nspec, nspat) where the first dimension (nspec) is spectral, and second dimension (nspat) is spatial. - std_trace : `numpy.ndarray`_, optional - This is a one dimensional float array with shape = (nspec,) containing the standard star - trace which is used as a crutch for tracing. If the no - standard star is provided the code uses the the slit - boundaries as the crutch. + std_trace : `astropy.table.Table`_:, optional + Table with the trace of the standard star on the input detector, + which is used as a crutch for tracing. For Echelle reduction, + the table has two columns: `ECH_ORDER` and `TRACE_SPAT`. + The shape of each row must be (nspec,). If None, + the slit boundaries are used as the crutch. manual_extract_dict : :obj:`dict`, optional Dict guiding the manual extraction show_peaks : :obj:`bool`, optional diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index ca5f7a7f81..4e33c12c81 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -1004,6 +1004,7 @@ def get_extraction_groups(self, model_full_slit=False) -> List[List[int]]: return groups + #TODO Should this be a classmethod on specobjs?? def get_std_trace(detname, std_outfile, chk_version=True): """ @@ -1016,9 +1017,11 @@ def get_std_trace(detname, std_outfile, chk_version=True): Filename with the standard star spec1d file. Can be None. Returns: - `numpy.ndarray`_: Trace of the standard star on input detector. Will - be None if ``std_outfile`` is None, or if the selected detector/mosaic - is not available in the provided spec1d file. + `astropy.table.Table`_: Table with the trace of the standard star on the input detector. + If this is a MultiSlit reduction, the table will have a single column: `TRACE_SPAT`. + If this is an Echelle reduction, the table will have two columns: `ECH_ORDER` and `TRACE_SPAT`. + Will be None if ``std_outfile`` is None, or if the selected detector/mosaic + is not available in the provided spec1d file, or for SlicerIFU reductions. """ sobjs = SpecObjs.from_fitsfile(std_outfile, chk_version=chk_version) @@ -1035,20 +1038,24 @@ def get_std_trace(detname, std_outfile, chk_version=True): # No standard extracted on this detector?? if sobjs_std is None: return None - std_trace = sobjs_std.TRACE_SPAT + + # create table that contains the trace of the standard + std_tab = Table() # flatten the array if this multislit if 'MultiSlit' in pypeline: - std_trace = std_trace.flatten() + std_tab['TRACE_SPAT'] = sobjs_std.TRACE_SPAT elif 'Echelle' in pypeline: - std_trace = std_trace.T + std_tab['ECH_ORDER'] = sobjs_std.ECH_ORDER + std_tab['TRACE_SPAT'] = sobjs_std.TRACE_SPAT elif 'SlicerIFU' in pypeline: - std_trace = None + std_tab = None else: msgs.error('Unrecognized pypeline') else: - std_trace = None + std_tab = None + + return std_tab - return std_trace def lst_to_array(lst, mask=None): """ From 707db079b066b01d29a50d9b6ee3d2183a8d3306 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Sun, 28 Jul 2024 11:10:44 -1000 Subject: [PATCH 07/17] add to release --- doc/releases/1.16.1dev.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/releases/1.16.1dev.rst b/doc/releases/1.16.1dev.rst index 31133d6cc2..86e1279582 100644 --- a/doc/releases/1.16.1dev.rst +++ b/doc/releases/1.16.1dev.rst @@ -21,6 +21,7 @@ Functionality/Performance Improvements and Additions - Added the possibility to decide if the extracted standard star spectrum should be used as a crutch for tracing the object in the science frame (before it was done as default). This is done by setting the parameter ``use_std_trace`` in FindObjPar. +- Now PypeIt can handle the case where "Standard star trace does not match the number of orders in the echelle data" Instrument-specific Updates --------------------------- From dcb8b3eaa0ef0282149d19261314864961002787 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Sun, 28 Jul 2024 11:23:35 -1000 Subject: [PATCH 08/17] fix bug --- pypeit/core/findobj_skymask.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index 8b0db5fb5e..58078c9cc9 100644 --- a/pypeit/core/findobj_skymask.py +++ b/pypeit/core/findobj_skymask.py @@ -620,11 +620,11 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs, idx = np.where(std_trace['ECH_ORDER'] == this_order)[0][0] # standard star trace in this order std_in = std_trace[idx]['TRACE_SPAT'] - x_trace = np.interp(slit_spec_pos, spec_vec, std_in[:,iord]) + x_trace = np.interp(slit_spec_pos, spec_vec, std_in) shift = np.interp( slit_spec_pos, spec_vec, slit_left[:,iord] + slit_width[:,iord]*frac_mean_new[iord]) - x_trace - thisobj.TRACE_SPAT = std_in[:,iord] + shift + thisobj.TRACE_SPAT = std_in + shift else: thisobj.TRACE_SPAT = slit_left[:, iord] + slit_width[:, iord] * frac_mean_new[iord] # new trace thisobj.trace_spec = spec_vec From fce85e6067c722b5d9e2e777f44c2540156afb9c Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Mon, 29 Jul 2024 08:11:58 -1000 Subject: [PATCH 09/17] fix bug --- pypeit/find_objects.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pypeit/find_objects.py b/pypeit/find_objects.py index a8975efc7b..dbb74596a8 100644 --- a/pypeit/find_objects.py +++ b/pypeit/find_objects.py @@ -800,13 +800,17 @@ def find_objects_pypeline(self, image, ivar, std_trace=None, maxnumber = self.par['reduce']['findobj']['maxnumber_std'] if self.std_redux \ else self.par['reduce']['findobj']['maxnumber_sci'] + # standard star + std_in = std_trace[0]['TRACE_SPAT'] if std_trace is not None else None + + # Find objects sobjs_slit = \ findobj_skymask.objs_in_slit(image, ivar, thismask, self.slits_left[:,slit_idx], self.slits_right[:,slit_idx], inmask=inmask, ncoeff=self.par['reduce']['findobj']['trace_npoly'], - std_trace=std_trace[0]['TRACE_SPAT'], + std_trace=std_in, snr_thresh=snr_thresh, hand_extract_dict=manual_extract_dict, specobj_dict=specobj_dict, show_peaks=show_peaks, From 64ca0283f0c17e268015461d3a577d1ed1de10f1 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Mon, 29 Jul 2024 19:36:55 -1000 Subject: [PATCH 10/17] fix mismatch order sensfunc data --- pypeit/coadd1d.py | 13 ++++--------- pypeit/sensfunc.py | 35 ++++++++++++++++++++++++++++------- 2 files changed, 32 insertions(+), 16 deletions(-) diff --git a/pypeit/coadd1d.py b/pypeit/coadd1d.py index fc6eafcbbf..d33dfe0cd4 100644 --- a/pypeit/coadd1d.py +++ b/pypeit/coadd1d.py @@ -202,7 +202,6 @@ def __init__(self, spec1dfiles, objids, spectrograph=None, par=None, sensfuncfil super().__init__(spec1dfiles, objids, spectrograph=spectrograph, par=par, sensfuncfile=sensfuncfile, setup_id=setup_id, debug=debug, show=show, chk_version=chk_version) - def load(self): """ Load the arrays we need for performing coadds. @@ -461,8 +460,7 @@ def coadd(self): lower=self.par['lower'], upper=self.par['upper'], maxrej=self.par['maxrej'], sn_clip=self.par['sn_clip'], debug=self.debug, show=self.show, show_exp=self.show) - - + return wave_grid_mid, wave_coadd, flux_coadd, ivar_coadd, gpm_coadd, order_stacks def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles): @@ -474,7 +472,7 @@ def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles): List of spec1d files for this setup. objids (list): List of objids. This is aligned with spec1dfiles - sensfuncfile (list): + sensfuncfiles (list): List of sensfuncfiles. This is aligned with spec1dfiles and objids Returns: @@ -488,7 +486,7 @@ def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles): indx = sobjs.name_indices(objids[iexp]) if not np.any(indx): msgs.error("No matching objects for {:s}. Odds are you input the wrong OBJID".format(objids[iexp])) - wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, _, _, _, header = \ + wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, _, _, meta_spec, header = \ sobjs[indx].unpack_object(ret_flam=self.par['flux_value'], extract_type=self.par['ex_value']) # This np.atleast2d hack deals with the situation where we are wave_iexp is actually Multislit data, i.e. we are treating # it like an echelle spectrograph with a single order. This usage case arises when we want to use the @@ -496,6 +494,7 @@ def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles): if wave_iexp.ndim == 1: wave_iexp, flux_iexp, ivar_iexp, gpm_iexp = np.atleast_2d(wave_iexp).T, np.atleast_2d(flux_iexp).T, np.atleast_2d(ivar_iexp).T, np.atleast_2d(gpm_iexp).T weights_sens_iexp = sensfunc.SensFunc.sensfunc_weights(sensfuncfiles[iexp], wave_iexp, + ech_order_vec=meta_spec['ECH_ORDERS'], debug=self.debug, chk_version=self.chk_version) # Allocate arrays on first iteration @@ -515,10 +514,8 @@ def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles): waves[...,iexp], fluxes[...,iexp], ivars[..., iexp], gpms[...,iexp], weights_sens[...,iexp] \ = wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, weights_sens_iexp - return waves, fluxes, ivars, gpms, weights_sens, header_out - def load(self): """ Load the arrays we need for performing echelle coadds. @@ -563,8 +560,6 @@ def load(self): for c, l in zip(combined, loaded): c.append(l) - - return waves, fluxes, ivars, gpms, weights_sens, headers diff --git a/pypeit/sensfunc.py b/pypeit/sensfunc.py index 1b5d548cb2..fdd7f4283c 100644 --- a/pypeit/sensfunc.py +++ b/pypeit/sensfunc.py @@ -807,7 +807,7 @@ def write_QA(self): @classmethod - def sensfunc_weights(cls, sensfile, waves, debug=False, extrap_sens=True, chk_version=True): + def sensfunc_weights(cls, sensfile, waves, ech_order_vec=None, debug=False, extrap_sens=True, chk_version=True): """ Get the weights based on the sensfunc @@ -817,8 +817,12 @@ def sensfunc_weights(cls, sensfile, waves, debug=False, extrap_sens=True, chk_ve waves (`numpy.ndarray`_): wavelength grid for your output weights. Shape is (nspec, norders, nexp) or (nspec, norders). + ech_order_vec (`numpy.ndarray`_, optional): + Vector of echelle orders. Only used for echelle data. debug (bool): default=False show the weights QA + extrap_sens (bool): default=True + Extrapolate the sensitivity function chk_version (:obj:`bool`, optional): When reading in existing files written by PypeIt, perform strict version checking to ensure a valid file. If False, the code @@ -833,6 +837,10 @@ def sensfunc_weights(cls, sensfile, waves, debug=False, extrap_sens=True, chk_ve if waves.ndim == 2: nspec, norder = waves.shape + if ech_order_vec.size != norder: + msgs.warn('The number of orders in the wave grid does not match the ' + 'number of orders in the unpacked sobjs. Echelle order vector not used.') + ech_order_vec = None nexp = 1 waves_stack = np.reshape(waves, (nspec, norder, 1)) elif waves.ndim == 3: @@ -845,16 +853,29 @@ def sensfunc_weights(cls, sensfile, waves, debug=False, extrap_sens=True, chk_ve else: msgs.error('Unrecognized dimensionality for waves') - weights_stack = np.zeros_like(waves_stack) + weights_stack = np.ones_like(waves_stack) - if norder != sens.zeropoint.shape[1]: + if norder != sens.zeropoint.shape[1] and ech_order_vec is None: msgs.error('The number of orders in {:} does not agree with your data. Wrong sensfile?'.format(sensfile)) - - for iord in range(norder): + elif norder != sens.zeropoint.shape[1] and ech_order_vec is not None: + msgs.warn('The number of orders in {:} does not match the number of orders in the data. ' + 'Using only the matching orders.'.format(sensfile)) + + # array of order to loop through + orders = np.arange(norder) if ech_order_vec is None else ech_order_vec + for iord,this_ord in enumerate(orders): + if ech_order_vec is None: + isens = iord + # find the index of the sensfunc for this order + elif ech_order_vec is not None and np.any(sens.sens['ECH_ORDERS'].value == this_ord): + isens = np.where(sens.sens['ECH_ORDERS'].value == this_ord)[0][0] + else: + # if the order is not in the sensfunc file, skip it + continue for iexp in range(nexp): sensfunc_iord = flux_calib.get_sensfunc_factor(waves_stack[:,iord,iexp], - sens.wave[:,iord], - sens.zeropoint[:,iord], 1.0, + sens.wave[:,isens], + sens.zeropoint[:,isens], 1.0, extrap_sens=extrap_sens) weights_stack[:,iord,iexp] = utils.inverse(sensfunc_iord) From ff3b4660cdfd0397f3366249d26cecad989e6995 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Mon, 29 Jul 2024 19:55:22 -1000 Subject: [PATCH 11/17] error and doc --- doc/releases/1.16.1dev.rst | 8 ++++---- pypeit/coadd1d.py | 12 +++++++++--- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/doc/releases/1.16.1dev.rst b/doc/releases/1.16.1dev.rst index 304a056372..dacc94368e 100644 --- a/doc/releases/1.16.1dev.rst +++ b/doc/releases/1.16.1dev.rst @@ -22,6 +22,7 @@ Functionality/Performance Improvements and Additions used as a crutch for tracing the object in the science frame (before it was done as default). This is done by setting the parameter ``use_std_trace`` in FindObjPar. - Now PypeIt can handle the case where "Standard star trace does not match the number of orders in the echelle data" + both in `run_pypeit` and in `pypeit_coadd_1dspec`. Instrument-specific Updates --------------------------- @@ -33,15 +34,14 @@ Script Changes expansion of and changes to the cache system. - Added ``pypeit_clean_cache`` script to facilitate both viewing and removing files in the cache. -- Added the ``--extr`` parameter in the ``pypeit_sensfunc`` script (also as a SensFuncPar) - to allow the user to specify the extraction method to use when computing the sensitivity - function (before only optimal extraction was used). - - A new script, called `pypeit_extract_datacube`, allows 1D spectra of point sources to be extracted from datacubes. - The sensitivity function is now generated outside of datacube generation. - The `grating_corr` column is now used to select the correct grating correction file for each spec2d file when generating the datacube. +- Added the ``--extr`` parameter in the ``pypeit_sensfunc`` script (also as a SensFuncPar) + to allow the user to specify the extraction method to use when computing the sensitivity + function (before only optimal extraction was used). Datamodel Changes ----------------- diff --git a/pypeit/coadd1d.py b/pypeit/coadd1d.py index d33dfe0cd4..e7ef072031 100644 --- a/pypeit/coadd1d.py +++ b/pypeit/coadd1d.py @@ -510,9 +510,15 @@ def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles): header_out['RA_OBJ'] = sobjs[indx][0]['RA'] header_out['DEC_OBJ'] = sobjs[indx][0]['DEC'] - # Store the information - waves[...,iexp], fluxes[...,iexp], ivars[..., iexp], gpms[...,iexp], weights_sens[...,iexp] \ - = wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, weights_sens_iexp + # TODO :: The error below can be removed if we refactor to use a list of numpy arrays. But, if we do that, + # we need to make several changes to the ech_combspec function. + try: + # Store the information + waves[...,iexp], fluxes[...,iexp], ivars[..., iexp], gpms[...,iexp], weights_sens[...,iexp] \ + = wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, weights_sens_iexp + except ValueError: + msgs.error('The shape (Nspec,Norder) of spectra is not consistent between exposures. ' + 'These spec1ds cannot be coadded at this time.') return waves, fluxes, ivars, gpms, weights_sens, header_out From 27410e4b5b1715be4c41ae23923cd371524d723e Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 31 Jul 2024 11:27:43 -1000 Subject: [PATCH 12/17] fix a bug --- pypeit/scripts/sensfunc.py | 10 ++++++++-- pypeit/spectrographs/keck_hires.py | 4 ++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/pypeit/scripts/sensfunc.py b/pypeit/scripts/sensfunc.py index 5477b84270..53ab98a3e0 100644 --- a/pypeit/scripts/sensfunc.py +++ b/pypeit/scripts/sensfunc.py @@ -21,8 +21,14 @@ def get_parser(cls, width=None): parser.add_argument("spec1dfile", type=str, help='spec1d file for the standard that will be used to compute ' 'the sensitivity function') - parser.add_argument("--extr", type=str, default='OPT', choices=['OPT', 'BOX'], - help='Extraction method to use. Default is OPT. Options are: OPT or BOX') + parser.add_argument("--extr", type=str, default=None, choices=['OPT', 'BOX'], + help="R|Override the default extraction method used for computing the sensitivity " + "function. Note that it is not possible to set --extr and " + "simultaneously use a .sens file with the --sens_file option. If " + "you are using a .sens file, set the algorithm there via:\n\n" + "F| [sensfunc]\n" + "F| extr = BOX\n" + "\nThe extraction options are: OPT or BOX") parser.add_argument("--algorithm", type=str, default=None, choices=['UVIS', 'IR'], help="R|Override the default algorithm for computing the sensitivity " "function. Note that it is not possible to set --algorithm and " diff --git a/pypeit/spectrographs/keck_hires.py b/pypeit/spectrographs/keck_hires.py index 13bf1fdcaf..22c92e639c 100644 --- a/pypeit/spectrographs/keck_hires.py +++ b/pypeit/spectrographs/keck_hires.py @@ -155,8 +155,8 @@ def default_pypeit_par(cls): # Mask 3 edges pixels since the slit is short, insted of default (5,5) par['reduce']['findobj']['find_trim_edge'] = [3, 3] # number of objects - # par['reduce']['findobj']['maxnumber_sci'] = 2 # Assume that there is only two object in each order. - # par['reduce']['findobj']['maxnumber_std'] = 1 # Assume that there is only one object in each order. + par['reduce']['findobj']['maxnumber_sci'] = 2 # Assume that there is max two object in each order. + par['reduce']['findobj']['maxnumber_std'] = 1 # Assume that there is only one object in each order. # Sensitivity function parameters par['sensfunc']['algorithm'] = 'IR' From 717bd6565af3c07699ab97b875023f07ff956c18 Mon Sep 17 00:00:00 2001 From: Kyle Westfall Date: Tue, 13 Aug 2024 15:57:13 -0700 Subject: [PATCH 13/17] PR comments --- pypeit/core/findobj_skymask.py | 6 +++--- pypeit/find_objects.py | 9 ++++----- pypeit/specobjs.py | 34 +++++++++++++++++++--------------- 3 files changed, 26 insertions(+), 23 deletions(-) diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index 58078c9cc9..f96f58e667 100644 --- a/pypeit/core/findobj_skymask.py +++ b/pypeit/core/findobj_skymask.py @@ -239,7 +239,7 @@ def ech_findobj_ineach_order( Good-pixel mask for input image. Must have the same shape as ``image``. If None, all pixels in ``slitmask`` with non-negative values are considered good. - std_trace (`astropy.table.Table`_:, optional): + std_trace (`astropy.table.Table`_, optional): Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. The table has two columns: `ECH_ORDER` and `TRACE_SPAT`. The shape of each row must be (nspec,). @@ -489,7 +489,7 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs, ``np.arange(norders)`` (but this is *not* recommended). obj_id (`numpy.ndarray`_): Object IDs of the objects linked together. - std_trace (`astropy.table.Table`_:, optional): + std_trace (`astropy.table.Table`_, optional): Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. The table has two columns: `ECH_ORDER` and `TRACE_SPAT`. The shape of each row must be (nspec,). @@ -1136,7 +1136,7 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order Plate scale in arcsec/pix. This can either be a single float for every order, or an array with shape (norders,) providing the plate scale of each order. - std_trace (`astropy.table.Table`_:, optional): + std_trace (`astropy.table.Table`_, optional): Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. The table has two columns: `ECH_ORDER` and `TRACE_SPAT`. The shape of each row must be (nspec,). diff --git a/pypeit/find_objects.py b/pypeit/find_objects.py index dbb74596a8..44ba731508 100644 --- a/pypeit/find_objects.py +++ b/pypeit/find_objects.py @@ -15,7 +15,6 @@ from pypeit import specobjs from pypeit import msgs, utils -from pypeit import flatfield from pypeit.display import display from pypeit.core import skysub, qa, parse, flat, flexure from pypeit.core import procimg @@ -326,7 +325,7 @@ def run(self, std_trace=None, show_peaks=False, show_skysub_fit=False): Parameters ---------- - std_trace : `astropy.table.Table`_:, optional + std_trace : `astropy.table.Table`_, optional Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. For MultiSlit reduction, the table has a single column: `TRACE_SPAT`. @@ -405,7 +404,7 @@ def find_objects(self, image, ivar, std_trace=None, Image to search for objects from. This floating-point image has shape (nspec, nspat) where the first dimension (nspec) is spectral, and second dimension (nspat) is spatial. - std_trace : `astropy.table.Table`_:, optional + std_trace : `astropy.table.Table`_, optional Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. For MultiSlit reduction, the table has a single column: `TRACE_SPAT`. @@ -731,7 +730,7 @@ def find_objects_pypeline(self, image, ivar, std_trace=None, Image to search for objects from. This floating-point image has shape (nspec, nspat) where the first dimension (nspec) is spectral, and second dimension (nspat) is spatial. - std_trace : `astropy.table.Table`_:, optional + std_trace : `astropy.table.Table`_, optional Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. For MultiSlit reduction, the table has a single column: `TRACE_SPAT`. @@ -886,7 +885,7 @@ def find_objects_pypeline(self, image, ivar, std_trace=None, Image to search for objects from. This floating-point image has shape (nspec, nspat) where the first dimension (nspec) is spectral, and second dimension (nspat) is spatial. - std_trace : `astropy.table.Table`_:, optional + std_trace : `astropy.table.Table`_, optional Table with the trace of the standard star on the input detector, which is used as a crutch for tracing. For Echelle reduction, the table has two columns: `ECH_ORDER` and `TRACE_SPAT`. diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index ad82265f9c..deff44be53 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -1070,7 +1070,7 @@ def lst_to_array(lst, mask=None): Allows for a list of Quantity objects Args: - lst : list + lst (:obj:`list`): Should be number or Quantities mask (`numpy.ndarray`_, optional): Boolean array used to limit to a subset of the list. True=good @@ -1078,20 +1078,24 @@ def lst_to_array(lst, mask=None): Returns: `numpy.ndarray`_, `astropy.units.Quantity`_: Converted list """ - if mask is None: - mask = np.array([True]*len(lst)) + _mask = np.ones(len(lst), dtype=bool) if mask is None else mask + + # Return a Quantity array if isinstance(lst[0], units.Quantity): - return units.Quantity(lst)[mask] - else: - # if all the elements of lst have the same type, np.array(lst)[mask] will work - if len(set(map(type, lst))) == 1: - return np.array(lst)[mask] - else: - # if not, we need to use dtype="object" - return np.array(lst, dtype="object")[mask] - # TODO: The dtype="object" is needed for the case where one element of lst is not a list but None. - # This would prevent the error "ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions..." - # However, this is may introduce other issues, where an array of objects creates problems in other parts of the code. - # I am using this workaround. Any suggestions/ideas? + return units.Quantity(lst)[_mask] + + # If all the elements of lst have the same type, np.array(lst)[mask] will work + if len(set(map(type, lst))) == 1: + return np.array(lst)[_mask] + + # Otherwise, we have to set the array type to object + return np.array(lst, dtype=object)[_mask] + + # TODO: The dtype="object" is needed for the case where one element of lst + # is not a list but None. This would prevent the error "ValueError: setting + # an array element with a sequence. The requested array has an inhomogeneous + # shape after 1 dimensions..." However, this is may introduce other issues, + # where an array of objects creates problems in other parts of the code. I + # am using this workaround. Any suggestions/ideas? From c26be35643f8cb75e39207deedead39b7f53f9cf Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Tue, 13 Aug 2024 14:29:22 -1000 Subject: [PATCH 14/17] move orig_ech_objfind to deprecated --- pypeit/core/findobj_skymask.py | 675 --------------------------------- 1 file changed, 675 deletions(-) diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index 9a6d34e889..e2d624a37f 100644 --- a/pypeit/core/findobj_skymask.py +++ b/pypeit/core/findobj_skymask.py @@ -1338,681 +1338,6 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order return sobjs_ech -# TODO: THIS IS DEPRECATED -def orig_ech_objfind(image, ivar, slitmask, slit_left, slit_righ, order_vec, maskslits, det='DET01', - inmask=None, spec_min_max=None, fof_link=1.5, plate_scale=0.2, - std_trace=None, ncoeff=5, npca=None, coeff_npoly=None, max_snr=2.0, min_snr=1.0, - nabove_min_snr=2, pca_explained_var=99.0, box_radius=2.0, fwhm=3.0, - use_user_fwhm=False, maxdev=2.0, hand_extract_dict=None, nperorder=2, - extract_maskwidth=3.0, snr_thresh=10.0, - specobj_dict=None, trim_edg=(5,5), - show_peaks=False, show_fits=False, show_single_fits=False, - show_trace=False, show_single_trace=False, show_pca=False, - debug_all=False, objfindQA_filename=None): - """ - Object finding routine for Echelle spectrographs. - - This routine: - - #. Runs object finding on each order individually - - #. Links the objects found together using a friends-of-friends algorithm - on fractional order position. - - #. For objects which were only found on some orders, the standard (or - the slit boundaries) are placed at the appropriate fractional - position along the order. - - #. A PCA fit to the traces is performed using the routine above pca_fit - - Args: - image (`numpy.ndarray`_): - (Floating-point) Image to use for object search with shape (nspec, - nspat). The first dimension (nspec) is spectral, and second - dimension (nspat) is spatial. Note this image can either have the - sky background in it, or have already been sky subtracted. Object - finding works best on sky-subtracted images. Ideally, object finding - is run in another routine, global sky-subtraction performed, and - then this code should be run. However, it is also possible to run - this code on non-sky-subtracted images. - ivar (`numpy.ndarray`_): - Floating-point inverse variance image for the input image. Shape - must match ``image``, (nspec, nspat). - slitmask (`numpy.ndarray`_): - Integer image indicating the pixels that belong to each order. - Pixels that are not on an order have value -1, and those that are on - an order have a value equal to the slit number (i.e. 0 to nslits-1 - from left to right on the image). Shape must match ``image``, - (nspec, nspat). - slit_left (`numpy.ndarray`_): - Left boundary of orders to be extracted (given as floating point - pixels). Shape is (nspec, norders), where norders is the total - number of traced echelle orders. - slit_righ (`numpy.ndarray`_): - Right boundary of orders to be extracted (given as floating point - pixels). Shape is (nspec, norders), where norders is the total - number of traced echelle orders. - order_vec (`numpy.ndarray`_): - Vector identifying the Echelle orders for each pair of order edges - found. This is saved to the output :class:`~pypeit.specobj.SpecObj` - objects. If the orders are not known, this can be - ``np.arange(norders)`` (but this is *not* recommended). - maskslits (`numpy.ndarray`_): - Boolean array selecting orders that should be ignored (i.e., good - orders are False, bad orders are True). Shape must be (norders,). - det (:obj:`str`, optional): - The name of the detector containing the object. Only used if - ``specobj_dict`` is None. - inmask (`numpy.ndarray`_, optional): - Good-pixel mask for input image. Must have the same shape as - ``image``. If None, all pixels in ``slitmask`` with non-negative - values are considered good. - spec_min_max (`numpy.ndarray`_, optional): - 2D array defining the minimum and maximum pixel in the spectral - direction with useable data for each order. Shape must be (2, - norders). This should only be used for echelle spectrographs for - which the orders do not entirely cover the detector. PCA tracing - will re-map the traces such that they all have the same length, - compute the PCA, and then re-map the orders back. This improves - performance for echelle spectrographs by removing the nonlinear - shrinking of the orders so that the linear pca operation can better - predict the traces. If None, the minimum and maximum values will be - determined automatically from ``slitmask``. - fof_link (:obj:`float`, optional): - Friends-of-friends linking length in arcseconds used to link - together traces across orders. The routine links together at - the same fractional slit position and links them together - with a friends-of-friends algorithm using this linking - length. - plate_scale (:obj:`float`, `numpy.ndarray`_, optional): - Plate scale in arcsec/pix. This can either be a single float for - every order, or an array with shape (norders,) providing the plate - scale of each order. - std_trace (`numpy.ndarray`_, optional): - Vector with the standard star trace, which is used as a crutch for - tracing. Shape must be (nspec,). If None, the slit boundaries are - used as the crutch. - ncoeff (:obj:`int`, optional): - Order of polynomial fit to traces. - npca (:obj:`int`, optional): - Number of PCA components to keep during PCA decomposition of the - object traces. If None, the number of components set by requiring - the PCA accounts for approximately 99% of the variance. - coeff_npoly (:obj:`int`, optional): - Order of polynomial used for PCA coefficients fitting. If None, - value set automatically, see - :func:`~pypeit.tracepca.pca_trace_object`. - max_snr (:obj:`float`, optional): - For an object to be included in the output object, it must have a - max S/N ratio above this value. - min_snr (:obj:`float`, optional): - For an object to be included in the output object, it must have a - a median S/N ratio above this value for at least - ``nabove_min_snr`` orders (see below). - nabove_min_snr (:obj:`int`, optional): - The required number of orders that an object must have with median - SNR greater than ``min_snr`` in order to be included in the output - object. - pca_explained_var (:obj:`float`, optional): - The percentage (i.e., not the fraction) of the variance in the data - accounted for by the PCA used to truncate the number of PCA - coefficients to keep (see ``npca``). Ignored if ``npca`` is provided - directly; see :func:`~pypeit.tracepca.pca_trace_object`. - box_radius (:obj:`float`, optional): - Box_car extraction radius in arcseconds to assign to each detected - object and to be used later for boxcar extraction. In this method - ``box_radius`` is converted into pixels using ``plate_scale``. - ``box_radius`` is also used for SNR calculation and trimming. - fwhm (:obj:`float`, optional): - Estimated fwhm of the objects in pixels - use_user_fwhm (:obj:`bool`, optional): - If True, ``PypeIt`` will use the spatial profile FWHM input by the - user (see ``fwhm``) rather than determine the spatial FWHM from the - smashed spatial profile via the automated algorithm. - maxdev (:obj:`float`, optional): - Maximum deviation of pixels from polynomial fit to trace - used to reject bad pixels in trace fitting. - hand_extract_dict (:obj:`dict`, optional): - Dictionary with info on manual extraction; see - :class:`~pypeit.manual_extract.ManualExtractionObj`. - nperorder (:obj:`int`, optional): - Maximum number of objects allowed per order. If there are more - detections than this number, the code will select the ``nperorder`` - most significant detections. However, hand apertures will always be - returned and do not count against this budget. - extract_maskwidth (:obj:`float`, optional): - Determines the initial size of the region in units of FWHM that will - be used for local sky subtraction; See :func:`objs_in_slit` and - :func:`~pypeit.core.skysub.local_skysub_extract`. - snr_thresh (:obj:`float`, optional): - SNR threshold for finding objects - specobj_dict (:obj:`dict`, optional): - Dictionary containing meta-data for the objects that will be - propagated into the :class:`~pypeit.specobj.SpecObj` objects. The - expected components are: - - - SLITID: The slit ID number - - DET: The detector identifier - - OBJTYPE: The object type - - PYPELINE: The class of pipeline algorithms applied - - If None, the dictionary is filled with the following placeholders:: - - specobj_dict = {'SLITID': 999, 'DET': 'DET01', - 'OBJTYPE': 'unknown', 'PYPELINE': 'unknown'} - - trim_edg (:obj:`tuple`, optional): - A two-tuple of integers or floats used to ignore objects within this - many pixels of the left and right slit boundaries, respectively. - show_peaks (:obj:`bool`, optional): - Plot the QA of the object peak finding in each order. - show_fits (:obj:`bool`, optional): - Plot trace fitting for final fits using PCA as crutch. - show_single_fits (:obj:`bool`, optional): - Plot trace fitting for single order fits. - show_trace (:obj:`bool`, optional): - Display the object traces on top of the image. - show_single_trace (:obj:`bool`, optional): - Display the object traces on top of the single order. - show_pca (:obj:`bool`, optional): - Display debugging plots for the PCA decomposition. - debug_all (:obj:`bool`, optional): - Show all the debugging plots. If True, this also overrides any - provided values for ``show_peaks``, ``show_trace``, and - ``show_pca``, setting them to True. - objfindQA_filename (:obj:`str`, optional): - Full path (directory and filename) for the object profile QA plot. - If None, not plot is produced and saved. - - Returns: - :class:`~pypeit.specobjs.SpecObjs`: Object containing the objects - detected. - """ - raise DeprecationWarning - msgs.error("This ginormous method as been Deprecated") - - #debug_all=True - if debug_all: - show_peaks = True - #show_fits = True - #show_single_fits = True - show_trace = True - show_pca = True - #show_single_trace = True - # TODO: This isn't used, right? - debug = True - - - if specobj_dict is None: - specobj_dict = {'SLITID': 999, 'ECH_ORDERINDX': 999, - 'DET': det, 'OBJTYPE': 'unknown', 'PYPELINE': 'Echelle'} - - # TODO Update FOF algorithm here with the one from scikit-learn. - - allmask = slitmask > -1 - if inmask is None: - inmask = allmask - - nspec, nspat = image.shape - norders = len(order_vec) - - # Find the spat IDs - gdslit_spat = np.unique(slitmask[slitmask >= 0]).astype(int) # Unique sorts - if gdslit_spat.size != np.sum(np.invert(maskslits)): - msgs.error('Masking of slitmask not in sync with that of maskslits. This is a bug') - #msgs.error('There is a mismatch between the number of valid orders found by PypeIt and ' - # 'the number expected for this spectrograph. Unable to continue. Please ' - # 'submit an issue on Github: https://github.com/pypeit/PypeIt/issues .') - - if spec_min_max is None: - spec_min_max = np.zeros((2,norders), dtype=int) - for iord in range(norders): - ispec, ispat = np.where(slitmask == gdslit_spat[iord]) - spec_min_max[:,iord] = ispec.min(), ispec.max() - - # Setup the plate scale - if isinstance(plate_scale,(float, int)): - plate_scale_ord = np.full(norders, plate_scale) - elif isinstance(plate_scale,(np.ndarray, list, tuple)): - if len(plate_scale) == norders: - plate_scale_ord = plate_scale - elif len(plate_scale) == 1: - plate_scale_ord = np.full(norders, plate_scale[0]) - else: - msgs.error('Invalid size for plate_scale. It must either have one element or norders elements') - else: - msgs.error('Invalid type for plate scale') - - specmid = nspec // 2 - spec_vec = np.arange(nspec) - slit_width = slit_righ - slit_left - slit_spec_pos = nspec/2.0 - - # TODO JFH This hand apertures in echelle needs to be completely refactored. - # Hand prep - # Determine the location of the source on *all* of the orders - if hand_extract_dict is not None: - f_spats = [] - for ss, spat, spec in zip(range(len(hand_extract_dict['spec'])), - hand_extract_dict['spat'], - hand_extract_dict['spec']): - # Find the input slit - ispec = int(np.clip(np.round(spec),0,nspec-1)) - ispat = int(np.clip(np.round(spat),0,nspat-1)) - slit = slitmask[ispec, ispat] - if slit == -1: - msgs.error('You are requesting a manual extraction at a position ' + - f'(spat, spec)={spat, spec} that is not on one of the echelle orders. Check your pypeit file.') - # Fractions - iord_hand = gdslit_spat.tolist().index(slit) - f_spat = (spat - slit_left[ispec, iord_hand]) / ( - slit_righ[ispec, iord_hand] - slit_left[ispec, iord_hand]) - f_spats.append(f_spat) - - # Loop over orders and find objects - sobjs = specobjs.SpecObjs() - # TODO: replace orderindx with the true order number here? Maybe not. Clean - # up SLITID and orderindx! - gdorders = np.arange(norders)[np.invert(maskslits)] - for iord in gdorders: #range(norders): - qa_title = 'Finding objects on order # {:d}'.format(order_vec[iord]) - msgs.info(qa_title) - thisslit_gpm = slitmask == gdslit_spat[iord] - inmask_iord = inmask & thisslit_gpm - specobj_dict['SLITID'] = gdslit_spat[iord] - specobj_dict['ECH_ORDERINDX'] = iord - specobj_dict['ECH_ORDER'] = order_vec[iord] - std_in = None if std_trace is None else std_trace[:, iord] - - # TODO JFH: Fix this. The way this code works, you should only need to create a single hand object, - # not one at every location on the order - if hand_extract_dict is not None: - new_hand_extract_dict = copy.deepcopy(hand_extract_dict) - for ss, spat, spec, f_spat in zip(range(len(hand_extract_dict['spec'])), - hand_extract_dict['spat'], - hand_extract_dict['spec'], f_spats): - ispec = int(spec) - new_hand_extract_dict['spec'][ss] = ispec - new_hand_extract_dict['spat'][ss] = slit_left[ispec,iord] + f_spat*( - slit_righ[ispec,iord]-slit_left[ispec,iord]) - else: - new_hand_extract_dict = None - - # Get SLTIORD_ID for the objfind QA - ech_objfindQA_filename = objfindQA_filename.replace('S0999', 'S{:04d}'.format(order_vec[iord])) \ - if objfindQA_filename is not None else None - # Run - sobjs_slit = \ - objs_in_slit(image, ivar, thisslit_gpm, slit_left[:,iord], slit_righ[:,iord], spec_min_max=spec_min_max[:,iord], - inmask=inmask_iord,std_trace=std_in, ncoeff=ncoeff, fwhm=fwhm, use_user_fwhm=use_user_fwhm, maxdev=maxdev, - hand_extract_dict=new_hand_extract_dict, nperslit=nperorder, extract_maskwidth=extract_maskwidth, - snr_thresh=snr_thresh, trim_edg=trim_edg, boxcar_rad=box_radius/plate_scale_ord[iord], - show_peaks=show_peaks, show_fits=show_single_fits, - show_trace=show_single_trace, qa_title=qa_title, specobj_dict=specobj_dict, - objfindQA_filename=ech_objfindQA_filename) - sobjs.add_sobj(sobjs_slit) - - nfound = len(sobjs) - - if nfound == 0: - msgs.warn('No objects found') - return sobjs - - FOF_frac = fof_link/(np.median(np.median(slit_width,axis=0)*plate_scale_ord)) - # Run the FOF. We use fake coordinates - fracpos = sobjs.SPAT_FRACPOS - ra_fake = fracpos/1000.0 # Divide all angles by 1000 to make geometry euclidian - dec_fake = np.zeros_like(fracpos) - if nfound>1: - inobj_id, multobj_id, firstobj_id, nextobj_id \ - = pydl.spheregroup(ra_fake, dec_fake, FOF_frac/1000.0) - # TODO spheregroup returns zero based indices but we use one based. We should probably add 1 to inobj_id here, - # i.e. obj_id_init = inobj_id + 1 - obj_id_init = inobj_id.copy() - elif nfound==1: - obj_id_init = np.zeros(1,dtype='int') - - uni_obj_id_init, uni_ind_init = np.unique(obj_id_init, return_index=True) - - # Now loop over the unique objects and check that there is only one object per order. If FOF - # grouped > 1 objects on the same order, then this will be popped out as its own unique object - obj_id = obj_id_init.copy() - nobj_init = len(uni_obj_id_init) - for iobj in range(nobj_init): - for iord in range(norders): - on_order = (obj_id_init == uni_obj_id_init[iobj]) & (sobjs.ECH_ORDERINDX == iord) - if (np.sum(on_order) > 1): - msgs.warn('Found multiple objects in a FOF group on order iord={:d}'.format(iord) + msgs.newline() + - 'Spawning new objects to maintain a single object per order.') - off_order = (obj_id_init == uni_obj_id_init[iobj]) & (sobjs.ECH_ORDERINDX != iord) - ind = np.where(on_order)[0] - if np.any(off_order): - # Keep the closest object to the location of the rest of the group (on other orders) - # as corresponding to this obj_id, and spawn new obj_ids for the others. - frac_mean = np.mean(fracpos[off_order]) - min_dist_ind = np.argmin(np.abs(fracpos[ind] - frac_mean)) - else: - # If there are no other objects with this obj_id to compare to, then we simply have multiple - # objects grouped together on the same order, so just spawn new object IDs for them to maintain - # one obj_id per order - min_dist_ind = 0 - ind_rest = np.setdiff1d(ind,ind[min_dist_ind]) - # JFH OLD LINE with bug - #obj_id[ind_rest] = (np.arange(len(ind_rest)) + 1) + obj_id_init.max() - obj_id[ind_rest] = (np.arange(len(ind_rest)) + 1) + obj_id.max() - - uni_obj_id, uni_ind = np.unique(obj_id, return_index=True) - nobj = len(uni_obj_id) - msgs.info('FOF matching found {:d}'.format(nobj) + ' unique objects') - - gfrac = np.zeros(nfound) - for jj in range(nobj): - this_obj_id = obj_id == uni_obj_id[jj] - gfrac[this_obj_id] = np.median(fracpos[this_obj_id]) - - uni_frac = gfrac[uni_ind] - - # Sort with respect to fractional slit location to guarantee that we have a similarly sorted list of objects later - isort_frac = uni_frac.argsort(kind='stable') - uni_obj_id = uni_obj_id[isort_frac] - uni_frac = uni_frac[isort_frac] - - sobjs_align = sobjs.copy() - # Loop over the orders and assign each specobj a fractional position and a obj_id number - for iobj in range(nobj): - for iord in range(norders): - on_order = (obj_id == uni_obj_id[iobj]) & (sobjs_align.ECH_ORDERINDX == iord) - sobjs_align[on_order].ECH_FRACPOS = uni_frac[iobj] - sobjs_align[on_order].ECH_OBJID = uni_obj_id[iobj] - sobjs_align[on_order].OBJID = uni_obj_id[iobj] - sobjs_align[on_order].ech_frac_was_fit = False - - # Reset names (just in case) - sobjs_align.set_names() - # Now loop over objects and fill in the missing objects and their traces. We will fit the fraction slit position of - # the good orders where an object was found and use that fit to predict the fractional slit position on the bad orders - # where no object was found - for iobj in range(nobj): - # Grab all the members of this obj_id from the object list - indx_obj_id = sobjs_align.ECH_OBJID == uni_obj_id[iobj] - nthisobj_id = np.sum(indx_obj_id) - # Perform the fit if this objects shows up on more than three orders - if (nthisobj_id > 3) and (nthisobj_id 1: - msgs.error('Problem in echelle object finding. The same objid={:d} appears {:d} times on echelle orderindx ={:d}' - ' even after duplicate obj_ids the orders were removed. ' - 'Report this bug to PypeIt developers'.format(uni_obj_id[iobj],num_on_order, iord)) - - - - # Loop over the objects and perform a quick and dirty extraction to assess S/N. - varimg = utils.calc_ivar(ivar) - flux_box = np.zeros((nspec, norders, nobj)) - ivar_box = np.zeros((nspec, norders, nobj)) - mask_box = np.zeros((nspec, norders, nobj)) - SNR_arr = np.zeros((norders, nobj)) - slitfracpos_arr = np.zeros((norders, nobj)) - for iobj in range(nobj): - for iord in range(norders): - iorder_vec = order_vec[iord] - indx = sobjs_align.slitorder_objid_indices(iorder_vec, uni_obj_id[iobj]) - #indx = (sobjs_align.ECH_OBJID == uni_obj_id[iobj]) & (sobjs_align.ECH_ORDERINDX == iord) - #spec = sobjs_align[indx][0] - inmask_iord = inmask & (slitmask == gdslit_spat[iord]) - # TODO make the snippet below its own function quick_extraction() - box_rad_pix = box_radius/plate_scale_ord[iord] - - # TODO -- We probably shouldn't be operating on a SpecObjs but instead a SpecObj - flux_tmp = moment1d(image*inmask_iord, sobjs_align[indx][0].TRACE_SPAT, 2*box_rad_pix, - row=sobjs_align[indx][0].trace_spec)[0] - var_tmp = moment1d(varimg*inmask_iord, sobjs_align[indx][0].TRACE_SPAT, 2*box_rad_pix, - row=sobjs_align[indx][0].trace_spec)[0] - ivar_tmp = utils.calc_ivar(var_tmp) - pixtot = moment1d(ivar*0 + 1.0, sobjs_align[indx][0].TRACE_SPAT, 2*box_rad_pix, - row=sobjs_align[indx][0].trace_spec)[0] - mask_tmp = moment1d(ivar*inmask_iord == 0.0, sobjs_align[indx][0].TRACE_SPAT, 2*box_rad_pix, - row=sobjs_align[indx][0].trace_spec)[0] != pixtot - - flux_box[:,iord,iobj] = flux_tmp*mask_tmp - ivar_box[:,iord,iobj] = np.fmax(ivar_tmp*mask_tmp,0.0) - mask_box[:,iord,iobj] = mask_tmp - mean, med_sn, stddev = astropy.stats.sigma_clipped_stats( - flux_box[mask_tmp,iord,iobj]*np.sqrt(ivar_box[mask_tmp,iord,iobj]), - sigma_lower=5.0,sigma_upper=5.0 - ) - # ToDO assign this to sobjs_align for use in the extraction - SNR_arr[iord,iobj] = med_sn - sobjs_align[indx][0].ech_snr = med_sn - # For hand extractions - slitfracpos_arr[iord,iobj] = sobjs_align[indx][0].SPAT_FRACPOS - - # Purge objects with low SNR that don't show up in enough orders, sort the list of objects with respect to obj_id - # and orderindx - keep_obj = np.zeros(nobj,dtype=bool) - sobjs_trim = specobjs.SpecObjs() - # objids are 1 based so that we can easily asign the negative to negative objects - iobj_keep = 1 - iobj_keep_not_hand = 1 - - # TODO JFH: Fix this ugly and dangerous hack that was added to accomodate hand apertures - hand_frac = [-1000] if hand_extract_dict is None else [int(np.round(ispat*1000)) for ispat in f_spats] - - ## Loop over objects from highest SNR to lowest SNR. Apply the S/N constraints. Once we hit the maximum number - # objects requested exit, except keep the hand apertures that were requested. - isort_SNR_max = np.argsort(np.median(SNR_arr,axis=0), kind='stable')[::-1] - for iobj in isort_SNR_max: - hand_ap_flag = int(np.round(slitfracpos_arr[0, iobj]*1000)) in hand_frac - SNR_constraint = (SNR_arr[:,iobj].max() > max_snr) or (np.sum(SNR_arr[:,iobj] > min_snr) >= nabove_min_snr) - nperorder_constraint = (iobj_keep-1) < nperorder - if (SNR_constraint and nperorder_constraint) or hand_ap_flag: - keep_obj[iobj] = True - ikeep = sobjs_align.ECH_OBJID == uni_obj_id[iobj] - sobjs_keep = sobjs_align[ikeep].copy() - sobjs_keep.ECH_OBJID = iobj_keep - sobjs_keep.OBJID = iobj_keep -# for spec in sobjs_keep: -# spec.ECH_OBJID = iobj_keep -# #spec.OBJID = iobj_keep - sobjs_trim.add_sobj(sobjs_keep[np.argsort(sobjs_keep.ECH_ORDERINDX, kind='stable')]) - iobj_keep += 1 - if not hand_ap_flag: - iobj_keep_not_hand += 1 - else: - if not nperorder_constraint: - msgs.info('Purging object #{:d}'.format(iobj) + - ' since there are already {:d} objects automatically identified ' - 'and you set nperorder={:d}'.format(iobj_keep_not_hand-1, nperorder)) - else: - msgs.info('Purging object #{:d}'.format(iobj) + ' which does not satisfy max_snr > {:5.2f} OR min_snr > {:5.2f}'.format(max_snr, min_snr) + - ' on at least nabove_min_snr >= {:d}'.format(nabove_min_snr) + ' orders') - - - nobj_trim = np.sum(keep_obj) - - if nobj_trim == 0: - msgs.warn('No objects found') - sobjs_final = specobjs.SpecObjs() - return sobjs_final - - # TODO JFH: We need to think about how to implement returning a maximum number of objects, where the objects - # returned are the highest S/N ones. It is a bit complicated with regards to the individual object finding and then - # the linking that is performed above, and also making sure the hand apertures don't get removed. - SNR_arr_trim = SNR_arr[:,keep_obj] - - - sobjs_final = sobjs_trim.copy() - # Loop over the objects one by one and adjust/predict the traces - pca_fits = np.zeros((nspec, norders, nobj_trim)) - - # Create the trc_inmask for iterative fitting below - trc_inmask = np.zeros((nspec, norders), dtype=bool) - for iord in range(norders): - trc_inmask[:,iord] = (spec_vec >= spec_min_max[0,iord]) & (spec_vec <= spec_min_max[1,iord]) - - for iobj in range(nobj_trim): - indx_obj_id = sobjs_final.ECH_OBJID == (iobj + 1) - # PCA predict all the orders now (where we have used the standard or slit boundary for the bad orders above) - msgs.info('Fitting echelle object finding PCA for object {:d}/{:d} with median SNR = {:5.3f}'.format( - iobj + 1,nobj_trim,np.median(sobjs_final[indx_obj_id].ech_snr))) - pca_fits[:,:,iobj] \ - = tracepca.pca_trace_object(sobjs_final[indx_obj_id].TRACE_SPAT.T, - order=coeff_npoly, npca=npca, - pca_explained_var=pca_explained_var, - trace_wgt=np.fmax(sobjs_final[indx_obj_id].ech_snr, 1.0)**2, - debug=show_pca) - - # Trial and error shows weighting by S/N instead of S/N^2 performs better - # JXP -- Updated to now be S/N**2, i.e. inverse variance, with fitting fit - - # Perform iterative flux weighted centroiding using new PCA predictions - xinit_fweight = pca_fits[:,:,iobj].copy() - inmask_now = inmask & allmask - xfit_fweight = fit_trace(image, xinit_fweight, ncoeff, bpm=np.invert(inmask_now), - trace_bpm=np.invert(trc_inmask), fwhm=fwhm, maxdev=maxdev, - debug=show_fits)[0] - - # Perform iterative Gaussian weighted centroiding - xinit_gweight = xfit_fweight.copy() - xfit_gweight = fit_trace(image, xinit_gweight, ncoeff, bpm=np.invert(inmask_now), - trace_bpm=np.invert(trc_inmask), weighting='gaussian', fwhm=fwhm, - maxdev=maxdev, debug=show_fits)[0] - - #TODO Assign the new traces. Only assign the orders that were not orginally detected and traced. If this works - # well, we will avoid doing all of the iter_tracefits above to make the code faster. - for iord, spec in enumerate(sobjs_final[indx_obj_id]): - # JFH added the condition on ech_frac_was_fit with S/N cut on 7-7-19. - # TODO is this robust against half the order being masked? - if spec.ech_frac_was_fit & (spec.ech_snr > 1.0): - spec.TRACE_SPAT = xfit_gweight[:,iord] - spec.SPAT_PIXPOS = spec.TRACE_SPAT[specmid] - - #TODO Put in some criterion here that does not let the fractional position change too much during the iterative - # tracefitting. The problem is spurious apertures identified on one slit can be pulled over to the center of flux - # resulting in a bunch of objects landing on top of each other. - - # Set the IDs - sobjs_final[:].ECH_ORDER = order_vec[sobjs_final[:].ECH_ORDERINDX] - #for spec in sobjs_final: - # spec.ech_order = order_vec[spec.ECH_ORDERINDX] - sobjs_final.set_names() - - if show_trace: - viewer, ch = display.show_image(image*allmask) - - for spec in sobjs_trim: - color = 'red' if spec.ech_frac_was_fit else 'magenta' - ## Showing the final flux weighted centroiding from PCA predictions - display.show_trace(viewer, ch, spec.TRACE_SPAT, spec.NAME, color=color) - - for iobj in range(nobj_trim): - for iord in range(norders): - ## Showing PCA predicted locations before recomputing flux/gaussian weighted centroiding - display.show_trace(viewer, ch, pca_fits[:,iord, iobj], str(uni_frac[iobj]), color='yellow') - ## Showing the final traces from this routine - display.show_trace(viewer, ch, sobjs_final.TRACE_SPAT[iord].T, sobjs_final.NAME, color='cyan') - - # Labels for the points - text_final = [dict(type='text', args=(nspat / 2 -40, nspec / 2, 'final trace'), - kwargs=dict(color='cyan', fontsize=20))] - - text_pca = [dict(type='text', args=(nspat / 2 -40, nspec / 2 - 30, 'PCA fit'),kwargs=dict(color='yellow', fontsize=20))] - - text_fit = [dict(type='text', args=(nspat / 2 -40, nspec / 2 - 60, 'predicted'),kwargs=dict(color='red', fontsize=20))] - - text_notfit = [dict(type='text', args=(nspat / 2 -40, nspec / 2 - 90, 'originally found'),kwargs=dict(color='magenta', fontsize=20))] - - canvas = viewer.canvas(ch._chname) - canvas_list = text_final + text_pca + text_fit + text_notfit - canvas.add('constructedcanvas', canvas_list) - # TODO two things need to be debugged. 1) For objects which were found and traced, i don't think we should be updating the tracing with - # the PCA. This just adds a failutre mode. 2) The PCA fit is going wild for X-shooter. Debug that. - # Vette - for sobj in sobjs_final: - if not sobj.ready_for_extraction(): - msgs.error("Bad SpecObj. Can't proceed") - - return sobjs_final - - def objfind_QA(spat_peaks, snr_peaks, spat_vector, snr_vector, snr_thresh, qa_title, peak_gpm, near_edge_bpm, nperslit_bpm, objfindQA_filename=None, show=False): From 0ac9fe3578ae329eb13711c0183166f8297a470a Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Tue, 13 Aug 2024 14:30:27 -1000 Subject: [PATCH 15/17] move orig_ech_objfind to deprecated2 --- deprecated/old_ech_objfind.py | 681 ++++++++++++++++++++++++++++++++++ 1 file changed, 681 insertions(+) create mode 100644 deprecated/old_ech_objfind.py diff --git a/deprecated/old_ech_objfind.py b/deprecated/old_ech_objfind.py new file mode 100644 index 0000000000..e0e7be9309 --- /dev/null +++ b/deprecated/old_ech_objfind.py @@ -0,0 +1,681 @@ + +def orig_ech_objfind(image, ivar, slitmask, slit_left, slit_righ, order_vec, maskslits, det='DET01', + inmask=None, spec_min_max=None, fof_link=1.5, plate_scale=0.2, + std_trace=None, ncoeff=5, npca=None, coeff_npoly=None, max_snr=2.0, min_snr=1.0, + nabove_min_snr=2, pca_explained_var=99.0, box_radius=2.0, fwhm=3.0, + use_user_fwhm=False, maxdev=2.0, hand_extract_dict=None, nperorder=2, + extract_maskwidth=3.0, snr_thresh=10.0, + specobj_dict=None, trim_edg=(5, 5), + show_peaks=False, show_fits=False, show_single_fits=False, + show_trace=False, show_single_trace=False, show_pca=False, + debug_all=False, objfindQA_filename=None): + """ + Object finding routine for Echelle spectrographs. + + This routine: + + #. Runs object finding on each order individually + + #. Links the objects found together using a friends-of-friends algorithm + on fractional order position. + + #. For objects which were only found on some orders, the standard (or + the slit boundaries) are placed at the appropriate fractional + position along the order. + + #. A PCA fit to the traces is performed using the routine above pca_fit + + Args: + image (`numpy.ndarray`_): + (Floating-point) Image to use for object search with shape (nspec, + nspat). The first dimension (nspec) is spectral, and second + dimension (nspat) is spatial. Note this image can either have the + sky background in it, or have already been sky subtracted. Object + finding works best on sky-subtracted images. Ideally, object finding + is run in another routine, global sky-subtraction performed, and + then this code should be run. However, it is also possible to run + this code on non-sky-subtracted images. + ivar (`numpy.ndarray`_): + Floating-point inverse variance image for the input image. Shape + must match ``image``, (nspec, nspat). + slitmask (`numpy.ndarray`_): + Integer image indicating the pixels that belong to each order. + Pixels that are not on an order have value -1, and those that are on + an order have a value equal to the slit number (i.e. 0 to nslits-1 + from left to right on the image). Shape must match ``image``, + (nspec, nspat). + slit_left (`numpy.ndarray`_): + Left boundary of orders to be extracted (given as floating point + pixels). Shape is (nspec, norders), where norders is the total + number of traced echelle orders. + slit_righ (`numpy.ndarray`_): + Right boundary of orders to be extracted (given as floating point + pixels). Shape is (nspec, norders), where norders is the total + number of traced echelle orders. + order_vec (`numpy.ndarray`_): + Vector identifying the Echelle orders for each pair of order edges + found. This is saved to the output :class:`~pypeit.specobj.SpecObj` + objects. If the orders are not known, this can be + ``np.arange(norders)`` (but this is *not* recommended). + maskslits (`numpy.ndarray`_): + Boolean array selecting orders that should be ignored (i.e., good + orders are False, bad orders are True). Shape must be (norders,). + det (:obj:`str`, optional): + The name of the detector containing the object. Only used if + ``specobj_dict`` is None. + inmask (`numpy.ndarray`_, optional): + Good-pixel mask for input image. Must have the same shape as + ``image``. If None, all pixels in ``slitmask`` with non-negative + values are considered good. + spec_min_max (`numpy.ndarray`_, optional): + 2D array defining the minimum and maximum pixel in the spectral + direction with useable data for each order. Shape must be (2, + norders). This should only be used for echelle spectrographs for + which the orders do not entirely cover the detector. PCA tracing + will re-map the traces such that they all have the same length, + compute the PCA, and then re-map the orders back. This improves + performance for echelle spectrographs by removing the nonlinear + shrinking of the orders so that the linear pca operation can better + predict the traces. If None, the minimum and maximum values will be + determined automatically from ``slitmask``. + fof_link (:obj:`float`, optional): + Friends-of-friends linking length in arcseconds used to link + together traces across orders. The routine links together at + the same fractional slit position and links them together + with a friends-of-friends algorithm using this linking + length. + plate_scale (:obj:`float`, `numpy.ndarray`_, optional): + Plate scale in arcsec/pix. This can either be a single float for + every order, or an array with shape (norders,) providing the plate + scale of each order. + std_trace (`numpy.ndarray`_, optional): + Vector with the standard star trace, which is used as a crutch for + tracing. Shape must be (nspec,). If None, the slit boundaries are + used as the crutch. + ncoeff (:obj:`int`, optional): + Order of polynomial fit to traces. + npca (:obj:`int`, optional): + Number of PCA components to keep during PCA decomposition of the + object traces. If None, the number of components set by requiring + the PCA accounts for approximately 99% of the variance. + coeff_npoly (:obj:`int`, optional): + Order of polynomial used for PCA coefficients fitting. If None, + value set automatically, see + :func:`~pypeit.tracepca.pca_trace_object`. + max_snr (:obj:`float`, optional): + For an object to be included in the output object, it must have a + max S/N ratio above this value. + min_snr (:obj:`float`, optional): + For an object to be included in the output object, it must have a + a median S/N ratio above this value for at least + ``nabove_min_snr`` orders (see below). + nabove_min_snr (:obj:`int`, optional): + The required number of orders that an object must have with median + SNR greater than ``min_snr`` in order to be included in the output + object. + pca_explained_var (:obj:`float`, optional): + The percentage (i.e., not the fraction) of the variance in the data + accounted for by the PCA used to truncate the number of PCA + coefficients to keep (see ``npca``). Ignored if ``npca`` is provided + directly; see :func:`~pypeit.tracepca.pca_trace_object`. + box_radius (:obj:`float`, optional): + Box_car extraction radius in arcseconds to assign to each detected + object and to be used later for boxcar extraction. In this method + ``box_radius`` is converted into pixels using ``plate_scale``. + ``box_radius`` is also used for SNR calculation and trimming. + fwhm (:obj:`float`, optional): + Estimated fwhm of the objects in pixels + use_user_fwhm (:obj:`bool`, optional): + If True, ``PypeIt`` will use the spatial profile FWHM input by the + user (see ``fwhm``) rather than determine the spatial FWHM from the + smashed spatial profile via the automated algorithm. + maxdev (:obj:`float`, optional): + Maximum deviation of pixels from polynomial fit to trace + used to reject bad pixels in trace fitting. + hand_extract_dict (:obj:`dict`, optional): + Dictionary with info on manual extraction; see + :class:`~pypeit.manual_extract.ManualExtractionObj`. + nperorder (:obj:`int`, optional): + Maximum number of objects allowed per order. If there are more + detections than this number, the code will select the ``nperorder`` + most significant detections. However, hand apertures will always be + returned and do not count against this budget. + extract_maskwidth (:obj:`float`, optional): + Determines the initial size of the region in units of FWHM that will + be used for local sky subtraction; See :func:`objs_in_slit` and + :func:`~pypeit.core.skysub.local_skysub_extract`. + snr_thresh (:obj:`float`, optional): + SNR threshold for finding objects + specobj_dict (:obj:`dict`, optional): + Dictionary containing meta-data for the objects that will be + propagated into the :class:`~pypeit.specobj.SpecObj` objects. The + expected components are: + + - SLITID: The slit ID number + - DET: The detector identifier + - OBJTYPE: The object type + - PYPELINE: The class of pipeline algorithms applied + + If None, the dictionary is filled with the following placeholders:: + + specobj_dict = {'SLITID': 999, 'DET': 'DET01', + 'OBJTYPE': 'unknown', 'PYPELINE': 'unknown'} + + trim_edg (:obj:`tuple`, optional): + A two-tuple of integers or floats used to ignore objects within this + many pixels of the left and right slit boundaries, respectively. + show_peaks (:obj:`bool`, optional): + Plot the QA of the object peak finding in each order. + show_fits (:obj:`bool`, optional): + Plot trace fitting for final fits using PCA as crutch. + show_single_fits (:obj:`bool`, optional): + Plot trace fitting for single order fits. + show_trace (:obj:`bool`, optional): + Display the object traces on top of the image. + show_single_trace (:obj:`bool`, optional): + Display the object traces on top of the single order. + show_pca (:obj:`bool`, optional): + Display debugging plots for the PCA decomposition. + debug_all (:obj:`bool`, optional): + Show all the debugging plots. If True, this also overrides any + provided values for ``show_peaks``, ``show_trace``, and + ``show_pca``, setting them to True. + objfindQA_filename (:obj:`str`, optional): + Full path (directory and filename) for the object profile QA plot. + If None, not plot is produced and saved. + + Returns: + :class:`~pypeit.specobjs.SpecObjs`: Object containing the objects + detected. + """ + raise DeprecationWarning + msgs.error("This ginormous method as been Deprecated") + + # debug_all=True + if debug_all: + show_peaks = True + # show_fits = True + # show_single_fits = True + show_trace = True + show_pca = True + # show_single_trace = True + # TODO: This isn't used, right? + debug = True + + if specobj_dict is None: + specobj_dict = {'SLITID': 999, 'ECH_ORDERINDX': 999, + 'DET': det, 'OBJTYPE': 'unknown', 'PYPELINE': 'Echelle'} + + # TODO Update FOF algorithm here with the one from scikit-learn. + + allmask = slitmask > -1 + if inmask is None: + inmask = allmask + + nspec, nspat = image.shape + norders = len(order_vec) + + # Find the spat IDs + gdslit_spat = np.unique(slitmask[slitmask >= 0]).astype(int) # Unique sorts + if gdslit_spat.size != np.sum(np.invert(maskslits)): + msgs.error('Masking of slitmask not in sync with that of maskslits. This is a bug') + # msgs.error('There is a mismatch between the number of valid orders found by PypeIt and ' + # 'the number expected for this spectrograph. Unable to continue. Please ' + # 'submit an issue on Github: https://github.com/pypeit/PypeIt/issues .') + + if spec_min_max is None: + spec_min_max = np.zeros((2, norders), dtype=int) + for iord in range(norders): + ispec, ispat = np.where(slitmask == gdslit_spat[iord]) + spec_min_max[:, iord] = ispec.min(), ispec.max() + + # Setup the plate scale + if isinstance(plate_scale, (float, int)): + plate_scale_ord = np.full(norders, plate_scale) + elif isinstance(plate_scale, (np.ndarray, list, tuple)): + if len(plate_scale) == norders: + plate_scale_ord = plate_scale + elif len(plate_scale) == 1: + plate_scale_ord = np.full(norders, plate_scale[0]) + else: + msgs.error('Invalid size for plate_scale. It must either have one element or norders elements') + else: + msgs.error('Invalid type for plate scale') + + specmid = nspec // 2 + spec_vec = np.arange(nspec) + slit_width = slit_righ - slit_left + slit_spec_pos = nspec / 2.0 + + # TODO JFH This hand apertures in echelle needs to be completely refactored. + # Hand prep + # Determine the location of the source on *all* of the orders + if hand_extract_dict is not None: + f_spats = [] + for ss, spat, spec in zip(range(len(hand_extract_dict['spec'])), + hand_extract_dict['spat'], + hand_extract_dict['spec']): + # Find the input slit + ispec = int(np.clip(np.round(spec), 0, nspec - 1)) + ispat = int(np.clip(np.round(spat), 0, nspat - 1)) + slit = slitmask[ispec, ispat] + if slit == -1: + msgs.error('You are requesting a manual extraction at a position ' + + f'(spat, spec)={spat, spec} that is not on one of the echelle orders. Check your pypeit file.') + # Fractions + iord_hand = gdslit_spat.tolist().index(slit) + f_spat = (spat - slit_left[ispec, iord_hand]) / ( + slit_righ[ispec, iord_hand] - slit_left[ispec, iord_hand]) + f_spats.append(f_spat) + + # Loop over orders and find objects + sobjs = specobjs.SpecObjs() + # TODO: replace orderindx with the true order number here? Maybe not. Clean + # up SLITID and orderindx! + gdorders = np.arange(norders)[np.invert(maskslits)] + for iord in gdorders: # range(norders): + qa_title = 'Finding objects on order # {:d}'.format(order_vec[iord]) + msgs.info(qa_title) + thisslit_gpm = slitmask == gdslit_spat[iord] + inmask_iord = inmask & thisslit_gpm + specobj_dict['SLITID'] = gdslit_spat[iord] + specobj_dict['ECH_ORDERINDX'] = iord + specobj_dict['ECH_ORDER'] = order_vec[iord] + std_in = None if std_trace is None else std_trace[:, iord] + + # TODO JFH: Fix this. The way this code works, you should only need to create a single hand object, + # not one at every location on the order + if hand_extract_dict is not None: + new_hand_extract_dict = copy.deepcopy(hand_extract_dict) + for ss, spat, spec, f_spat in zip(range(len(hand_extract_dict['spec'])), + hand_extract_dict['spat'], + hand_extract_dict['spec'], f_spats): + ispec = int(spec) + new_hand_extract_dict['spec'][ss] = ispec + new_hand_extract_dict['spat'][ss] = slit_left[ispec, iord] + f_spat * ( + slit_righ[ispec, iord] - slit_left[ispec, iord]) + else: + new_hand_extract_dict = None + + # Get SLTIORD_ID for the objfind QA + ech_objfindQA_filename = objfindQA_filename.replace('S0999', 'S{:04d}'.format(order_vec[iord])) \ + if objfindQA_filename is not None else None + # Run + sobjs_slit = \ + objs_in_slit(image, ivar, thisslit_gpm, slit_left[:, iord], slit_righ[:, iord], + spec_min_max=spec_min_max[:, iord], + inmask=inmask_iord, std_trace=std_in, ncoeff=ncoeff, fwhm=fwhm, use_user_fwhm=use_user_fwhm, + maxdev=maxdev, + hand_extract_dict=new_hand_extract_dict, nperslit=nperorder, + extract_maskwidth=extract_maskwidth, + snr_thresh=snr_thresh, trim_edg=trim_edg, boxcar_rad=box_radius / plate_scale_ord[iord], + show_peaks=show_peaks, show_fits=show_single_fits, + show_trace=show_single_trace, qa_title=qa_title, specobj_dict=specobj_dict, + objfindQA_filename=ech_objfindQA_filename) + sobjs.add_sobj(sobjs_slit) + + nfound = len(sobjs) + + if nfound == 0: + msgs.warn('No objects found') + return sobjs + + FOF_frac = fof_link / (np.median(np.median(slit_width, axis=0) * plate_scale_ord)) + # Run the FOF. We use fake coordinates + fracpos = sobjs.SPAT_FRACPOS + ra_fake = fracpos / 1000.0 # Divide all angles by 1000 to make geometry euclidian + dec_fake = np.zeros_like(fracpos) + if nfound > 1: + inobj_id, multobj_id, firstobj_id, nextobj_id \ + = pydl.spheregroup(ra_fake, dec_fake, FOF_frac / 1000.0) + # TODO spheregroup returns zero based indices but we use one based. We should probably add 1 to inobj_id here, + # i.e. obj_id_init = inobj_id + 1 + obj_id_init = inobj_id.copy() + elif nfound == 1: + obj_id_init = np.zeros(1, dtype='int') + + uni_obj_id_init, uni_ind_init = np.unique(obj_id_init, return_index=True) + + # Now loop over the unique objects and check that there is only one object per order. If FOF + # grouped > 1 objects on the same order, then this will be popped out as its own unique object + obj_id = obj_id_init.copy() + nobj_init = len(uni_obj_id_init) + for iobj in range(nobj_init): + for iord in range(norders): + on_order = (obj_id_init == uni_obj_id_init[iobj]) & (sobjs.ECH_ORDERINDX == iord) + if (np.sum(on_order) > 1): + msgs.warn('Found multiple objects in a FOF group on order iord={:d}'.format(iord) + msgs.newline() + + 'Spawning new objects to maintain a single object per order.') + off_order = (obj_id_init == uni_obj_id_init[iobj]) & (sobjs.ECH_ORDERINDX != iord) + ind = np.where(on_order)[0] + if np.any(off_order): + # Keep the closest object to the location of the rest of the group (on other orders) + # as corresponding to this obj_id, and spawn new obj_ids for the others. + frac_mean = np.mean(fracpos[off_order]) + min_dist_ind = np.argmin(np.abs(fracpos[ind] - frac_mean)) + else: + # If there are no other objects with this obj_id to compare to, then we simply have multiple + # objects grouped together on the same order, so just spawn new object IDs for them to maintain + # one obj_id per order + min_dist_ind = 0 + ind_rest = np.setdiff1d(ind, ind[min_dist_ind]) + # JFH OLD LINE with bug + # obj_id[ind_rest] = (np.arange(len(ind_rest)) + 1) + obj_id_init.max() + obj_id[ind_rest] = (np.arange(len(ind_rest)) + 1) + obj_id.max() + + uni_obj_id, uni_ind = np.unique(obj_id, return_index=True) + nobj = len(uni_obj_id) + msgs.info('FOF matching found {:d}'.format(nobj) + ' unique objects') + + gfrac = np.zeros(nfound) + for jj in range(nobj): + this_obj_id = obj_id == uni_obj_id[jj] + gfrac[this_obj_id] = np.median(fracpos[this_obj_id]) + + uni_frac = gfrac[uni_ind] + + # Sort with respect to fractional slit location to guarantee that we have a similarly sorted list of objects later + isort_frac = uni_frac.argsort(kind='stable') + uni_obj_id = uni_obj_id[isort_frac] + uni_frac = uni_frac[isort_frac] + + sobjs_align = sobjs.copy() + # Loop over the orders and assign each specobj a fractional position and a obj_id number + for iobj in range(nobj): + for iord in range(norders): + on_order = (obj_id == uni_obj_id[iobj]) & (sobjs_align.ECH_ORDERINDX == iord) + sobjs_align[on_order].ECH_FRACPOS = uni_frac[iobj] + sobjs_align[on_order].ECH_OBJID = uni_obj_id[iobj] + sobjs_align[on_order].OBJID = uni_obj_id[iobj] + sobjs_align[on_order].ech_frac_was_fit = False + + # Reset names (just in case) + sobjs_align.set_names() + # Now loop over objects and fill in the missing objects and their traces. We will fit the fraction slit position of + # the good orders where an object was found and use that fit to predict the fractional slit position on the bad orders + # where no object was found + for iobj in range(nobj): + # Grab all the members of this obj_id from the object list + indx_obj_id = sobjs_align.ECH_OBJID == uni_obj_id[iobj] + nthisobj_id = np.sum(indx_obj_id) + # Perform the fit if this objects shows up on more than three orders + if (nthisobj_id > 3) and (nthisobj_id < norders): + thisorderindx = sobjs_align[indx_obj_id].ECH_ORDERINDX + goodorder = np.zeros(norders, dtype=bool) + goodorder[thisorderindx] = True + badorder = np.invert(goodorder) + xcen_good = (sobjs_align[indx_obj_id].TRACE_SPAT).T + slit_frac_good = (xcen_good - slit_left[:, goodorder]) / slit_width[:, goodorder] + # Fractional slit position averaged across the spectral direction for each order + frac_mean_good = np.mean(slit_frac_good, 0) + # Perform a linear fit to fractional slit position + # TODO Do this as a S/N weighted fit similar to what is now in the pca_trace algorithm? + # msk_frac, poly_coeff_frac = fitting.robust_fit(order_vec[goodorder], frac_mean_good, 1, + pypeitFit = fitting.robust_fit(order_vec[goodorder], frac_mean_good, 1, + function='polynomial', maxiter=20, lower=2, upper=2, + use_mad=True, sticky=False, + minx=order_vec.min(), maxx=order_vec.max()) + frac_mean_new = np.zeros(norders) + frac_mean_new[badorder] = pypeitFit.eval( + order_vec[badorder]) # , minx = order_vec.min(),maxx=order_vec.max()) + frac_mean_new[goodorder] = frac_mean_good + # TODO This QA needs some work + if show_pca: + frac_mean_fit = pypeitFit.eval(order_vec) + plt.plot(order_vec[goodorder][pypeitFit.bool_gpm], frac_mean_new[goodorder][pypeitFit.bool_gpm], 'ko', + mfc='k', markersize=8.0, label='Good Orders Kept') + plt.plot(order_vec[goodorder][np.invert(pypeitFit.bool_gpm)], + frac_mean_new[goodorder][np.invert(pypeitFit.bool_gpm)], 'ro', mfc='k', markersize=8.0, + label='Good Orders Rejected') + plt.plot(order_vec[badorder], frac_mean_new[badorder], 'ko', mfc='None', markersize=8.0, + label='Predicted Bad Orders') + plt.plot(order_vec, frac_mean_new, '+', color='cyan', markersize=12.0, label='Final Order Fraction') + plt.plot(order_vec, frac_mean_fit, 'r-', label='Fractional Order Position Fit') + plt.xlabel('Order Index', fontsize=14) + plt.ylabel('Fractional Slit Position', fontsize=14) + plt.title('Fractional Slit Position Fit') + plt.legend() + plt.show() + else: + frac_mean_new = np.full(norders, uni_frac[iobj]) + + # Now loop over the orders and add objects on the ordrers for which the current object was not found + for iord in range(norders): + # Is the current object detected on this order? + on_order = (sobjs_align.ECH_OBJID == uni_obj_id[iobj]) & (sobjs_align.ECH_ORDERINDX == iord) + num_on_order = np.sum(on_order) + if num_on_order == 0: + # If it is not, create a new sobjs and add to sobjs_align and assign required tags + thisobj = specobj.SpecObj('Echelle', sobjs_align[0].DET, + OBJTYPE=sobjs_align[0].OBJTYPE, + ECH_ORDERINDX=iord, + ECH_ORDER=order_vec[iord]) + # thisobj.ECH_ORDERINDX = iord + # thisobj.ech_order = order_vec[iord] + thisobj.SPAT_FRACPOS = uni_frac[iobj] + # Assign traces using the fractional position fit above + if std_trace is not None: + x_trace = np.interp(slit_spec_pos, spec_vec, std_trace[:, iord]) + shift = np.interp(slit_spec_pos, spec_vec, + slit_left[:, iord] + slit_width[:, iord] * frac_mean_new[iord]) - x_trace + thisobj.TRACE_SPAT = std_trace[:, iord] + shift + else: + thisobj.TRACE_SPAT = slit_left[:, iord] + slit_width[:, iord] * frac_mean_new[iord] # new trace + thisobj.trace_spec = spec_vec + thisobj.SPAT_PIXPOS = thisobj.TRACE_SPAT[specmid] + # Use the real detections of this objects for the FWHM + this_obj_id = obj_id == uni_obj_id[iobj] + # Assign to the fwhm of the nearest detected order + imin = np.argmin(np.abs(sobjs_align[this_obj_id].ECH_ORDERINDX - iord)) + thisobj.FWHM = sobjs_align[imin].FWHM + thisobj.maskwidth = sobjs_align[imin].maskwidth + thisobj.smash_peakflux = sobjs_align[imin].smash_peakflux + thisobj.smash_snr = sobjs_align[imin].smash_snr + thisobj.BOX_RADIUS = sobjs_align[imin].BOX_RADIUS + thisobj.ECH_FRACPOS = uni_frac[iobj] + thisobj.ECH_OBJID = uni_obj_id[iobj] + thisobj.OBJID = uni_obj_id[iobj] + thisobj.SLITID = gdslit_spat[iord] + thisobj.ech_frac_was_fit = True + thisobj.set_name() + sobjs_align.add_sobj(thisobj) + obj_id = np.append(obj_id, uni_obj_id[iobj]) + gfrac = np.append(gfrac, uni_frac[iobj]) + elif num_on_order == 1: + # Object is already on this order so no need to do anything + pass + elif num_on_order > 1: + msgs.error( + 'Problem in echelle object finding. The same objid={:d} appears {:d} times on echelle orderindx ={:d}' + ' even after duplicate obj_ids the orders were removed. ' + 'Report this bug to PypeIt developers'.format(uni_obj_id[iobj], num_on_order, iord)) + + # Loop over the objects and perform a quick and dirty extraction to assess S/N. + varimg = utils.calc_ivar(ivar) + flux_box = np.zeros((nspec, norders, nobj)) + ivar_box = np.zeros((nspec, norders, nobj)) + mask_box = np.zeros((nspec, norders, nobj)) + SNR_arr = np.zeros((norders, nobj)) + slitfracpos_arr = np.zeros((norders, nobj)) + for iobj in range(nobj): + for iord in range(norders): + iorder_vec = order_vec[iord] + indx = sobjs_align.slitorder_objid_indices(iorder_vec, uni_obj_id[iobj]) + # indx = (sobjs_align.ECH_OBJID == uni_obj_id[iobj]) & (sobjs_align.ECH_ORDERINDX == iord) + # spec = sobjs_align[indx][0] + inmask_iord = inmask & (slitmask == gdslit_spat[iord]) + # TODO make the snippet below its own function quick_extraction() + box_rad_pix = box_radius / plate_scale_ord[iord] + + # TODO -- We probably shouldn't be operating on a SpecObjs but instead a SpecObj + flux_tmp = moment1d(image * inmask_iord, sobjs_align[indx][0].TRACE_SPAT, 2 * box_rad_pix, + row=sobjs_align[indx][0].trace_spec)[0] + var_tmp = moment1d(varimg * inmask_iord, sobjs_align[indx][0].TRACE_SPAT, 2 * box_rad_pix, + row=sobjs_align[indx][0].trace_spec)[0] + ivar_tmp = utils.calc_ivar(var_tmp) + pixtot = moment1d(ivar * 0 + 1.0, sobjs_align[indx][0].TRACE_SPAT, 2 * box_rad_pix, + row=sobjs_align[indx][0].trace_spec)[0] + mask_tmp = moment1d(ivar * inmask_iord == 0.0, sobjs_align[indx][0].TRACE_SPAT, 2 * box_rad_pix, + row=sobjs_align[indx][0].trace_spec)[0] != pixtot + + flux_box[:, iord, iobj] = flux_tmp * mask_tmp + ivar_box[:, iord, iobj] = np.fmax(ivar_tmp * mask_tmp, 0.0) + mask_box[:, iord, iobj] = mask_tmp + mean, med_sn, stddev = astropy.stats.sigma_clipped_stats( + flux_box[mask_tmp, iord, iobj] * np.sqrt(ivar_box[mask_tmp, iord, iobj]), + sigma_lower=5.0, sigma_upper=5.0 + ) + # ToDO assign this to sobjs_align for use in the extraction + SNR_arr[iord, iobj] = med_sn + sobjs_align[indx][0].ech_snr = med_sn + # For hand extractions + slitfracpos_arr[iord, iobj] = sobjs_align[indx][0].SPAT_FRACPOS + + # Purge objects with low SNR that don't show up in enough orders, sort the list of objects with respect to obj_id + # and orderindx + keep_obj = np.zeros(nobj, dtype=bool) + sobjs_trim = specobjs.SpecObjs() + # objids are 1 based so that we can easily asign the negative to negative objects + iobj_keep = 1 + iobj_keep_not_hand = 1 + + # TODO JFH: Fix this ugly and dangerous hack that was added to accomodate hand apertures + hand_frac = [-1000] if hand_extract_dict is None else [int(np.round(ispat * 1000)) for ispat in f_spats] + + ## Loop over objects from highest SNR to lowest SNR. Apply the S/N constraints. Once we hit the maximum number + # objects requested exit, except keep the hand apertures that were requested. + isort_SNR_max = np.argsort(np.median(SNR_arr, axis=0), kind='stable')[::-1] + for iobj in isort_SNR_max: + hand_ap_flag = int(np.round(slitfracpos_arr[0, iobj] * 1000)) in hand_frac + SNR_constraint = (SNR_arr[:, iobj].max() > max_snr) or (np.sum(SNR_arr[:, iobj] > min_snr) >= nabove_min_snr) + nperorder_constraint = (iobj_keep - 1) < nperorder + if (SNR_constraint and nperorder_constraint) or hand_ap_flag: + keep_obj[iobj] = True + ikeep = sobjs_align.ECH_OBJID == uni_obj_id[iobj] + sobjs_keep = sobjs_align[ikeep].copy() + sobjs_keep.ECH_OBJID = iobj_keep + sobjs_keep.OBJID = iobj_keep + # for spec in sobjs_keep: + # spec.ECH_OBJID = iobj_keep + # #spec.OBJID = iobj_keep + sobjs_trim.add_sobj(sobjs_keep[np.argsort(sobjs_keep.ECH_ORDERINDX, kind='stable')]) + iobj_keep += 1 + if not hand_ap_flag: + iobj_keep_not_hand += 1 + else: + if not nperorder_constraint: + msgs.info('Purging object #{:d}'.format(iobj) + + ' since there are already {:d} objects automatically identified ' + 'and you set nperorder={:d}'.format(iobj_keep_not_hand - 1, nperorder)) + else: + msgs.info('Purging object #{:d}'.format( + iobj) + ' which does not satisfy max_snr > {:5.2f} OR min_snr > {:5.2f}'.format(max_snr, min_snr) + + ' on at least nabove_min_snr >= {:d}'.format(nabove_min_snr) + ' orders') + + nobj_trim = np.sum(keep_obj) + + if nobj_trim == 0: + msgs.warn('No objects found') + sobjs_final = specobjs.SpecObjs() + return sobjs_final + + # TODO JFH: We need to think about how to implement returning a maximum number of objects, where the objects + # returned are the highest S/N ones. It is a bit complicated with regards to the individual object finding and then + # the linking that is performed above, and also making sure the hand apertures don't get removed. + SNR_arr_trim = SNR_arr[:, keep_obj] + + sobjs_final = sobjs_trim.copy() + # Loop over the objects one by one and adjust/predict the traces + pca_fits = np.zeros((nspec, norders, nobj_trim)) + + # Create the trc_inmask for iterative fitting below + trc_inmask = np.zeros((nspec, norders), dtype=bool) + for iord in range(norders): + trc_inmask[:, iord] = (spec_vec >= spec_min_max[0, iord]) & (spec_vec <= spec_min_max[1, iord]) + + for iobj in range(nobj_trim): + indx_obj_id = sobjs_final.ECH_OBJID == (iobj + 1) + # PCA predict all the orders now (where we have used the standard or slit boundary for the bad orders above) + msgs.info('Fitting echelle object finding PCA for object {:d}/{:d} with median SNR = {:5.3f}'.format( + iobj + 1, nobj_trim, np.median(sobjs_final[indx_obj_id].ech_snr))) + pca_fits[:, :, iobj] \ + = tracepca.pca_trace_object(sobjs_final[indx_obj_id].TRACE_SPAT.T, + order=coeff_npoly, npca=npca, + pca_explained_var=pca_explained_var, + trace_wgt=np.fmax(sobjs_final[indx_obj_id].ech_snr, 1.0) ** 2, + debug=show_pca) + + # Trial and error shows weighting by S/N instead of S/N^2 performs better + # JXP -- Updated to now be S/N**2, i.e. inverse variance, with fitting fit + + # Perform iterative flux weighted centroiding using new PCA predictions + xinit_fweight = pca_fits[:, :, iobj].copy() + inmask_now = inmask & allmask + xfit_fweight = fit_trace(image, xinit_fweight, ncoeff, bpm=np.invert(inmask_now), + trace_bpm=np.invert(trc_inmask), fwhm=fwhm, maxdev=maxdev, + debug=show_fits)[0] + + # Perform iterative Gaussian weighted centroiding + xinit_gweight = xfit_fweight.copy() + xfit_gweight = fit_trace(image, xinit_gweight, ncoeff, bpm=np.invert(inmask_now), + trace_bpm=np.invert(trc_inmask), weighting='gaussian', fwhm=fwhm, + maxdev=maxdev, debug=show_fits)[0] + + # TODO Assign the new traces. Only assign the orders that were not orginally detected and traced. If this works + # well, we will avoid doing all of the iter_tracefits above to make the code faster. + for iord, spec in enumerate(sobjs_final[indx_obj_id]): + # JFH added the condition on ech_frac_was_fit with S/N cut on 7-7-19. + # TODO is this robust against half the order being masked? + if spec.ech_frac_was_fit & (spec.ech_snr > 1.0): + spec.TRACE_SPAT = xfit_gweight[:, iord] + spec.SPAT_PIXPOS = spec.TRACE_SPAT[specmid] + + # TODO Put in some criterion here that does not let the fractional position change too much during the iterative + # tracefitting. The problem is spurious apertures identified on one slit can be pulled over to the center of flux + # resulting in a bunch of objects landing on top of each other. + + # Set the IDs + sobjs_final[:].ECH_ORDER = order_vec[sobjs_final[:].ECH_ORDERINDX] + # for spec in sobjs_final: + # spec.ech_order = order_vec[spec.ECH_ORDERINDX] + sobjs_final.set_names() + + if show_trace: + viewer, ch = display.show_image(image * allmask) + + for spec in sobjs_trim: + color = 'red' if spec.ech_frac_was_fit else 'magenta' + ## Showing the final flux weighted centroiding from PCA predictions + display.show_trace(viewer, ch, spec.TRACE_SPAT, spec.NAME, color=color) + + for iobj in range(nobj_trim): + for iord in range(norders): + ## Showing PCA predicted locations before recomputing flux/gaussian weighted centroiding + display.show_trace(viewer, ch, pca_fits[:, iord, iobj], str(uni_frac[iobj]), color='yellow') + ## Showing the final traces from this routine + display.show_trace(viewer, ch, sobjs_final.TRACE_SPAT[iord].T, sobjs_final.NAME, color='cyan') + + # Labels for the points + text_final = [dict(type='text', args=(nspat / 2 - 40, nspec / 2, 'final trace'), + kwargs=dict(color='cyan', fontsize=20))] + + text_pca = [dict(type='text', args=(nspat / 2 - 40, nspec / 2 - 30, 'PCA fit'), + kwargs=dict(color='yellow', fontsize=20))] + + text_fit = [dict(type='text', args=(nspat / 2 - 40, nspec / 2 - 60, 'predicted'), + kwargs=dict(color='red', fontsize=20))] + + text_notfit = [dict(type='text', args=(nspat / 2 - 40, nspec / 2 - 90, 'originally found'), + kwargs=dict(color='magenta', fontsize=20))] + + canvas = viewer.canvas(ch._chname) + canvas_list = text_final + text_pca + text_fit + text_notfit + canvas.add('constructedcanvas', canvas_list) + # TODO two things need to be debugged. 1) For objects which were found and traced, i don't think we should be updating the tracing with + # the PCA. This just adds a failutre mode. 2) The PCA fit is going wild for X-shooter. Debug that. + # Vette + for sobj in sobjs_final: + if not sobj.ready_for_extraction(): + msgs.error("Bad SpecObj. Can't proceed") + + return sobjs_final \ No newline at end of file From f8afdbd31b494a1f5941effd007faad5f99fbada Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Thu, 15 Aug 2024 12:49:07 -1000 Subject: [PATCH 16/17] remove None fluxes when creating sensfunc --- pypeit/sensfunc.py | 2 +- pypeit/specobjs.py | 25 +++++++++++++++++++++---- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/pypeit/sensfunc.py b/pypeit/sensfunc.py index fdd7f4283c..9b3760385c 100644 --- a/pypeit/sensfunc.py +++ b/pypeit/sensfunc.py @@ -255,7 +255,7 @@ def __init__(self, spec1dfile, sensfile, par, par_fluxcalib=None, debug=False, # Unpack standard wave, counts, counts_ivar, counts_mask, trace_spec, trace_spat, self.meta_spec, header \ - = self.sobjs_std.unpack_object(ret_flam=False, extract_type=self.extr) + = self.sobjs_std.unpack_object(ret_flam=False, extract_type=self.extr, remove_missing=True) # Compute the blaze function # TODO Make the blaze function optional diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index deff44be53..965c414610 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -188,7 +188,7 @@ def nobj(self): """ return len(self.specobjs) - def unpack_object(self, ret_flam=False, extract_type='OPT'): + def unpack_object(self, ret_flam=False, extract_type='OPT', remove_missing=False): """ Utility function to unpack the sobjs for one object and return various numpy arrays describing the spectrum and meta @@ -198,6 +198,11 @@ def unpack_object(self, ret_flam=False, extract_type='OPT'): Args: ret_flam (:obj:`bool`, optional): If True return the FLAM, otherwise return COUNTS. + extract_type (:obj:`str`, optional): + Extraction type to use. Default is 'OPT'. + remove_missing (:obj:`bool`, optional): + If True, remove any missing data (i.e. where the flux is None). + Default is False. Returns: tuple: Returns the following where all numpy arrays @@ -217,15 +222,27 @@ def unpack_object(self, ret_flam=False, extract_type='OPT'): spec1d file """ # Prep - norddet = self.nobj flux_attr = 'FLAM' if ret_flam else 'COUNTS' flux_key = '{}_{}'.format(extract_type, flux_attr) wave_key = '{}_WAVE'.format(extract_type) - if np.any([f is None for f in getattr(self, flux_key)]): + # Check for missing data + none_flux = [f is None for f in getattr(self, flux_key)] + if np.any(none_flux): other = 'OPT' if extract_type == 'BOX' else 'BOX' - msgs.error(f"Some or all fluxes are not available for {extract_type} extraction. Try {other} extraction.") + msg = f"{extract_type} extracted flux is not available for all slits/orders. " \ + f"Consider trying the {other} extraction." + if not remove_missing: + msgs.error(msg) + else: + msg += f"{msgs.newline()}-- The missing data will be removed --" + msgs.warn(msg) + # Remove missing data + r_indx = np.where(none_flux)[0] + self.remove_sobj(r_indx) + # + norddet = self.nobj nspec = getattr(self, flux_key)[0].size # Allocate arrays and unpack spectrum wave = np.zeros((nspec, norddet)) From 4daf23206de0053c089d86b3bf89951ce8701b50 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Mon, 19 Aug 2024 11:39:59 -1000 Subject: [PATCH 17/17] added a note --- pypeit/specobjs.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/pypeit/specobjs.py b/pypeit/specobjs.py index 965c414610..4b6da7ef36 100644 --- a/pypeit/specobjs.py +++ b/pypeit/specobjs.py @@ -1108,11 +1108,12 @@ def lst_to_array(lst, mask=None): # Otherwise, we have to set the array type to object return np.array(lst, dtype=object)[_mask] - # TODO: The dtype="object" is needed for the case where one element of lst - # is not a list but None. This would prevent the error "ValueError: setting - # an array element with a sequence. The requested array has an inhomogeneous - # shape after 1 dimensions..." However, this is may introduce other issues, - # where an array of objects creates problems in other parts of the code. I - # am using this workaround. Any suggestions/ideas? + # NOTE: The dtype="object" is needed for the case where one element of lst + # is not a list but None. For example, if trying to unpack SpecObjs OPT fluxes + # and for one slit/order the OPT extraction failed (but not the BOX extraction), + # OPT_COUNTS is None for that slit/order, and lst would be something like + # [array, array, array, None, array], which makes np.array to fail and give the error + # "ValueError: setting an array element with a sequence. The requested array has an + # inhomogeneous shape after 1 dimensions..."